[R] Help: Maximum likelihood estimation

[Ravi Varadhan]

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: rvaradhan at jhmi.edu


----- Original Message -----
From: roach <roachyang at gmail.com>
Date: Saturday, October 23, 2010 4:41 am
Subject: Re: [R] Help: Maximum likelihood estimation
To: r-help at 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 at r-project.org mailing list
>  
>  PLEASE do read the posting guide 
>  and provide commented, minimal, self-contained, reproducible code.