Can you provide a reproducible code?
Ravi.
____________________________________________________________________
Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University
Ph. (410) 502-2619
email: rvarad...@jhmi.edu
----- Original Message -----
From: roach <roachy...@gmail.com>
Date: Saturday, October 23, 2010 4:41 am
Subject: Re: [R] Help: Maximum likelihood estimation
To: r-help@r-project.org
> 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:
>
> 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:
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> R-help@r-project.org mailing list
>
> PLEASE do read the posting guide
> and provide commented, minimal, self-contained, reproducible code.
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
