[R] compute maximum likelihood estimator for a multinomial function
    Benedikt Gehr 
    Benedikt.Gehr at access.uzh.ch
       
    Wed Nov  4 13:28:50 CET 2009
    
    
  
Hi there
I am trying to learn how to compute mle in R for a multinomial negative 
log likelihood function.
I am using for this the book by B. Bolker "Ecological models and data in 
R", chapter 6: "Likelihood an all that". But he has no example for 
multinomial functions.
What I did is the following:
I first defined a function for the negative log likelihood:
#######################
log.multi<-function(d,p){
-sum(dmultinom(d,prob=p,log=T))
}
##then defined the parameters
d<-data ## (in my case counts of animals found dead in different 
age-classes 1:19)
p<-d/sum(d) ##or some values close to the analytical solution n(i)/N for 
the  probabilities
##then I used the optim function like this:
opt1<-optim(par=p,fn=log.multi,d=d)
#########################
This works perfectly when I use some simple imaginary numbers instead of 
my real data; e.g. d<-c(3,6,9)
However for my real data this doesnt work because the differences 
between the different counts seem to be to large.
here are my data:
###########################
d<-count
 > d
 [1]  7 15  9  3  6 13 11 17 15 36 38 52 49 53 33 20  6  2  1
###########################
I get the error message:  "probabilities can't neither be negative nor 
all 0"
So I wanted to try the "mle" or the "mle2" functions of the packages 
"stats4" and "bbmle" respectively. But I just can't get them running.
I tried it like this:
###########################
library(bbmle)
log.multi<-function(d,p.est){
-sum(dmultinom(x=d,prob=p.est,log=T))
}
d<-c(3,5,7)
p.est<-d/sum(d)
m1<-mle2(minuslogl=log.multi,start=list(p.est=p.est),data=list(N=sum(d),d=d))
###########################
I get the error:
"Error in dmultinom(x = d, prob = p.est, log = T) :
  x[] and prob[] must be equal length vectors."
And no matter how I tried (and I tried for a looooong time!!) I cant't 
figure out what the problem is. I'm sure it a very silly mistake but I 
don't know what it is.
I think it must have something to do with how I define the parameters.
Could anyone help me with this?
Thank you so much in advance
Benedikt
    
    
More information about the R-help
mailing list