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
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, 1, 3, 4, 1, 1, 1,
3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1,
1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 2, 2, 4, 1, 1,
1, 2, 1, 1, 2, 1, 4, 1, 1, 1, 1, 1, 1, 3, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 4, 1, 1, 1, 1, 1, 1, 1, 4, 5, 2, 1, 1,
3, 1, 5, 1, 3, 1, 5, 4, 1, 4, 5, 4, 1, 5, 2, 1, 5, 1, 1, 4, 3,
3, 1, 2, 5, 4, 3, 5, 1, 2, 1, 1, 1, 1, 3, 3, 1, 1, 4, 1, 2, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 1, 1, 1, 1, 1, 4, 2, 1, 5, 1, 2,
1, 1, 3, 3, 4, 4, 1, 4, 4, 1, 1, 3, 4, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 4, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 1,
4, 1, 3, 3, 1, 2, 1, 1, 1, 1, 2, 2, 2, 1, 1, 1, 1, 4, 4, 3, 5,
1, 1, 3, 4, 5, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 4, 5, 5, 1, 1,
1, 1, 1, 1, 1, 1, 2, 2, 3, 5, 4, 5, 1, 1, 5, 1, 4, 4, 4, 2, 5,
1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 3, 5, 4, 3, 2, 1, 4, 1,
5, 1, 1, 2, 1, 1, 1, 5, 2, 4, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 4, 4, 3, 1, 5, 4, 1, 1, 2, 1, 1, 1, 1, 3, 4, 2, 1, 3, 5,
1, 1, 1, 1, 1, 1, 1, 3, 4, 3, 2, 5, 1, 1, 1, 5, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 5, 4, 4, 5, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
5, 1, 5, 5, 1, 5, 1, 1, 1, 1, 4, 1, 3, 3, 1, 1, 1, 1, 1, 1, 1,
1, 1, 3, 1, 3, 4, 2, 3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 1,
1, 5, 4, 1, 1, 1, 1, 1, 1, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 2, 1, 1, 3, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 3,
3, 5, 3, 5, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 2,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 2, 1, 1, 1, 1, 1, 1,
1, 1, 3, 1, 1, 1, 2, 2, 1, 1, 1, 1, 5, 4, 1, 1, 4, 2, 2, 4, 1,
1, 3, 1, 1, 1, 2, 1, 4, 1, 3, 4, 2, 1, 1, 4, 3, 1, 1, 2, 3, 4,
1, 2, 3, 1, 1, 1, 1, 1, 1, 4, 1, 1, 4, 3, 1, 4, 3, 1, 1, 1, 4,
4, 3, 4, 5, 2, 1, 1, 1, 1, 1, 2, 1, 1, 4, 3, 1, 2, 1, 3, 1, 3,
2, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
2, 1, 1, 1, 1, 1, 1, 1, 5, 1, 2, 1, 4, 1, 4, 1, 4, 1, 2, 4, 1,
4, 4, 4, 3, 1, 1, 3, 4, 1, 1, 1, 1, 1, 1, 1, 3, 1, 1, 2, 3, 1,
2, 5, 1, 1, 1, 1, 2, 1, 1, 2, 2, 5, 5, 5, 3, 4, 4, 1, 1, 2, 2,
3, 3, 1, 4, 3, 1, 2, 1, 1, 3, 1, 3, 1, 1, 4, 3, 5, 1, 4, 3, 1,
1, 4, 5, 1, 1, 5, 1, 1, 1, 5, 1, 1, 1, 5, 1, 2, 1, 1, 1, 1, 1,
1, 1, 2, 5, 5, 2, 5, 5, 5, 2, 1, 1, 1, 1, 3, 3, 1, 1, 2, 2, 4,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 1, 1, 1, 5, 1, 1, 1, 1, 1,
1, 1, 2, 2, 4, 1, 5, 4, 5, 2, 1, 1, 1, 1, 2, 3, 1, 1, 2, 5, 5,
5, 5, 5, 1, 1, 1, 1, 4, 1, 2, 1, 1, 1, 2, 1, 1, 1, 2, 1, 3, 1,
1, 4, 1, 4, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 4, 4, 2, 5, 4,
1, 4, 1, 1, 1, 1, 1, 1, 1, 5, 1, 1, 2, 2, 1, 1, 1, 1, 1, 1, 1,
1, 3, 1, 1, 1, 1, 2, 1, 1, 2, 1, 1, 1, 2, 1, 1, 1, 1, 1, 2, 1,
4, 4, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 4, 1, 1, 1, 3, 4, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 2, 3, 1, 1, 1, 1, 1, 3, 3, 1, 1, 3, 2, 2, 4, 2, 1,
1, 1, 3, 4, 4, 1, 3, 1, 5, 1, 1, 3, 1, 3, 1, 1, 1, 4, 1, 3, 1,
3, 1, 1, 1, 1, 1, 3, 4, 5, 1, 2, 1, 4, 2, 4, 1, 1, 1, 1, 3, 3,
1, 1, 3, 3, 1, 1, 1, 2, 1, 4, 4, 2, 2, 3, 1, 3, 1, 1, 1, 2, 4,
1, 3, 1, 1, 1, 1, 3, 1, 1, 1, 2, 2, 1, 4, 1, 1, 1, 5, 1, 1, 3,
5, 1, 2, 2, 1, 5, 1, 3, 1, 3, 4, 3, 5, 2, 3, 2, 4, 4, 3, 5, 4,
4, 4, 1, 1, 3, 4, 1, 1, 3, 4, 1, 1, 5, 1, 4, 4, 3, 5, 3, 5, 2,
5, 5, 5, 2, 3, 5, 5, 5, 5, 5, 5, 4, 4, 1, 2, 1, 2, 5, 1, 2, 1,
2, 1, 1, 1, 4, 1, 2, 1, 1, 2, 4, 4, 4, 1, 3, 1, 2, 4, 3, 5, 1,
2, 5, 1, 2, 1, 1, 5, 5, 2, 2, 1, 2, 5, 5, 5, 5, 5, 4, 3, 5, 1,
1, 2, 1, 3, 2, 4, 1, 3, 1, 1, 2, 2, 3, 3, 3, 1, 1, 4, 4, 4, 1,
1, 1, 5, 1, 1, 1, 3, 1, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
5, 5, 5, 5, 5, 4, 5, 5, 4, 4, 1, 1, 1, 1, 3, 2, 1, 1, 1, 1, 2,
3, 1, 1, 4, 1, 1, 3, 3, 1, 2, 1, 1, 3, 1, 4, 5, 5, 1, 5, 4, 1,
4, 1, 1, 1, 2, 1, 2, 4, 1, 1, 4, 1, 1, 1, 1, 1, 1, 2, 1, 2, 3,
1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 1,
1, 1, 1, 1, 4, 1, 1, 1, 5, 2, 3, 5, 1, 1, 2, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 3, 1, 2, 1, 1, 1, 1, 1, 1, 3, 1, 1, 1, 4, 1, 2,
1, 2, 1, 3, 4, 1, 1, 2, 3, 3, 1, 1, 3, 1, 5, 1, 3, 4, 1, 3, 1,
1, 1, 1, 1, 1, 2, 1, 5, 1, 1, 2, 1, 1, 1, 1, 4, 3, 3, 2, 5, 1,
4, 2, 1, 1, 1, 1, 3, 2, 2, 1, 1, 1, 1, 4, 1, 2, 5, 2, 3, 3, 1,
2, 3, 1, 5, 1, 1, 1, 1, 3, 1, 4, 1, 1, 1, 1, 1, 3, 1, 1, 1, 2,
4, 1, 4, 2, 1, 4, 1, 3, 1, 4, 2, 3, 4, 1, 3, 4, 2, 1, 3, 1, 3,
3, 1, 4, 3, 1, 3, 4, 1, 1, 5, 4, 4, 5, 4, 5, 5, 5, 5, 4, 5, 3,
5, 4, 2, 4, 4, 4, 5, 1, 1, 3, 2, 4, 1, 1, 1, 1, 1, 1, 4, 1, 1,
2, 2, 2, 1, 1, 1, 4, 1, 4, 1, 2, 1, 5, 5, 1, 5, 3, 5, 1, 5, 5,
5, 5, 5, 5, 5, 5, 5, 5, 5, 3, 5, 5, 5, 1, 2, 2, 2, 3, 1, 3, 2,
3, 1, 4, 1, 1, 3, 4, 2, 1, 1, 4, 3, 5, 5, 1, 1, 4, 1, 1, 1, 5,
5, 5, 4, 5, 2, 5, 5, 5, 5, 3, 3, 5, 4, 5, 2, 4, 5, 5)
______________________________________________
R-help@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.
