[R] Problem with MASS::fitdistr().

[J C Nash]

Peter is correct. I was about to reply when I saw his post.

It should be possible to suppress the Hessian call. I try to do this
generally in my optimx package as computing the Hessian by finite differences
uses a lot more compute-time than solving the optimization problem that
precedes the usual Hessian calculation.

It may be time to consider some update to optim() in that regard, or to
put in some exception handling. It isn't likely in any main-line part
of optim() but in the post-solution phase of the routines. I'd welcome
some discussion on that with a view to improvement, but am not sure if
it should be R-devel or R-package-dev lists or off-list.

John Nash




On 2020-04-26 4:53 a.m., peter dalgaard wrote:
> The optim() function has trouble calculating derivatives at/near the boundary, because it is using a simplistic finite-difference formula centered on the parameter. optimx::optimr() may work better.
> 
> -pd
> 
>> On 26 Apr 2020, at 09:02 , Abby Spurdle <spurdle.a using gmail.com> wrote:
>>
>> I ran your example.
>> It's possible that it's another bug in the optim function.
>>
>> Here's the optim call (from within fitdistr):
>>
>> stats::optim(x = c(1, 4, 1, 2, 3, 1, 1, 1, 2, 2, 2, 2, 1, 1,
>> 1, 1, 1, 4, 4, 3, 1, 2, 2, 1, 1, 3, 1, 1, 1, 4, 1, 1, 1, 1, 1, #more lines...
>> 1, 4, 1, 1, 1, 5, 5, 5, 4, 5, 2, 5, 5, 5, 5, 3, 3, 5, 4, 5, 2, #removed...
>> 4, 5, 5), topn = 5, lower = lwr, par = list(alpha = 1.010652,
>>    beta = 1.929018), fn = function(parm, ...) -sum(log(dens(parm, ...))),
>>    hessian = TRUE, method = "L-BFGS-B")
>>
>> And here's dens:
>>
>> function (parm, x, ...)
>> densfun(x, parm[1], parm[2], ...)
>>
>> I can't see any reason why it should call dens with parm[1] < lower[1].
>>
>> On Sun, Apr 26, 2020 at 5:50 PM Abby Spurdle <spurdle.a using gmail.com> wrote:
>>>
>>> I haven't run your example.
>>> I may try tomorrow-ish if no one else answers.
>>>
>>> But one question: Are you sure the "x" and "i" are correct in your function?
>>> It looks like a typo...
>>>
>>>
>>> On Sun, Apr 26, 2020 at 2:14 PM Rolf Turner <r.turner using auckland.ac.nz> wrote:
>>>>
>>>>
>>>> For some reason fitdistr() does not seem to be passing on the "..."
>>>> argument "lower" to optim() in the proper manner, and as result
>>>> falls over.
>>>>
>>>> Here is my example; note that data are attached in the file "x.txt".
>>>>
>>>> dhse <- function(i,alpha,beta,topn) {
>>>>    x <- seq(0,1,length=topn+2)[-c(1,topn+2)]
>>>>    p <- dbeta(x,alpha,beta)
>>>>    if(any(!is.finite(p))) browser()
>>>>    (p/sum(p))[i]
>>>> }
>>>>
>>>> lwr  <- rep(sqrt(.Machine$double.eps),2)
>>>> par0 <- c(alpha=1.010652,beta=1.929018)
>>>> x    <- dget("x.txt")
>>>> fit  <- MASS::fitdistr(x,densfun=dhse,topn=5,start=as.list(par0),
>>>>                       lower=lwr)
>>>>
>>>> The browser() in dhse() allows you to see that alpha has gone negative,
>>>> taking a value:
>>>>
>>>>>       alpha
>>>>> -0.001999985
>>>>
>>>> Continuing causes fitdistr() to fall over with the error message:
>>>>
>>>>> Error in stats::optim(x = c(1, 4, 1, 2, 3, 1, 1, 1, 2, 2, 2, 2, 1, 1,  :
>>>>>  non-finite finite-difference value [1]
>>>>
>>>> If I eschew using fitdistr() and "roll-my-own" as follows:
>>>>
>>>> foo <- function(par,x,topn){-sum(log(dhse(i=x,alpha=par[1],
>>>>                                           beta=par[2],
>>>>                                           topn=topn)))}
>>>>
>>>> fit <- optim(par0,fn=foo,method="L-BFGS-B",lower=lwr,topn=5,x=x)
>>>>
>>>> then optim() returns a result without complaint.
>>>>
>>>> Am I somehow messing up the syntax for fitdistr()?
>>>>
>>>> cheers,
>>>>
>>>> Rolf Turner
>>>>
>>>> P. S. I've tried supplying the "method" argument, method="L-BFGS-B"
>>>> explicitly to fitdistr(); doesn't seem to help.
>>>>
>>>> R.T.
>>>>
>>>> --
>>>> Honorary Research Fellow
>>>> Department of Statistics
>>>> University of Auckland
>>>> Phone: +64-9-373-7599 ext. 88276
>>>> ______________________________________________
>>>> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
>>>> 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 using r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> 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.
>