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@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
~~~~~~~~~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907
______________________________________________
R-help@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.
