【问题标题】:Zero Inflated Poisson Distribution: Fail to estimate the parameter with error code 100零膨胀泊松分布:无法估计参数,错误代码为 100
【发布时间】:2021-06-01 04:10:42
【问题描述】:

以下是我正在处理的一种数据集:

data <- c(0, 1, 0, 11, 2, 0, 3, 0, 0, 2, 1, 3, 1, 0, 1, 0, 0, 0, 2, 3, 
0, 0, 0, 8, 1, 1, 1, 0, 1, 1, 2, 7, 0, 0, 0, 5, 2, 3, 6, 1, 1, 
5, 2, 9, 0, 0, 1, 21, 16, 2, 9, 6, 25, 2, 1, 12, 16, 14, 15, 
15, 6, 1, 12, 12, 13, 5, 5, 6, 4, 7, 11, 8, 4, 5, 8, 3, 8, 4, 
7, 4, 7, 2, 5, 6, 4, 5, 1, 0, 8, 5, 6, 8, 9, 8, 9, 7, 7, 9, 8, 
9, 4, 4, 7, 13, 9, 13, 12, 10, 9, 8, 7, 11, 5, 5, 0, 1, 33, 4, 
22, 19, 22, 9, 5, 4, 17, 7, 7, 4, 5, 3, 0, 0, 9, 3, 0, 0, 36, 
40, 5, 4, 0, 11, 0, 7, 5, 25, 39, 26, 4, 20, 12, 4, 17, 3, 22, 
12, 14, 8, 9, 11, 7, 11, 10, 9, 16, 6, 24, 8, 5, 6, 14, 3, 9, 
4, 1, 20, 0, 1, 7, 9, 0, 12, 2, 29, 56, 16, 8, 28, 0, 19, 25, 
35, 87, 56, 66, 60, 58, 14, 10, 12, 13, 13, 34, 26, 18, 13, 22, 
13, 12, 15, 41, 11, 11, 11, 5, 6, 7, 8, 8, 17, 16, 12, 21, 38, 
34, 10, 77, 41, 7, 12, 1, 16, 20, 8, 5, 2, 20, 7, 16, 12, 6, 
10, 31, 12)

我使用了fitdistrplus 包来使用最大似然拟合分布。
同样,gamlss 包定义了零膨胀泊松分布的 pdf、CDF。

library(fitdistrplus)   
library(gamlss)  

这是数据集的均值和标准差。

mu=mean(data)
sigma=sd(data)

数据集的平均值为 10.75,标准差为 13.050。我尝试使用 fitdist 函数拟合零膨胀泊松分布,如下所示。

fit_zip = fitdist(data, 'ZIP', start = list(mu = mu, sigma = sigma))

一旦指定的函数被编译,编译器就会抛出以下错误

Error in fitdist(data, "ZIP", start = list(mu = mu, sigma = sigma)) : 
  the function mle failed to estimate the parameters, 
                with the error code 100

很明显mu和sigma的起始值初始化有问题。我不确定在什么基础上初始化起始值。
任何帮助表示赞赏。

【问题讨论】:

  • 1) 学习help("zip")。参数sigma 不是标准差。 2) 回想一下泊松分布 E(x) == var(x)。然后比较mean(data[data != 0])var(data[data != 0])。您的数据不太可能来自零膨胀泊松分布。 3) 我认为您不能在此发行版中使用fitdist。首先,fitdist 与您在假设 sigma 是标准差时犯了同样的错误。那是可以解决的。但是,当我使用模拟数据时,我会收到错误,因为 d/pZIP 函数的行为与错误输入的预期不同。

标签: r distribution


【解决方案1】:

您可以使用VGAM::*zipois() 代替gamlss.dist::*ZIP()。起始参数有点不同,见?dzipois,不一定是你输入的平均值。

