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