2015-06-17 23 views
1

我試圖找到Beta分佈的參數。它是一個更大的循環的一部分(在Stackoverflow中可以找到它),它可以毫無困難地適合其他幾個發行版,但是,這個特殊的循環給了我很大的困難。問題是與fitdistr函數。沒有零值,我仍然得到:stats - 'vmmin'中的初始值不是有限的

我沒有零值。在統計軟件包的文檔中,它指出使用參數ncp = 0可以幫助解決太接近於零的值的問題。我嘗試過,但它也不起作用。

我在做什麼錯?

Tks,

P.D,我用其他方法估計參數。事情是這樣的,因爲它是一個函數的一部分,所以它更加方便和靈活。

... 
 
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) 
 
    }

以下是錯誤代碼時,我運行fitdistr功能:

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

這裏是包含在可變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

+0

數據中的零和數值都可能導致此錯誤,分配的選擇可能不正確。嘗試省略你的1,看看是否修復它。 – jeremycg

回答

0

你可以先用最大似然嘗試,如果失敗,嘗試匹配的時刻,

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") 

enter image description here

+0

用statd包中的原始函數代替fitdistrplus工作。謝謝。 –