library('fitdistrplus')
library('VGAM')  
fitdist(x, 'zipois', start=list(lambda=1, pstr0=.1))
# $start.arg
# $start.arg$lambda
# [1] 1
# 
# $start.arg$pstr0
# [1] 0.1
# 
# 
# $fix.arg
# NULL
# 
# Fitting of the distribution ' zipois ' by maximum likelihood 
# Parameters:
#   estimate Std. Error
# lambda 12.2009631 0.23823650
# pstr0   0.1187056 0.02069464
# Warning messages:
#   1: In fitdist(x, "zipois", start = list(lambda = 1, pstr0 = 0.1)) :
#   The dzipois function should return a zero-length vector when input has length zero
# 2: In fitdist(x, "zipois", start = list(lambda = 1, pstr0 = 0.1)) :
#   The pzipois function should return a zero-length vector when input has length zero

不过,我不确定这个警告。

您可以使用不同的向量x1(如下)来比较它,gamlss.dist 不会失败。


数据:

x <- c(0, 1, 0, 11, 2, 0, 3, 0, 0, 2, 1, 3, 1, 0, 1, 0, 0, 0, 2, 3, 
0, 0, 0, 8, 1, 1, 1, 0, 1, 1, 2, 7, 0, 0, 0, 5, 2, 3, 6, 1, 1, 
5, 2, 9, 0, 0, 1, 21, 16, 2, 9, 6, 25, 2, 1, 12, 16, 14, 15, 
15, 6, 1, 12, 12, 13, 5, 5, 6, 4, 7, 11, 8, 4, 5, 8, 3, 8, 4, 
7, 4, 7, 2, 5, 6, 4, 5, 1, 0, 8, 5, 6, 8, 9, 8, 9, 7, 7, 9, 8, 
9, 4, 4, 7, 13, 9, 13, 12, 10, 9, 8, 7, 11, 5, 5, 0, 1, 33, 4, 
22, 19, 22, 9, 5, 4, 17, 7, 7, 4, 5, 3, 0, 0, 9, 3, 0, 0, 36, 
40, 5, 4, 0, 11, 0, 7, 5, 25, 39, 26, 4, 20, 12, 4, 17, 3, 22, 
12, 14, 8, 9, 11, 7, 11, 10, 9, 16, 6, 24, 8, 5, 6, 14, 3, 9, 
4, 1, 20, 0, 1, 7, 9, 0, 12, 2, 29, 56, 16, 8, 28, 0, 19, 25, 
35, 87, 56, 66, 60, 58, 14, 10, 12, 13, 13, 34, 26, 18, 13, 22, 
13, 12, 15, 41, 11, 11, 11, 5, 6, 7, 8, 8, 17, 16, 12, 21, 38, 
34, 10, 77, 41, 7, 12, 1, 16, 20, 8, 5, 2, 20, 7, 16, 12, 6, 
10, 31, 12)

x1 <- c(0, 63, 1, 4, 1, 44, 2, 2, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 3, 
0, 0, 2, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
6, 1, 11, 1, 1, 0, 0, 0, 2)

【讨论】:

  • 您确定这实际上适合零膨胀分布吗?我得到与print(fitdist(data, 'pois', start = list(lambda = 10), discrete = TRUE, method = "mle")) 完全相同的 lambda 和不确定性。但当然,OP 的数据可能不是零膨胀的。
  • @Roland 你说得对,谢谢,pstr0= 也需要优化,目前它只适合泊松,我会编辑它。
  • 似乎可以工作,但与 OP 的数据非常不匹配。也许添加一些情节? hist(data, freq = FALSE, breaks = 0:90 - 0.5); curve(dzipois(x, lambda = 12.2, pstr0 = 0.12), add = TRUE, col = "blue", from = 0, to = 90, type = "p", n = 91)
  • @Roland 当然,没有询问是否适合:) 另请参阅other question 中的图。
猜你喜欢
  • 1970-01-01
  • 2011-11-01
  • 2023-04-04
  • 1970-01-01
  • 2018-02-10
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多