[R] Likelihood returning inf values to optim(L-BFGS-B) other options?
Prof Brian Ripley
ripley at stats.ox.ac.uk
Sat Apr 7 09:07:23 CEST 2007
On Sat, 7 Apr 2007, Michael Jungbluth wrote:
> Thank you very much for your postings. Rewriting the likelihood with lgamma
> helps a lot and the mistake with "fnscale" was quite stupid (Sorry for
> that!).
>
> The model is working for most of the parameter sets but I am still facing
> some inf-returns on my (with lgamma updated and negative) loglikelihood if I
> am putting in some extreme parameter values (e.g. u-shaped beta densities
> for the simulation data which generates x,t_x and T). Actually these are the
> ones which are really interesting for my project. So, is there another
> similar optimization algorithm which can deal with inf-returns?
Please read all the hints I gave you. If you transform the parameter
space, you can use optim(method="BFGS"). You could also try nlminb, but I
don't find that anything like as reliable.
>
> Thanks a lot!
>
> Best regards,
> Michael
>
> Zitat von Prof Brian Ripley <ripley at stats.ox.ac.uk>:
>
>> On Thu, 5 Apr 2007, Ravi Varadhan wrote:
>>
>>> In your code, the variables x (which I assume is the observed data),
>>> Tvec,
>>> and flag are not passed to the function as arguments. This could be a
>>> potential problem.
>>
>> I think scoping will probably find them.
>>
>>> Another problem could be that you have to use "negative"
>>> log-likelihood function as input to optim, since by default it
>>> "minimizes"
>>> the function, whereas you are interested in finding the argmax of
>>> log-likelihood. So, in your function you should return (-ll) instead of
>>> ll.
>>
>> OR set fnscale. This is the most serious problem.
>>
>>> If the above strategies don't work, I would try different initial values
>>> (it
>>> would be best if you have a data-driven strategy for picking a starting
>>> value) and different optimization methods (e.g. conjugate gradient with
>>> "Polak-Ribiere" steplength option, Nelder-Mead, etc.).
>>
>> It looks to me as if the calculations are very vulnerable to
>> overflow/underflow, as they use gamma and not lgamma. They could be
>> rearranged to be much stabler by computing the sum of logs for each
>> sub-expression.
>>
>> There were over 50 warnings, which we were not shown. They probably
>> explained the problem.
>>
>> Beyond that, the feasible region seems to be the interior of the
>> positive orthant, in which case transforming the parameters (e.g.
>> working with their logs) would be a good idea.
>>
>> Finally, always supply analytical gradients when you can (as would be
>> easy here).
>>
>>
>>> -----Original Message-----
>>> From: r-help-bounces at stat.math.ethz.ch
>>> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of
>>> r at micom-solutions.de
>>> Sent: Thursday, April 05, 2007 6:12 AM
>>> To: r-help at stat.math.ethz.ch
>>> Subject: [R] Likelihood returning inf values to optim(L-BFGS-B) other
>>> options?
>>>
>>> Dear R-help list,
>>>
>>> I am working on an optimization with R by evaluating a likelihood
>>> function that contains lots of Gamma calculations (BGNBD: Hardie Fader
>>> Lee 2005 Management Science). Since I am forced to implement lower
>>> bounds for the four parameters included in the model, I chose the
>>> optim() function mith L-BFGS-B as method. But the likelihood often
>>> returns inf-values which L-BFGS-B can't deal with.
>>>
>>> Are there any other options to implement an optimization algorithm
>>> with R accounting for lower bounds and a four parameter-space?
>>>
>>> Here is the error message I receive (german):
>>> --
>>>>
>>> out=optim(c(.1,.1,.1,.1),fn,method="L-BFGS-B",lower=c(.0001,.0001,.0001,.000
>>> 1,.0001))
>>> Fehler in optim(c(0.1, 0.1, 0.1, 0.1), fn, method = "L-BFGS-B", lower
>>> = c(1e-04, :
>>> L-BFGS-B benötigt endliche Werte von 'fn'
>>> Zusätzlich: Es gab 50 oder mehr Warnungen (Anzeige der ersten 50 mit
>>> warnings())
>>> --
>>> And this is the likelihood function:
>>> --
>>> fn<-function(p) {
>>> A1=(gamma(p[1]+x)*p[2]^p[1])/(gamma(p[1]))
>>> A2=(gamma(p[3]+p[4])*gamma(p[4]+x))/(gamma(p[4])*gamma(p[3]+p[4]+x))
>>> A3=(1/(p[2]+Tvec))^(p[1]+x)
>>> A4=(p[3]/(p[4]+x-1))*((1/(p[2]+t_x))^(p[1]+x))
>>> ll=sum(log(A1*A2*(A3+flag*A4)))
>>> return(ll)
>>> }
>>>
>>> Thank you very much for your help in advance!
>>>
>>> Best regards,
>>>
>>> Michael
>>>
>>> ______________________________________________
>>> R-help at 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.
>>>
>>> ______________________________________________
>>> R-help at 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.
>>>
>>
>> --
>> 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
>
>
>
>
--
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