Problem with MASS::fitdistr().

str (fit)

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 at 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 at 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