使用R难以拟合伽玛分布

Question

我试图估计符合生态密度(即每面积生物量)数据的伽马分布的参数.我一直在使用R中MASS包中的fitdistr()命令(版本3.0.0:x86_64-w64-mingw32/x64(64位)).这是分布参数的最大似然估计命令.

数据向量非常大,但汇总统计数据如下:

闵.= 0; 第一曲.= 87.67; 中位数= 199.5; 平均值= 1255; 差异= 2.79E + 07; 第三曲.= 385.6; 最大.= 33880

我用来运行MLE过程的代码是:

gdist <- fitdistr(data, dgamma, 
                  start=list(shape=1, scale=1/(mean(data))),lower=c(1,0.1))

R给我以下错误:

optim中的错误(x = c(6.46791148085828,4060.54750836902,99.6201565968665,:非有限的有限差分值[1]

其他遇到此类问题并转向stackoverflow寻求帮助的人似乎已经找到了将"lower ="参数添加到其代码中和/或删除零的解决方案.如果我删除零观测值,我发现R将提供拟合参数,但我的印象是伽玛分布涵盖范围0 <= x> inf(Forbes et al.2011.Statistical Distributions)？

我是否对伽马分布的范围产生了错误的印象？或者是否有一些我在MLE上缺少的其他问题(我是新手).

Answer 1

通过矩量法得到粗略估计(匹配均值=形状*比例和方差=形状*比例^ 2)我们有

mean <- 1255
var <- 2.79e7
shape = mean^2/var   ## 0.056
scale = var/mean     ## 22231

现在从此分发中生成一些数据:

set.seed(101)
x = rgamma(1e4,shape,scale=scale)
summary(x)
##     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##     0.00      0.00      0.06   1258.00     98.26 110600.00 

MASS::fitdistr(x,"gamma")  ## error

我强烈怀疑问题是底层optim调用假定参数(形状和比例,或形状和速率)大致相同,但它们并非如此.您可以通过扩展数据来解决这个问题:

(m <- MASS::fitdistr(x/2e4,"gamma"))  ## works fine
##      shape           rate    
##  0.0570282411   0.9067274280 
## (0.0005855527) (0.0390939393)

fitdistr 给出一个速率参数而不是一个比例参数:回到你想要的形状参数,反转并重新缩放......

1/coef(m)["rate"]*2e4  ## 22057

顺便说一下,模拟数据的分位数与您的数据不匹配的事实(例如,中位数x= 0.06,数据中位数为199)表明您的数据可能不适合Gamma -例如,可能有一些异常值影响均值和方差多于分位数？

我上面的PS使用了内置的'gamma'估算器,fitdistr而不是使用dgamma:基于矩量法的起始值,以及将数据缩放2e4,这是有效的(尽管NaNs produced除非我们指定,否则会发出警告lower)

 m2 <- MASS::fitdistr(x/2e4,dgamma,
        start=list(shape=0.05,scale=1), lower=c(0.01,0.01))