11

我试图估计适合生态密度（即每面积生物量）数据的伽马分布参数。我一直在使用 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给我以下错误：

优化错误（x = c(6.46791148085828, 4060.54750836902, 99.6201565968665, : 非有限差分值 [1]

其他经历过此类问题并转向 stackoverflow 寻求帮助的人似乎已经找到了在他们的代码中添加“lower=”参数和/或删除零的解决方案。我发现如果我删除零观测值，R 将提供拟合参数，但我的印象是 gamma 分布覆盖了 0 <= x > inf 的范围（Forbes et al. 2011. Statistical Distributions）？

我对伽马分布的范围有错误的印象吗？或者是否还有其他关于 MLE 的问题（我是新手）。

4

1 回答 1

22

通过矩的方法得到一个粗略的估计（匹配均值=形状*比例和方差=形状*比例^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))
于 2013-07-12T14:16:31.480 回答