Re: [R] maximum likelihood estimation of 5 parameters

From: Charles C. Berry <cberry_at_tajo.ucsd.edu>
Date: Sat 06 Jan 2007 - 01:06:54 GMT

Franco,

It is up to the user of mle() to define a function that is numerically well behaved - even near the boundaries - and/or to provide suitable boundary constraints.

Using gamma() in situations in which the argument could be 500 (say) will cause problems. It appears that the maximum of fr is attained for values for which gamma(a2+n) cannot be evaluated using machine arithmetic (and for very small values of 'p', too). Does this function truly have a finite maximum for these data?

Dividing by very small numbers is also a potential killer.

Trying to exponentiate potentially negative numbers is also asking for trouble.

So, you need to write a more robust version of 'fr'.

Some suggestions:

 	Use lgamma() or choose() or lchoose() rather than gamma() or
 	factorial()

 	Try to add and subtract logarithms in preference to
 	multiplying and dividing as cancellations usually are more
 	accurate and you are less likely to run into machine accuracy
 	issues.

 	Replace p with plogis(p.x) and consider bounding p.x to keep
 	plogis(x) away from 0 and 1

 	Establish boundary constraints to keep
 	things like (1+(e/b1)^a1)*(1+(b1/e)^n)
 	from causing problems. Consider log1p(), etc or writing this
 	out longhand and doing some cancellations.

 	Worry about fnscale. It seems unlikely that all args are indeed on
 	the same scale.

 	Use stingy boundary constrants to get a better set of starting
 	values. Progressively relax the constraints and watch how the
 	parameter changes.

If after all that, your new fr function still hands you non-finite values, then you need to investigate to find out where they are coming from.

There are many ways to do this, but broadly speaking you are 'debugging' your code, so use the manuals and RSiteSearch to figure out what you need to do to become skilled in debugging.

On Fri, 5 Jan 2007, francogrex wrote:

>
> No I have not forgotten to use a negative fnscale to optimize, so as you
> suggest I will post some parts of the code I am running to show you the
> errors:
>
>> n
> [1] 3 1 4 54 6 58 20 14 3 14 4 65 1 7 9 10 2 4
> 66
> [20] 5 9 7 12 7 55 105 2 5 10 55 5 28 1 1 6 2 1
> 30
> [39] 6 49 7 21 8 7
>> e
> [1] 21.763201 1.209070 4.836270 32.644798 19.546600 24.584400 30.226700
> [8] 6.045340 14.010100 3.113350 21.015100 12.583100 15.826200 19.458401
> [15] 3.891690 1.329970 0.241814 3.143580 13.057900 0.725441 18.136000
> [22] 2.187660 6.319900 1.701510 29.654900 36.460999 7.292190 1.215370
> [29] 3.209070 19.995001 11.972300 3.455920 0.138539 0.113350 1.360200
> [36] 1.889170 1.518890 18.226700 4.050380 27.340099 1.181360 16.370300
> [43] 20.589399 25.314899
>
>> fr<-function(a1,b1,a2,b2,p){
> +
> + w<-((gamma(a1+n)))/((gamma(a1)*factorial(n))*(1+(e/b1)^a1)*(1+(b1/e)^n))
> + z<-((gamma(a2+n)))/((gamma(a2)*factorial(n))*(1+(e/b2)^a2)*(1+(b2/e)^n))
> +
> + sum (log( (p*w)+ ((1-p)*z) ))
> +
> + }
>>
>> mle((fr),
>> start=list(a1=0.2,b1=0.1,a2=2,b2=4,p=0.33),method="BFGS",control=list(fnscale=-1))
> Error in optim(start, f, method = method, hessian = TRUE, ...) :
> non-finite finite-difference value [2]
>
> And with the L-BFGS-B:
>
> Error in optim(start, f, method = method, hessian = TRUE, ...) :
> L-BFGS-B needs finite values of 'fn'
>
> AND WITH NELDER-MEAD it doesn't work either (same error), but when I change
> intial parameters (though I shouldn't, it gives something very weird
> (negatives or sometimes huge values).
>
> Call:
> mle(minuslogl = (fr), start = list(a1 = 1, b1 = 1, a2 = 10, b2 = 10,
> p = 0.9), method = "Nelder-Mead", control = list(fnscale = -1))
>
> Coefficients:
> a1 b1 a2 b2 p
> -2.5035823 0.6236359 26.5562988 12.9604112 -0.1383767
>
> Thanks
>
>
>
> Ravi Varadhan wrote:
>>
>> Franco,
>> Is it possible that you have failed to provide the negative of
>> loglikelihood
>> to "optim", since optim, by default, minimizes a function? If you want to
>> do this withput redefining the log-likelihood, you should set fnscale= -1
>> (as hinted by Prof. Ripley). This would turn the problem into a
>> maximization problem.
>>
>> If this doesn't work, you should provide more details (a reproducible code
>> with actual error message).
>>
>> Ravi.
>>
>> ----------------------------------------------------------------------------
>> -------
>>
>> Ravi Varadhan, Ph.D.
>>
>> Assistant Professor, The Center on Aging and Health
>>
>> Division of Geriatric Medicine and Gerontology
>>
>> Johns Hopkins University
>>
>> Ph: (410) 502-2619
>>
>> Fax: (410) 614-9625
>>
>> Email: rvaradhan@jhmi.edu
>>
>> Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
>>
>>
>>
>> ----------------------------------------------------------------------------
>> --------
>>
>> -----Original Message-----
>> From: r-help-bounces@stat.math.ethz.ch
>> [mailto:r-help-bounces@stat.math.ethz.ch] On Behalf Of francogrex
>> Sent: Friday, January 05, 2007 10:42 AM
>> To: r-help@stat.math.ethz.ch
>> Subject: Re: [R] maximum likelihood estimation of 5 parameters
>>
>>
>>
>> Franco,
>> You can provide lower and upper bounds on the parameters if you use optim
>> with method="L-BFGS-B".
>> Hth, Ingmar
>>
>> Thanks, but when I use L-BFGS-B it tells me that there is an error in
>> optim(start, f, method = method, hessian = TRUE, ...) : L-BFGS-B needs
>> finite values of 'fn'
>>
>> --
>> View this message in context:
>> http://www.nabble.com/maximum-likelihood-estimation-of-5-parameters-tf292536
>> 4.html#a8180120
>> Sent from the R help mailing list archive at Nabble.com.
>>
>> ______________________________________________
>> 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.
>>
>>
>
> --
> View this message in context: http://www.nabble.com/maximum-likelihood-estimation-of-5-parameters-tf2925364.html#a8186869
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> 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.
>

Charles C. Berry                        (858) 534-2098
                                          Dept of Family/Preventive Medicine
E mailto:cberry@tajo.ucsd.edu	         UC San Diego
http://biostat.ucsd.edu/~cberry/         La Jolla, San Diego 92093-0717

______________________________________________
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. Received on Sat Jan 06 12:25:38 2007

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 06 Jan 2007 - 01:30:26 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.