% ITEST file itest.m  Tests P for independence
% Version of 5/9/95
% Tests for independence the matrix of joint 
% probabilities for a simple pair {X,Y}
pt = input('Enter matrix of joint probabilities  ');
disp(' ')
px = sum(pt);                  % Marginal probabilities for X
py = sum(pt');                 % Marginal probabilities for Y (reversed)
[a,b] = meshgrid(px,py); 
PT = a.*b;                     % Joint independent probabilities
D  = abs(pt - PT) > 1e-9;      % Threshold set above roundoff
if total(D) > 0
  disp('The pair {X,Y} is NOT independent')
  disp('To see where the product rule fails, call for D')
else
  disp('The pair {X,Y} is independent')
end