r - 将泊松二项式分布拟合到数据

Question

当我尝试使用 fitdistrplus 来适应数据分布时，我总是遇到同样的错误。 MWE 如下。简而言之，我想对某些数据拟合泊松二项式分布。我正在使用 poisbin R 包来实现泊松二项式 p、d、q、r 函数(我也尝试过 poisbinom，但出现了同样的错误)。在 MWE 中，我创建了 dd，成功的向量。然后，我尝试使用 fitdist 来拟合给定 start 列表中起始值的分布。该错误表明(我认为)我给它的起始值的名称不在 dpoibin 函数中，这就是我陷入困境的地方。

library(fitdistrplus)
library(poibin)
set.seed(123)
dd <- rpoibin(10, pp=seq(0.1, 0.5, length.out=10))
ppp <- runif(10)
ret <- try(fitdistrplus::fitdist(dd, distr=dpoibin,
    start=list(pp = ppp)))

错误消息:

Error in checkparamlist(arg_startfix$start.arg, arg_startfix$fix.arg, : 'start' must specify names which are arguments to 'distr'.

Answer 1

该错误来自函数 fitdistrplus:::checkparamlist，该函数由 fitdist 调用，以确保列表中的名称传递给 start 与传递给 distr 的函数中的参数名称匹配。当您将 ppp 这样的向量作为 start 中的参数传递时，checkparamlist 通过附加一个整数来重命名向量的每个元素。这意味着参数名称将变为“pp1”、“pp2”、“pp3”，依此类推，直到“pp10”。由于没有传递名为 pp 的参数，因此会引发错误。

由于这个问题，我不确定是否有办法在 fitdist 中估计矢量化参数，但幸运的是，在这种情况下，我们可以轻松地自己拟合分布。

因为我们知道分布的平均值是

方差为

(Reference)

然后我们知道，如果我们有一个示例 dd，如果 pp 完全符合分布，则以下函数将返回 0:

objective <- function(pp) {
  abs(mean(dd) - sum(pp)) + abs(sum(pp * (1 - pp)) - var(dd))
}

为了演示这一点，让我们从 rpoibin 中获取一个更大的样本

set.seed(123)

dd  <- poibin::rpoibin(100000, pp=seq(0.1, 0.5, length.out=10))
ppp <- runif(10)

现在我们找到优化目标函数的一组值:

pp_opt <- optim(par = ppp, objective)$par

pp_opt
#>  [1] 0.45594175 0.08754997 0.54250499 0.28056432 0.30363457 0.28354584
#>  [7] 0.17861750 0.21109410 0.41562763 0.23920435

我们可以通过绘制直方图并将 dpoibin 的输出与我们为 pp 参数计算的值叠加来确认这是一个很好的拟合:

hist(dd, freq = FALSE, breaks = 0:11 - 0.5)
points(0:10, poibin::dpoibin(0:10, pp = pp_opt), col = "red")

请注意，对于 pp 的最佳值可能有很多解决方案，我们不应该期望得到 seq(0.1, 0.5, length.out = 10) 。首先，顺序并不重要。我们可以看到我们的 pp_opt 与 seq(0.1, 0.5, length.out = 10) 具有非常相似的均值和方差，这对于拟合分布

mean(seq(0.1, 0.5, length.out = 10))
#> [1] 0.3
mean(pp_opt)
#> [1] 0.2998285

sum((1 - pp_opt) * pp_opt)
#> [1] 1.930687
sum((1 - seq(0.1, 0.5, length.out = 10)) * seq(0.1, 0.5, length.out = 10))
#> [1] 1.937037

一般来说，由于排序以及无数组具有相同分布和计算方差的事实，不可能从给定样本中准确恢复 pp。

^{创建于 2023 年 7 月 18 日 reprex v2.0.2}