[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