在 pois 和 nbinom 的 fitdistrplus 包中使用 fitdist 时出现错误代码 100

Error code 100 using fitdist in the fitdistrplus package for pois and nbinom

我正在尝试对我的数据进行分布(ssdistssi 是值),但是我一直收到下面的错误代码,我不知道为什么。值得注意的是,错误仅针对 poisnbinom 抛出,我认为这些分布最适合数据。

Error in fitdistrplus::fitdist(ss$ssi[ss$name == s], distr = "nbinom") : the function mle failed to estimate the parameters, with the error code 100

这是我的代码:

ssdist <- read_excel(here("Data/Raw/ssdist.xlsx")) %>%
  verify(has_all_names("name", "ssi")) %>% 
  assert(name) %>% 
  assert(is.numeric, ssi)

library(fitdistrplus)

par0 <- par(mfrow = c(3,2))
for(s in unique(ssdist$name)) {
  fitdistrplus::descdist(ss$ssi[ss$name == s], discrete = TRUE)
  title(sub = s)
}
par(par0)

这是给我问题的代码块

    par0 <- par(mfrow = c(2,2))
for(s in unique(ssdist$name)) {
  f <- fitdistrplus::fitdist(ss$ssi[ss$name == s], distr = "nbinom")
  fitdistrplus::cdfcomp(f, main = "Cumulative distribution: Data vs. theoretical")
  title(sub = s)
  mtext(side = 3, line = 0,
        text = paste0("p-value: ", round(fitdistrplus::gofstat(f)$chisqpvalue, 4)))
}
par(par0)

这是我的数据的 link(我不确定如何在此处 post):https://docs.google.com/spreadsheets/d/1B-9sygnfDd8rUjyN3ZO6DwpvMO95FLMi02RlKpDC7Ig/edit?usp=sharing

评论太长了。

这里有很多问题。您得到的错误只是冰山一角,可以这么说。

fitdist(...) 中的默认方法是最大似然估计 (mle)。 fitdist(...) 使用 optim(...) 函数来执行此操作(实际上是为了最小化负对数似然,这相当于同一件事)。此错误来自 optim(...),表示最小化失败。当分布选择不当或分布参数的初始值选择不当时,这种情况很常见。在你的情况下,两者兼而有之。

您的数据有 non-integer 个值。负二项分布支持正整数(包括 0),因此根据定义 non-integer x 的概率密度为 0。如果您阅读有关 dnbinom(...) 的文档,您会看到这一点。

第二个问题是你的数据有很长的尾巴:

library(data.table)
library(ggplot2)
setDT(ss)
ggplot(ss, aes(x=ssi)) + 
  geom_histogram(aes(fill=name), color='grey80') + 
  scale_y_continuous(breaks=2*(0:5))+
  scale_fill_discrete(guide='none')+
  facet_wrap(~name)

因此,指数族中的任何分布能否提供足够的拟合值得怀疑。一种可能性是对数正态分布,但此分布支持正实数,并且您的数据为零。

另一种可能是威布尔分布:

get.params <- function(x) {
  start  <- list(shape=1, scale=1)
  f      <- fitdist(x, distr = dweibull, start=start, method='mge')
  return(c(as.list(f$estimate)))
}
params <- ss[, get.params(ssi), by=.(name)]

这似乎在处理长尾巴方面做得不错:

ss[params, c('shape', 'scale'):=.(i.shape, i.scale), on=.(name)]
setorder(ss, name, ssi)
ss[, sample:=quantile(ssi, probs = ppoints(.N)), by=.(name)]
ss[, theoretical:=qweibull(ppoints(.N), shape, scale), by=.(name)]
ss[, dens:=dweibull(ssi, shape, scale), by=.(name)]
ss[, smpl.CDF:=ecdf(ssi)(ssi), by=.(name)]
ss[, theor.CDF:=pweibull(ssi, shape, scale), by=.(name)]
ggplot(ss, aes(x=ssi))+
  geom_line(aes(y=theor.CDF, color=name))+
  geom_point(aes(y=smpl.CDF), shape=21)+
  scale_color_discrete(guide='none')+
  labs(title='CDF Plots')+
  facet_wrap(~name)

但即使在这里你真的应该看看 q-q 图:

ggplot(ss, aes(x=theoretical, y=sample))+
  geom_point(aes(color=name))+
  geom_abline(color='blue', linetype='dashed')+
  scale_color_discrete(guide='none')+
  labs(title='QQ Plots')+
  facet_wrap(~name, scales = 'free_x')

真正突出了极端的尾巴。 IMO chi-sq gof 指标在这种情况下毫无用处。