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
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: x.txt
URL:
<https://stat.ethz.ch/pipermail/r-help/attachments/20200426/f142865e/attachment.txt>
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.
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.