function psi=emexample(y1, y2, y3, y4, tol, start)
% %%Fisher's Example
% 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 = (.5*y1)/(.5+psi_current/4);
 y12 = y1 - y11;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% >> format long
% >> emexample(125, 18, 20, 34, 10^-6, 0.5)
% psi_current = 0.60824742268041
% psi_current = 0.62432105036927
% psi_current = 0.62648887907967
% psi_current = 0.62677732234731
% psi_current = 0.62681563211004
% psi_current = 0.62682071901931
% psi_current = 0.62682139445598
% 
% ans = 0.62682139445598
% 
% >> ans - (15+sqrt(53809))/394
% 
% ans = -1.034149983425436e-007