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

[peter dalgaard]

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.

-- 
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Office: A 4.23
Email: pd.mes using cbs.dk  Priv: PDalgd using gmail.com