%Rao's EM Example
function psi=emexample(y1, y2, y3, y4, tol, start) 
% function psi=emexample(y1, y2, y3, y4, tol, start) 
%
% Use: emexample(125, 18, 20, 34, 10^-6, 0.5)
%
psi_last = 0; 
n = y1 + y2 + y3 + y4; 
psi_current = start; 
psi = psi_current; 
%------------
while (abs(psi_last-psi)>tol )
   [y11, y12] = estep(psi_current, y1);
   psi = mstep(y12, y11, y4, n);
   psi_last = psi_current;
   psi_current = psi;
end %while

%-----------
function psi_new = mstep(y12, y11, y4, n)
 psi_new = (y12+y4)/(n-y11);
%-----------
function [y11, y12] = estep(psi_current, y1)
 y11 = y1 * (1/2)/(1/2 + psi_current/4);
 y12 = y1 - y11;
%------------

