compute maximum likelihood estimator for a multinomial function
Benedikt Gehr wrote:
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
______________________________________________ R-help at 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.
There are a few issues here.
The main one is that the parameters of a multinomial
have particular constraints (0<=p_i, sum(p_i)=1) which
are somewhat challenging to implement -- the 0<=p_i<=1
part is easy, but the sum constraint is hard (even
if you say that you're only going to specify the first N-1
parameters and let the last one be 1-sum_{i=1}^{N-1} p_i,
that still doesn't prevent the sum of the other parameters
from being >1). One approach is to use a transformation
like the additive log-ratio transform (described elsewhere
in the book).
By the way, you can show analytically that the MLEs
for the multinomial parameters are the observed frequencies ...
I'm not sure about confidence intervals off the top of
my head, but perhaps someone else will chime in.
The other issue is that mle() and mle2() require the
parameters to be specified separately (i.e. fun(p1,p2,p3,p4)
rather than fun(pvec)), although mle2() allows you to work
around this constraint.
View this message in context: http://old.nabble.com/compute-maximum-likelihood-estimator-for-a-multinomial-function-tp26195447p26196094.html Sent from the R help mailing list archive at Nabble.com.