Matlab Implementation

Here is the little matlab functions that I used to check the computations:

function psi=emalgo(y1,y2,y3,y4,tol,start)
%EM algorithm for Rao's multinomial model
%
%y1,y2,y_3,y_4 are the observed frequencies
%

%Initializations

psilast=0;
n=y1+y2+y3+y4;
psicurrent=start;
psi=psicurrent;
while (abs(psilast-psi)>tol )
   [y11, y12]=estep(psicurrent,y1);
   psi=mstep(y12,y11,y4,n);
   psilast=psicurrent;
   psicurrent=psi
end %while

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function psinew=mstep(y12,y11,y4,n)
psinew=(y12+y4)/(n-y11);
%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function [y11, y12]=estep(psicurrent,y1)
 y11=(.5*y1)/(.5+psicurrent/4);
 y12=y1-y11;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Results:

 tol=10^(-7);
 start=.05;    
emalgo(y1,y2,y3,y4,tol,start)
psicurrent =
   0.04268630916291
psicurrent =
   0.03914124282012
psicurrent =
   0.03740420858837
psicurrent =
   0.03654857635551
psicurrent =
   0.03612601017717
psicurrent =
   0.03591705195885
psicurrent =
   0.03581365694687
psicurrent =
   0.03576247981151
psicurrent =
   0.03573714487659
psicurrent =
   0.03572460200443
psicurrent =
   0.03571839201694
psicurrent =
   0.03571531738852
psicurrent =
   0.03571379509429
psicurrent =
   0.03571304138032
psicurrent =
   0.03571266820275
psicurrent =
   0.03571248343550
psicurrent =
   0.03571239195371
ans =
   0.03571239195371