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
> > ______________________________________________
> > R-help at 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.
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. > > > > > > -- > > > Honorary Research Fellow > > > Department of Statistics > > > University of Auckland > > > Phone: +64-9-373-7599 ext. 88276 > > > ______________________________________________ > > > R-help at 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.
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.
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 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. >>> >>> -- >>> Honorary Research Fellow >>> Department of Statistics >>> University of Auckland >>> Phone: +64-9-373-7599 ext. 88276 >>> ______________________________________________ >>> R-help at 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 at 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 at cbs.dk Priv: PDalgd at gmail.com
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 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. >>>> >>>> -- >>>> Honorary Research Fellow >>>> Department of Statistics >>>> University of Auckland >>>> Phone: +64-9-373-7599 ext. 88276 >>>> ______________________________________________ >>>> R-help at 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 at 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. >