r - 带有 Gamma 分布的fitdist中的错误

Question

以下是我的代码:

library(fitdistrplus)

s <- c(11, 4, 2, 9, 3, 1, 2, 2, 3, 2, 2, 5, 8,3, 15, 3, 9, 22, 0, 4, 10, 1, 9, 10, 11,
 2, 8, 2, 6, 0, 15, 0 , 2, 11, 0, 6, 3, 5, 0, 7, 6, 0, 7, 1, 0, 6, 4, 1, 3, 5,
 2, 6, 0, 10, 6, 4, 1, 17, 0, 1, 0, 6, 6, 1, 5, 4, 8, 0, 1, 1, 5, 15, 14, 8, 1,
 3, 2, 9, 4, 4, 1, 2, 18, 0, 0, 10, 5, 0, 5, 0, 1, 2, 0, 5, 1, 1, 2, 3, 7)

o <- fitdist(s, "gamma", method = "mle")
summary(o)
plot(o)

错误显示:

Error in fitdist(s, "gamma", method = "mle") : the function mle failed to estimate the parameters, with the error code 100

Answer 1

Gamma分布不允许零值(对于0的响应，似然将评估为零，对数似然将是无限的)，除非shape参数正好是1.0(即指数分布-见下文)。 ..这是一个统计/数学问题，而不是编程问题。您将不得不对零值做出明智的选择。根据对您的应用程序有意义的内容，您可以(例如)

选择其他发行版以测试

排除零值(fitdist(s[s>0], ...))

将零值设置为合理的非零值(fitdist(replace(s,which(s==0),0.1),...)

其中最好的(如果有的话)取决于您的应用程序。

@Sandipan Dey的第一个答案(在数据集中保留零) 看起来很有意义，但实际上它停留在形状参数等于1的位置。

o <- fitdist(s, "exp", method = "mle")

给出与@Sandipan的代码相同的答案(除了它估计rate = 0.2161572，是针对Gamma分布估计的比例参数= 4.626262的倒数，这只是参数设置的变化)。如果您选择适合指数而不是Gamma，那很好-但是您应该有目的地这样做，而不是偶然地...

为了说明包含零的拟合可能无法按预期方式工作，我将构造自己的负对数似然函数并显示每种情况的似然面。

mfun <- function(sh,sc,dd=s) {
  -sum(dgamma(dd,shape=sh,scale=sc,log=TRUE))
}

library(emdbook)  ## for curve3d() helper function

含零的表面:

cc1 <- curve3d(mfun(x,y),
      ## set up "shape" limits" so we evaluate
      ## exactly shape=1.000 ...
        xlim=c(0.55,3.55),
        n=c(41,41),
        ylim=c(2,5),
        sys3d="none")


png("gammazero1.png")
with(cc1,image(x,y,z))
dev.off()

在这种情况下，仅在shape = 1(即指数分布)处将表面定义为；白色区域代表无限对数似然。不是说shape = 1是最合适的，而是那是唯一合适的...

零值排除的表面:

cc2 <- curve3d(mfun(x,y,dd=s[s>0]),
               ## set up "shape" limits" so we evaluate
               ## exactly shape=1.000 ...
               xlim=c(0.55,3.55),
               n=c(41,41),
               ylim=c(2,5),
               sys3d="none")


png("gammazero2.png")
with(cc2,image(x,y,z))
with(cc2,contour(x,y,z,add=TRUE))
abline(v=1.0,lwd=2,lty=2)
dev.off()