function [TARGET APPROX] = p22() 

function pmf = makeRandPmf()
  pmf = rand([2 2 2]);
  pmf = pmf / sum(sum(sum(pmf)));
end;

function marg = marginalX(jointpmf)
   marg = sum(sum(jointpmf, 2),3);
end

function marg = marginalY(jointpmf)
   marg = shiftdim(sum(sum(jointpmf,1),3));
end

function marg = marginalZ(jointpmf)
   marg = shiftdim(sum(sum(jointpmf,1),2),2);
end;

function nextapprox = ipfIter(truemargX, truemargY, truemargZ, ...
    approxpmf)
   nextapprox = [];

    amX = marginalX(approxpmf);
    for i=1:size(approxpmf,1),
      for j=1:size(approxpmf,2),
	for k=1:size(approxpmf,3),
	  nextapprox(i,j,k) = ...  
	      approxpmf(i,j,k)*(truemargX(i)/amX(i)); 
	end;
      end;
    end;

    approxpmf = nextapprox;
    amY = marginalY(approxpmf);
    
    for i=1:size(approxpmf,1),
      for j=1:size(approxpmf,2),
	for k=1:size(approxpmf,3),
	  nextapprox(i,j,k) = ...  
	      approxpmf(i,j,k)*(truemargY(j)/amY(j)); 
	end;
      end;
    end;
    
    approxpmf = nextapprox;
    amZ = marginalZ(approxpmf);
    
    for i=1:size(approxpmf,1),
      for j=1:size(approxpmf,2),
	for k=1:size(approxpmf,3),
	  nextapprox(i,j,k) = ...  
	      approxpmf(i,j,k)*(truemargZ(k)/amZ(k)); 
	end;
      end;
    end;   
    

end;

function approximation = doipf(truepmf, iterations)
   curapprox = makeRandPmf;
   for i=1:iterations,
     truepmf;
     curapprox = ipfIter(marginalX(truepmf), ...
	 marginalY(truepmf), marginalZ(truepmf), curapprox);
   end;
   approximation = curapprox;
end;

%% first come up with a random 5x5x5 distribution 
TARGET = makeRandPmf
MARGINAL_X = marginalX(TARGET)
MARGINAL_Y = marginalY(TARGET)
MARGINAL_Z = marginalZ(TARGET)

APPROX = doipf(TARGET,120);

MARGINAL_X_DIFF = MARGINAL_X - marginalX(APPROX)
MARGINAL_Y_DIFF = MARGINAL_Y - marginalY(APPROX)
MARGINAL_Z_DIFF = MARGINAL_Z - marginalZ(APPROX)


end


