function mi = mutual_information(pjoint)
% MUTUAL_INFORMATION Compute the mutual information between X and Y.
% mi = mutual_information(pjoint)
%
% pjoint(i,j) = Pr(X=i, Y=j)

pj = sum(pjoint,2);
pV = sum(pjoint,1);
pprod = repmat(pj, 1, length(pV)) .* repmat(pV, length(pj), 1);
%pprod = pprod + (pprod==0);
mi = sum(sum(pjoint .* log2( (pjoint+eps) ./ (pprod+eps) )));
