No zero values and I still get: stats - initial value in 'vmmin' is not finite

Question

I am trying to find parameters of Beta distribution. It is part of a larger loop (found it here on Stackoverflow) that fits several other distributions without difficulties, however, this particular one is giving me hard time. The problem is with fitdistr function.

I have no zero values. In the documentation of stats package it states that the use of the parameter ncp=0 could help solve the problem of values that are too close to zero. I tried but it doesn't work either.

What am I doing wrong?

Tks,

P.D, I estimated the parameters using other way. The thing is, this one is by far more convenient and flexible since it is part of a function.

...
else if(distrib[[i]] == "beta"){
      gf_shape = "beta"
      fd_b <- fitdistr(data1,densfun=dbeta , start=list(shape1=0.5,shape2=0.5))
      est_shape1 = fd_b$estimate[[1]]
      est_shape2 = fd_b$estimate[[2]]
      
      ks = ks.test(data1, "pbeta", shape1=est_shape1, shape2=est_shape2)
      # add to results
      results[i,] = c(gf_shape, est_shape1, est_shape2, ks$statistic, ks$p.value) 
    }

Here is the error code when I run the fitdistr function:

Error in stats::optim(x = c(0.8606, 0.6833, 0.7697, 0.8082, 0.9736, 0.3547,  : 
  initial value in 'vmmin' is not finite

Here is the data contained in the variable data1.

[1] 0.8606 0.6833 0.7697 0.8082 0.9736 0.3547 0.5396 0.6068 0.9673 0.1232 0.1230 0.9062
 [13] 0.8581 0.9195 0.8718 0.2995 0.9588 0.3186 0.3841 0.4958 0.2380 0.1067 0.6461 0.4009
 [25] 0.9560 0.3274 0.7687 0.2633 0.0471 0.0697 0.2624 0.5604 0.9674 0.2222 0.7014 0.4497
 [37] 0.9338 0.6897 0.7774 0.1213 0.5236 0.9935 0.9979 0.5728 0.8919 0.9738 0.4148 0.7511
 [49] 0.6380 0.6931 0.3289 0.8363 0.6721 0.9462 0.5552 0.6020 0.9277 0.5347 0.7296 0.2748
 [61] 0.2404 0.3865 0.8381 0.4434 0.7989 0.4724 0.1292 0.2984 0.6156 0.6434 0.9051 0.9274
 [73] 0.3271 0.5559 0.4633 0.5042 0.2636 0.3432 0.5293 0.5310 0.8245 0.8818 0.7003 0.5456
 [85] 0.7186 0.6345 0.4206 0.1410 0.9882 0.9615 0.5349 0.0406 0.4825 0.8020 0.2100 0.0412
 [97] 0.7347 0.6702 0.8903 0.6576 0.6181 0.5562 0.2012 0.7089 0.1347 0.8188 0.7731 0.5648
[109] 1.0000 0.9612 0.5467 0.3996 0.8933 0.6008 0.9471 0.7940 0.6513 0.9878 0.6316 0.9577
[121] 0.8700 0.9101 1.0000 0.8311 0.8943 0.7236 0.7754 0.1535 0.9665 0.8847 0.6436 0.8146
[133] 0.9004 0.8190 0.5694 0.2346 0.0975 0.9178 0.4591 0.5278 0.7607 0.7919 0.8128 0.5495
[145] 0.8940 0.9158 0.4819 0.9168 0.7689 0.5053 0.9889 0.9656 0.4006 0.7682 0.7685 0.9740
[157] 0.6632 0.1044 0.8774 0.9524 0.6153 0.1046 0.1184 0.8995 0.6515 0.1793 0.7582 0.6055
[169] 0.8830 0.5494 0.6533 0.5356 0.9444 0.8056 0.5432 0.7286 0.9980 0.5776 0.8602 0.7249
[181] 0.5489 0.9401 0.8686 0.5991 0.8852 0.3737 0.6637 0.9300 0.4563 0.2764 0.2781 0.9860
[193] 0.7199 0.7748 0.3296 0.9003 0.6646 0.8454 0.8783 0.1340 0.4407 0.7323 0.9592 0.7804
[205] 0.5309 0.5338 0.2196 0.4862

> min(data1)
[1] 0.0406
> max(data1)
[1] 1

Answer 1

You could first try with maximum likelihood, if that fails, try matching moments,

library(fitdistrplus)

fitFunc <- function(dat) {
    require(fitdistrplus)
    res <- tryCatch({
        fitdist(dat, "beta", start=list(shape1=0.5, shape2=0.5), method="mle")
    }, error=function(e) return ( NULL ) )
    if (is.null(res))
        res <- fitdist(dat, "beta", start=list(shape1=0.5, shape2=0.5), method="mme")
    return(res)
}

res <- fitFunc(data1)
res
# .. some error output ...
# Parameters:
#         estimate
# shape1 1.4276936
# shape2 0.8329679

hist(data1, freq=F)
tst <- rbeta(100, coef(res)[1], coef(res)[2])
points(density(tst), type="l")