[R] Estimating parameters of 3 parameters lognormal distribution

Prof Brian Ripley ripley at stats.ox.ac.uk
Sat Jan 18 10:22:28 CET 2014


On 18/01/2014 09:06, Frede Aakmann Tøgersen wrote:
> It has to do with setting a to large upper bound for the gamma
> parameter. However, we are estimating the value of gamma because we
> do not know it so how can we set an upper bound for gamma??

You have starting values, so have some idea of the MLE ....

> Perhaps other optimization algorithms in optim() can be used
> (Nelder-Mead works for this x at least).

The parameters could be transformed here, and probably should.  But most 
optimization methods are going to cope with such simple constraints 
provided foo() returns NA outside the feasible region.

> Here is my findings:
>
>> N <- 100
>> mu <-  1
>> sigma <- 2
>> gamma <- 3
>>
>> set.seed(42)
>> x <- rlnorm(N, mu, sigma) + gamma
>>
>> min(x)
> [1] 3.006832
>>
>> foo <- function(x, mu, sigma, gamma)
> +                 {dlnorm(x - gamma, mu, sigma)}
>>
>> ## min(x) is not the right upper bound to use for gamma
>> fit  <- fitdistr(x, densfun = foo,
> +                 start = list(mu = 0.9, sigma = 1.9, gamma = 2.9),
> +                 lower = c(-Inf, 0, -Inf),
> +                 upper=c(Inf, Inf, min(x)))
> Error in stats::optim(x = c(45.178764955739, 3.87862565957867, 8.61957940802123,  :
>    L-BFGS-B needs finite values of 'fn'
>>
>> ## now upper bound for gamma is not to large
>> fit  <- fitdistr(x, densfun = foo,
> +                 start = list(mu = 0.9, sigma = 1.9, gamma = 2.9),
> +                 lower = c(-Inf, 0, -Inf),
> +                 upper=c(Inf, Inf, 3))
>> fit
>         mu         sigma        gamma
>    1.08694600   2.02546254   2.99393215
>   (0.20902998) (0.17736860) (0.01663611)
>>
>> ## now upper bound for gamma is not to large
>> ## still able to find correct estimates for parameters for other start values
>> fit  <- fitdistr(x, densfun = foo,
> +                 start = list(mu = 0, sigma = 1, gamma = 0),
> +                 lower = c(-Inf, 0.0001, -Inf),
> +                 upper=c(Inf, Inf, 3))
>> fit
>         mu         sigma        gamma
>    1.08671881   2.02579806   2.99394094
>   (0.20907178) (0.17745706) (0.01663591)
>>
>> ## Nelder-Mead method seems to do the job
>> fit  <- fitdistr(x, densfun = foo, method = "Nelder-Mead",
>>                 start = list(mu = 0, sigma = 1, gamma = 0))
>> fit
>         mu         sigma        gamma
>    1.08678997   2.02571374   2.99393924
>   (0.20906041) (0.17743258) (0.01663503)
>>
>> ## these two functions are defined and used in the code for fitdistr
>> ## myfn is the objective function to be optimized in optim()
>> dens <- function(parm, ...) foo(x, parm, ...)
>> myfn <- function(parm, ...) -sum(log(dens(parm, ...)))
>>
>> ## for gamma > 3. Here min(x) = 3.006832
>> ## have at least one illegal value of log-normal
>> table(dens(1, 2, min(x)) > 0)
>
> FALSE  TRUE
>      1    99
>>
>> ## resulting in
>> myfn(0.9, 1.9, min(x))
> [1] Inf
>> ## but for gamma = 3 we get finite likelihood
>> myfn(0.9, 1.9, 3.000000)
> [1] 322.4375
>> ## also finite likelihood
>> myfn(0.9, 1.9, -3.000000)
> [1] 456.7857
>>
>
> Yours sincerely / Med venlig hilsen
>
>
> Frede Aakmann Tøgersen
> Specialist, M.Sc., Ph.D.
> Plant Performance & Modeling
>
> Technology & Service Solutions
> T +45 9730 5135
> M +45 2547 6050
> frtog at vestas.com
> http://www.vestas.com
>
> Company reg. name: Vestas Wind Systems A/S
> This e-mail is subject to our e-mail disclaimer statement.
> Please refer to www.vestas.com/legal/notice
> If you have received this e-mail in error please contact the sender.
>
>
>> -----Original Message-----
>> From: Rolf Turner [mailto:r.turner at auckland.ac.nz]
>> Sent: 17. januar 2014 22:47
>> To: Vito Ricci
>> Cc: Frede Aakmann Tøgersen; Göran Broström; r-help at stat.math.ethz.ch
>> Subject: Re: [R] Estimating parameters of 3 parameters lognormal
>> distribution
>>
>>
>> Can you please tell us (me!) how you chose starting values?
>>
>> Out of curiosity I tried the following:
>>
>> set.seed(42)
>> x <- rlnorm(100,1,2) + 3
>> require(MASS)
>> strt <- list(mu=1,sigma=2,gamma=3)
>> fit  <- fitdistr(x,densfun=function(x,mu,sigma,gamma)
>>                                     {dlnorm(x-gamma,mu,sigma)
>>                             },
>>                  start=strt,lower=c(0,0,-Inf),
>>                  upper=c(Inf,Inf,min(x)))
>>
>> and it ran just fine and gave sensible answers.  But when I took:
>>
>> strt <- list(mu=0.9,sigma=1.9,gamma=2.9)
>>
>> (not very different from the previous starting values)
>>
>> I got an error:
>>
>>> Error in stats::optim(x = c(45.178764955739, 3.87862565957867,
>> 8.61957940802123,  :
>>>    L-BFGS-B needs finite values of 'fn'
>>
>>
>> On 18/01/14 01:00, Vito Ricci wrote:
>>> OK. It runs fine! Many thanks Frede.
>>> Regards.
>>> Vito
>>>
>>>
>>>
>>> Se non ora, quando?
>>> Se non qui, dove?
>>> Se non tu, chi?
>>>
>>>
>>>
>>> Il Venerdì 17 Gennaio 2014 12:38, Frede Aakmann Tøgersen
>> <frtog at vestas.com> ha scritto:
>>>
>>>
>>>> In package MASS there is the fitdistr function using maximum likelihood
>> estimation to infer on the parameters of distributions based on observed
>> data.
>>>>
>>>>
>>>> One of the arguments of fitdistr () allows you to specify the probability
>> density function.
>>>>
>>>>
>>>> You only know how the parametrization of your three parameters
>> lognormal distribution is defined since you really haven't told us much.
>>>>
>>>>
>>>> Please have a look at fitdistr() and tell us if fit your needs.
>
> ______________________________________________
> 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




More information about the R-help mailing list