[R] MLE Function

[Peter Dalgaard]

Terence Broderick wrote:
> I am just trying to teach myself how to use the mle function in R because it is much better than what is provided in MATLAB. I am following tutorial material from the internet, however, it gives the following errors, does anybody know what is happening to cause such errors, or does anybody know any better tutorial material on this particular subject.
>   
>   
>> x.gam<-rgamma(200,rate=0.5,shape=3.5)
>> x<-x.gam
>> library(stats4)
>> ll<-function(lambda,alfa){n<-200;x<-x.gam -n*alfa*log(lambda)+n*log(gamma(alfa))-9alfa-1)*sum(log(x))+lambda*sum(x)}
>>     
> Error: syntax error, unexpected SYMBOL, expecting '\n' or ';' or '}' in "ll<-function(lambda,alfa){n<-200;x<-x.gam -n*alfa*log(lambda)+n*log(gamma(alfa))-9alfa"
>   
>> ll<-function(lambda,alfa){n<-200;x<-x.gam -n*alfa*log(lambda)+n*log(gamma(alfa))-(alfa-1)*sum(log(x))+lambda*sum(x)}
>> est<-mle(minuslog=ll,start=list(lambda=2,alfa=1))
>>     
> Error in optim(start, f, method = method, hessian = TRUE, ...) : 
>         objective function in optim evaluates to length 200 not 1
>
>    
>   
Er, not what I get. Did your version have that linefeed after x <- x.gam 
? If not, then you'll get your negative log-likelihood added to x.gam 
and the resulting "likelihood" becomes a vector of length 200 instead of 
a scalar.

In general, the first piece of advice for mle() is to check that the 
likelihood function really is what it should be. Otherwise there is no 
telling what the result might mean...

Secondly, watch out for parameter constraints. With your function, it 
very easily happens that "alfa" tries to go negative in which case the 
gamma function in the likelihood will do crazy things.
A common trick in such cases is to reparametrize by log-parameters, i.e.

ll <- function(lambda,alfa){n<-200; x<-x.gam
-n*alfa*log(lambda)+n*lgamma(alfa)-(alfa-1)*sum(log(x))+lambda*sum(x)}

ll2 <- function(llam, lalf) ll(exp(llam),exp(lalf))
est <- mle(minuslog=ll2,start=list(llam=log(2),lalf=log(1)))

par(mfrow=c(2,1))
plot(profile(est))

Notice, incidentally, the use of lgamma rather than log(gamma(.)), which 
is prone to overflow.

In fact, you could also write this likelihood directly  as

-sum(dgamma(x, rate=lambda, shape=alfa, log=T))


>    
>
>
> audaces fortuna iuvat
>        
> ---------------------------------
>
> 	[[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at stat.math.ethz.ch 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.
>   


-- 
   O__  ---- Peter Dalgaard             Øster Farimagsgade 5, Entr.B
  c/ /'_ --- Dept. of Biostatistics     PO Box 2099, 1014 Cph. K
 (*) \(*) -- University of Copenhagen   Denmark          Ph:  (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)                  FAX: (+45) 35327907