R help - Feb 2010 - using mle2 for multinomial model optimization

Hi there
I'm trying to find the mle fo a multinomial model ->*L(N,h,S?x)*. There 
is only *N* I want to estimate, which is used in the number of successes 
for the last cell probability. These successes are given by: 
p^(N-x1-x2-...xi)
 All the other parameters (i.e. h and S) I know from somewhere else. 
Here is what I've tried to do so far for a imaginary data set:
#######################################################

cohort<-c(50,54,50)     #imaginary harvest numbers of a cohort harvested 
over 3 years

l<-length(cohort)+1    #nr of cell probabilities
h<-c(0.2,0.3,1)         #harvest rates for the 3 diferent years
S<-c(0.9,0.8)         #survival rates

mfun <- function(d) {
  S<-S            #survival rate
  h<-h            #harvest rate
  l<-length(d)
    p<-rep(0,l)         #Cell probabilities
    p[1]<-h[1]
  
    prod<-(1-h[1])*S[1]
    for (i in 2:(l-1)){
        p[i]<-prod*h[i]
        prod<-prod*(1-h[i])*S[i]           

    }
   p[l]<-1-sum(p[-l])    #last cell probability   
  -dmultinom(exp(d),prob=p,log=TRUE)    #exp(d)->backtransform the estimates
}

c<-c(cohort,100)            #untransformed initial values for the 
optimization,->100=N-x1-x2-x3
nvec<-c(rep("x",l-1),"N")     #names for the initial
vector
nrs<-c(1:(l-1),1)            #nrs for the initial vector
svec = log(c)             #transformation of the data to avoid 
constraints (x>0)
names(svec) <- paste(nvec,nrs,sep="")
parnames(mfun) <- names(svec)

m1 = mle2(mfun,start=svec,
     vecpar=TRUE,fixed=svec[1:(l-1)]) #only estimate 
"N-x1-x2-x3",everything else is fixed

    coef(m1)
############################################################################

The function "mfun" is working, however the mle2 won't work. How
can I
fix this?

Thanks so much for your help

Beni

Use "Nelder-Mead" as the method in optim.  This will give you the
correct solution.

m1 <- mle2(mfun, start=svec, vecpar=TRUE, fixed=svec[1:(l-1)],
method="Nelder")

Hope this is helpful,
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: Benedikt Gehr <benedikt.gehr at ieu.uzh.ch>
Date: Friday, February 12, 2010 8:26 am
Subject: [R] using mle2 for multinomial model optimization
To: r-help at r-project.org

> Hi there
>  I'm trying to find the mle fo a multinomial model ->*L(N,h,S?x)*. 
> There 
>  is only *N* I want to estimate, which is used in the number of 
> successes 
>  for the last cell probability. These successes are given by: 
>  p^(N-x1-x2-...xi)
>   All the other parameters (i.e. h and S) I know from somewhere else. 
> 
>  Here is what I've tried to do so far for a imaginary data set:
>  #######################################################
>  
>  cohort<-c(50,54,50)     #imaginary harvest numbers of a cohort 
> harvested 
>  over 3 years
>  
>  l<-length(cohort)+1    #nr of cell probabilities
>  h<-c(0.2,0.3,1)         #harvest rates for the 3 diferent years
>  S<-c(0.9,0.8)         #survival rates
>  
>  mfun <- function(d) {
>    S<-S            #survival rate
>    h<-h            #harvest rate
>    l<-length(d)
>      p<-rep(0,l)         #Cell probabilities
>      p[1]<-h[1]
>    
>      prod<-(1-h[1])*S[1]
>      for (i in 2:(l-1)){
>          p[i]<-prod*h[i]
>          prod<-prod*(1-h[i])*S[i]           
>  
>      }
>     p[l]<-1-sum(p[-l])    #last cell probability   
>    -dmultinom(exp(d),prob=p,log=TRUE)    #exp(d)->backtransform the
estimates
>  }
>  
>  c<-c(cohort,100)            #untransformed initial values for the 
>  optimization,->100=N-x1-x2-x3
>  nvec<-c(rep("x",l-1),"N")     #names for the
initial vector
>  nrs<-c(1:(l-1),1)            #nrs for the initial vector
>  svec = log(c)             #transformation of the data to avoid 
>  constraints (x>0)
>  names(svec) <- paste(nvec,nrs,sep="")
>  parnames(mfun) <- names(svec)
>  
>  m1 = mle2(mfun,start=svec,
>       vecpar=TRUE,fixed=svec[1:(l-1)]) #only estimate 
>  "N-x1-x2-x3",everything else is fixed
>  
>      coef(m1)
> 
############################################################################
>  
>  The function "mfun" is working, however the mle2 won't work.
How can
> I 
>  fix this?
>  
>  Thanks so much for your help
>  
>  Beni
>  
>  ______________________________________________
>  R-help at r-project.org mailing list
>  
>  PLEASE do read the posting guide 
>  and provide commented, minimal, self-contained, reproducible code.