R help - Nov 2009 - compute maximum likelihood estimator for a multinomial function

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

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.