Problem with MASS::fitdistr().
Dear Ms. Spurdle , Thanks for looking into this. I think that you are correct in that it is a problem with the hessian calculation. It seems that fitdistr() explicitly sets hessian=TRUE, with no possibility of opting out. It also seems that optim() ignores the "lower" argument when computing the hessian. I tried setting hessian=TRUE in my roll-your-own call to optim() and got the same crash, with alpha = -0.001999985, which of course results in disaster. Prof. Nash from time to time points out, on this and similar lists, that optim() --- which I believe he wrote --- is subject to problems. Perhaps this is one of them. It makes sense that there would be difficulties with computing a hessian when the parameters are subject to constraints. It's not at all clear to me how or if these difficulties can be circumvented. I figured that it was kind of nice that fitdistr() provides standard errors for the parameter estimates, but this isn't really vital for what I am trying to do --- and if trying to find these standard errors makes the estimation procedure fall over, then clearly standard errors have to be ditched. Thanks again for looking into my problem. cheers, Rolf
On 26/04/20 7:29 pm, Abby Spurdle wrote:
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 at 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 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.