[R] Help: Maximum likelihood estimation

[roach]

I'm not quite familiar with E-M algorithm, but I think what I did was the
first step of the iteration. 
The method used in the original article is as follow:
http://r.789695.n4.nabble.com/file/n3008241/untitled.png 
I gave lamda an initial value, and maximized the likelihood function.
This is the complete chunk of my code after using "alabama" package. 
The first iteration had no problem, but after a few iterations, I again got
warnings and the result is not good.
 Is it possible that it's because of some computational problems? because
there are too many log and exp in the functions? Or is there anything I
missed?


library("alabama")
# n=number of observation
w=seq(0.05,0.9,length.out=n)
# iteration
repeat{
lamda=mean(w)
## -log likelihood function
log.L=function(parm){
alpha0=parm[1]
alpha1=parm[2]
alpha2=parm[3]
beta0=parm[4]
beta1=parm[5]
beta2=parm[6]
beta3=parm[7]
# here sigma is actually sigma inverse
sigma11=parm[8]
sigma12=parm[9]
sigma21=parm[10]
sigma22=parm[11]

u1=-alpha0-alpha1*logp-alpha2*lakes+logq
u21=-beta0-beta1*logq-beta2*s-beta3+logp
u22=-beta0-beta1*logq-beta2*s+logp
expon1=u1^2*sigma11+u1*u21*sigma12+u1*u21*sigma21+u21^2*sigma22
expon2=u1^2*sigma11+u1*u22*sigma12+u1*u22*sigma21+u22^2*sigma22
const=-log(2*pi)+.5*log(sigma11*sigma22-sigma12*sigma21)+log(abs(1-alpha1*beta1))
logf=const+log(lamda*exp(-0.5*expon1)+(1-lamda)*exp(-0.5*expon2))
log.L=-sum(logf)
return(log.L)
}

## estimate with nonlinear constraint
hin=function(parm){
h=rep(NA,1)
h[1]=parm[8]*parm[11]-parm[9]*parm[10]
h[2]=parm[8]
h[3]=parm[11]
h
}

heq=function(parm){
h=rep(NA,1)
h[1]=parm[9]-parm[10]
h
}
max.like=constrOptim.nl(par=c(-0.5,-0.5,-0.5,-0.5,0.02,-0.02,0.02,1.9,-1.1,-1.1,1.9),fn=log.L,
hin=hin,heq=heq)
max.like$par


######
parm=max.like$par
alpha0=parm[1]
alpha1=parm[2]
alpha2=parm[3]
beta0=parm[4]
beta1=parm[5]
beta2=parm[6]
beta3=parm[7]
sigma11=parm[8]
sigma12=parm[9]
sigma21=parm[10]
sigma22=parm[11]
u1=-alpha0-alpha1*logp-alpha2*lakes+logq
u21=-beta0-beta1*logq-beta2*s-beta3+logp
u22=-beta0-beta1*logq-beta2*s+logp
expon1=u1^2*sigma11+u1*u21*sigma12+u1*u21*sigma21+u21^2*sigma22
expon2=u1^2*sigma11+u1*u22*sigma12+u1*u22*sigma21+u22^2*sigma22

h1_log=(-log(2*pi)+0.5*log(sigma11*sigma22-sigma12*sigma21))+(log(abs(1-alpha1*beta1))-0.5*expon1)
h2_log=(-log(2*pi)+0.5*log(sigma11*sigma22-sigma12*sigma21))+(log(abs(1-alpha1*beta1))-0.5*expon2)
w1=w

w=1/(1+(1-lamda)/lamda*exp(h2_log-h1_log))
if(cor(w,w1)>0.999) break
}


-- 
View this message in context: http://r.789695.n4.nabble.com/Help-Maximum-likelihood-estimation-tp3006584p3008241.html
Sent from the R help mailing list archive at Nabble.com.