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

[Abby Spurdle]

fitdistr computes a Hessian matrix.
I think optim ignores the lower value computing the Hessian.

Here's the result after removing the Hessian and Hessian-dependent info:

> str (fit)
List of 3
 $ estimate: Named num [1:2] 0.0000000149 1.0797684972
  ..- attr(*, "names")= chr [1:2] "alpha" "beta"
 $ loglik  : num -2039
 $ n       : int 1529
 - attr(*, "class")= chr "fitdistr"


On Sun, Apr 26, 2020 at 7:02 PM 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.