[R] Example on Maximum likelihood estimation using R
VincentLee
vincentlee_sg at yahoo.com
Fri Sep 5 16:10:46 CEST 2008
Hi all,
I am very new to R too, but I read that R is powerful.
May I know given a set of data, are there any simple examples on using mle()
to estimate parameters of a lognormal and weibull distribution ?
Hope to hear from you soon.
Thank you
Vincent
Ravi Varadhan wrote:
>
> Hi,
>
>
> You should re-write your function `fr' as follows:
>
>
> fr <- function(par){
> a <- par[1]
> b <- par[2]
> p <- par[3]
> lambda <- par[4]
>
> 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)
> }
>
> However, I don't understand what the following fragment is doing in `fr':
>
> l^2 > lambda*b*p
>
> Can you clarfy that?
>
> Ravi.
>
> ----------------------------------------------------------------------------
> -------
>
> Ravi Varadhan, Ph.D.
>
> Assistant Professor, The Center on Aging and Health
>
> Division of Geriatric Medicine and Gerontology
>
> Johns Hopkins University
>
> Ph: (410) 502-2619
>
> Fax: (410) 614-9625
>
> Email: rvaradhan at jhmi.edu
>
> Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
>
>
>
> ----------------------------------------------------------------------------
> --------
>
>
> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
> On
> Behalf Of toh
> Sent: Thursday, September 04, 2008 9:15 PM
> To: r-help at r-project.org
> Subject: Re: [R] Maximum likelihood estimation
>
>
> 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,2
>> 17,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,3
>>> 67,375,381,
>>>
>>> 401,411,414,417,425,430,431,433,435,437,444,446,446,448,451,453,460,4
>>> 63,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-tp19304249p193042
>>> 49.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.
>
> ______________________________________________
> 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.
>
> ______________________________________________
> 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-tp19304249p19332104.html
Sent from the R help mailing list archive at Nabble.com.
More information about the R-help
mailing list