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

[Abby Spurdle]

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.