使用R进行伽马分布拟合的困难

Question

我正在尝试估计适用于生态密度（即单位面积的生物量）数据的伽玛分布参数。我一直在使用 R 中 MASS 包中的 fitdistr() 命令进行最大似然估计。数据向量很大，但摘要统计信息如下:

Min. = 0; 1st Qu. = 87.67; Median = 199.5; Mean = 1255; Variance = 2.79E+07; 3rd Qu. = 385.6; Max. = 33880

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

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

R 给了我以下错误：

Error in optim(x = c(6.46791148085828, 4060.54750836902, 99.6201565968665, : non-finite finite-difference value [1]

其他遇到此类问题并在 stackoverflow 上寻求帮助的人似乎已经通过在代码中添加“lower=”参数和/或删除零来找到了解决方案。我发现如果删除零观察值，R将为拟合提供参数，但我一直认为伽马分布覆盖了0 <= x > inf 的范围（Forbes 等人，2011年。统计分布）？

我对伽马分布范围的印象是错误的吗？还是我在MLE方面存在其他问题（我是初学者）？

- niafall

Answer 1

通过矩估计法 (匹配均值=形状×比例尺和方差=形状×比例尺的平方) 得到一个粗略的估计，我们有：

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），这表明您的数据可能不太适合伽马分布——例如，可能有一些异常值影响了均值和方差，而不是分位数？

PS：上面我使用了fitdistr中的内置“gamma”估计器，而不是使用dgamma。使用矩法确定起始值，并将数据缩放2e4倍可以实现这一点（尽管它会给出一个关于产生NaNs的警告，除非我们指定lower参数）。

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

- Ben Bolker

---

如果形状参数小于1，我会怀疑数据是否适合伽马分布。虽然伽马分布可以允许这种情况，但是根据我的经验，特别是当与大规模数据一起使用时，这种值意味着数据可能尾部太重了，因此广义帕累托或极值分布可能更适合。 - Hong Ooi