错误代码100 :在pois和nbinom的适配包中使用fitdist

Question

我正在尝试将一个发行版与我的数据(ssdist；ssi是值)相匹配，但是我一直被抛出下面的错误代码，我不知道为什么。值得注意的是，错误只抛给pois和nbinom，它们是我认为最适合数据的发行版。

错误：：fitdist(ss$ssiss$name == s，distr = "nbinom")：函数mle无法估计参数，错误代码为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)

这里有一个指向我的数据的链接(我不确定如何在这里发布)：https://docs.google.com/spreadsheets/d/1B-9sygnfDd8rUjyN3ZO6DwpvMO95FLMi02RlKpDC7Ig/edit?usp=sharing

Answer 1

时间太长了，不能发表评论。

这里有很多问题。可以说，你所犯的错误只是冰山一角。

fitdist(...)的默认方法是最大似然估计(mle)。fitdist(...)使用optim(...)函数来做到这一点(实际上是为了最小化负日志的可能性，这相当于相同的事情)。此错误来自optim(...)，意味着最小化失败。当分发的选择很差，或者分发参数的初始值选择得很差时，这是非常常见的。在你的情况下，两者都有一点。

您的数据具有非整数值。负二项分布对正整数(包括0)有支持作用，因此根据定义，非整数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')

​

​

这确实突出了极端的尾巴。在这种情况下，国际海事组织的衡量标准是相当无用的。