在 pois 和 nbinom 的 fitdistrplus 包中使用 fitdist 时出现错误代码 100
Error code 100 using fitdist in the fitdistrplus package for pois and nbinom
我正在尝试对我的数据进行分布(ssdist
;ssi
是值),但是我一直收到下面的错误代码,我不知道为什么。值得注意的是,错误仅针对 pois
和 nbinom
抛出,我认为这些分布最适合数据。
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 指标在这种情况下毫无用处。
我正在尝试对我的数据进行分布(ssdist
;ssi
是值),但是我一直收到下面的错误代码,我不知道为什么。值得注意的是,错误仅针对 pois
和 nbinom
抛出,我认为这些分布最适合数据。
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 指标在这种情况下毫无用处。