[R] Maximum likelihood estimation

toh tohbin at yahoo.com
Fri Sep 5 03:15:06 CEST 2008


Yes I'm trying to optimize the parameters a, b, p and lambda where a > 0, b >
0 and 0 < p < 1. I attached the error message that I got when I run mle. 


> t <- c(1:90)
> y <-
> c(5,10,15,20,26,34,36,43,47,49,80,84,108,157,171,183,191,200,204,211,217,226,230,
+
234,236,240,243,252,254,259,263,264,268,271,277,290,309,324,331,346,367,375,381,
+
401,411,414,417,425,430,431,433,435,437,444,446,446,448,451,453,460,463,463,464,
+
464,465,465,465,466,467,467,467,468,469,469,469,469,470,472,472,473,473,473,473,
+ 473,473,473,475,475,475,475)
> 
> library(stats4)
> fr <- function(a, b, p, lambda){
+ l <- 0.5*(lambda + b*p + (1-p)*(lambda-b))
+ l^2 > lambda*b*p
+ delta <- sqrt(abs(l^2 - b*p*lambda))
+ mt <- a/p*(1-exp(-l*t)*cosh(delta*t)-(l-b*p)*exp(-t)*sinh(delta*t)/delta)
+ logl <- sum(diff(y)*log(diff(mt))-diff(mt)-lfactorial(diff(y)))
+ return(-logl)
+ }
> 
> mle(start=list(a=12,b=0.01,p=0.5,lambda=0.01),fr, method="L-BFGS-B", 
+ lower = c(0.002, 0.002, 0.002, 0.002), upper = c(Inf, Inf, 0.999,
Inf),control=list(fnscale=-1))
Error in optim(start, f, method = method, hessian = TRUE, ...) : 
  non-finite finite-difference value [3]












Prof Brian Ripley wrote:
> 
>>From ?optim
> 
>        fn: A function to be minimized (or maximized), with first
>            argument the vector of parameters over which minimization is
>            to take place.  It should return a scalar result.
> 
> I think you intended to optimize over c(a,b,p, lambda), so you need to 
> specify them as a single vector.
> 
> You may be making life unnecessarily hard for yourself: see function mle() 
> in package stats4.
> 
> Showing your code without a verbal description of what you are doing nor 
> the error message you got is less helpful than we need.
> 
> On Wed, 3 Sep 2008, toh wrote:
> 
>>
>> Hi R-experts,
>> I'm new to R in mle. I tried to do the following but just couldn't get it
>> right. Hope someone can point out the mistakes. thanks a lot.
>>
>> t <- c(1:90)
>> y <-
>> c(5,10,15,20,26,34,36,43,47,49,80,84,108,157,171,183,191,200,204,211,217,226,230,
>>
>> 234,236,240,243,252,254,259,263,264,268,271,277,290,309,324,331,346,367,375,381,
>>
>> 401,411,414,417,425,430,431,433,435,437,444,446,446,448,451,453,460,463,463,464,
>>
>> 464,465,465,465,466,467,467,467,468,469,469,469,469,470,472,472,473,473,473,473,
>> 	473,473,473,475,475,475,475)
>> fr <- function(a, b, p, lambda){
>> l <- 0.5*(lambda + b*p + (1-p)*(lambda-b))
>> l^2 > lambda*b*p
>> delta <- sqrt(abs(l^2 - b*p*lambda))
>> mt <- a/p*(1-exp(-l*t)*cosh(delta*t)-(l-b*p)*exp(-t)*sinh(delta*t)/delta)
>> logl <- sum(diff(y)*log(diff(mt))-diff(mt)-lfactorial(diff(y)))
>> return(-logl)
>> }
>> optim(c(15,0.01,0.5,0.01),fr, method="L-BFGS-B",
>> lower = c(0.002, 0.002, 0.002, 0.002), upper = c(Inf, Inf, 0.999,
>> Inf),control=list(fnscale=-1))
>>
>> -- 
>> View this message in context:
>> http://www.nabble.com/Maximum-likelihood-estimation-tp19304249p19304249.html
>> Sent from the R help mailing list archive at Nabble.com.
>>
>> ______________________________________________
>> 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.
>>
> 
> -- 
> Brian D. Ripley,                  ripley at stats.ox.ac.uk
> Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
> University of Oxford,             Tel:  +44 1865 272861 (self)
> 1 South Parks Road,                     +44 1865 272866 (PA)
> Oxford OX1 3TG, UK                Fax:  +44 1865 272595
> 
> ______________________________________________
> 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.
> 
> 

-- 
View this message in context: http://www.nabble.com/Maximum-likelihood-estimation-tp19304249p19323140.html
Sent from the R help mailing list archive at Nabble.com.



More information about the R-help mailing list