[R] Likelihood returning inf values to optim(L-BFGS-B) other options?

[Prof Brian Ripley]

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