From: Prof Brian Ripley <ripley_at_stats.ox.ac.uk>

Date: Sat 07 Apr 2007 - 07:07:23 GMT

Date: Sat 07 Apr 2007 - 07:07:23 GMT

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@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@stat.math.ethz.ch
**>>> [mailto:r-help-bounces@stat.math.ethz.ch] On Behalf Of
**>>> r@micom-solutions.de
**>>> Sent: Thursday, April 05, 2007 6:12 AM
**>>> To: r-help@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@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@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@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@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 272595Received on Sat Apr 07 17:14:50 2007______________________________________________ R-help@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.

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.1.8, at Sat 07 Apr 2007 - 11:30:59 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*