[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