[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