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: rvarad...@jhmi.edu
----- Original Message -----
From: Benedikt Gehr <benedikt.g...@ieu.uzh.ch>
Date: Friday, February 12, 2010 8:26 am
Subject: [R] using mle2 for multinomial model optimization
To: r-help@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@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.
