%% EM example, simple case
clear

%% Data
NA = 725; 
NB = 258;
NAB = 72;
NO = 1073;

N=2128;

if (NA+NB+NAB+NO~=N), warning('wierdness'); end
    
%% EM

%Initialise
p=1/3; q=1/3; r=1-p-q;

%EM loop
notdone=true;
while (notdone),
    %E
    NAA=NA*p^2/(p^2+2*p*r);
    NAO=NA-NAA;
    NBB=NB*q^2/(q^2+2*q*r);
    NBO=NB-NBB;

    %M
    pp=(2*NAA+NAO+NAB)/(2*N);
    qq=(2*NBB+NBO+NAB)/(2*N);
    rr=1-pp-qq;

    notdone=any(abs([(p-pp)/p,(q-qq)/q])>eps);
    p=pp; q=qq; r=rr;
    display(sprintf(...
      '%f %f %f %f %f %f %f',...
      NAA,NAO,NBB,NBO,p,q,r...
    ));
end
display(sprintf('Converged to MP\n'));