function [TARGET APPROX] = p2() 

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);
    amY = marginalY(approxpmf);
    amZ = marginalZ(approxpmf);

    for i=1:size(approxpmf,1),
      for j=1:size(approxpmf,2),
	for k=1:size(approxpmf,3),
	  xratio = truemargX(i)/amX(i)
	  yratio = truemargY(j)/amY(j)
	  zratio = truemargZ(k)/amZ(k)

	  nextapprox(i,j,k) = ...  
	      approxpmf(i,j,k)*(truemargX(i)/amX(i))*(truemargY(j)/amY(j))*(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,20);

end


