使用泊松似然函数的期望最大化答案

【问题标题】：Expectation Maximization using a Poisson likelihood function使用泊松似然函数的期望最大化
【发布时间】：2022-02-11 23:57:08
【问题描述】：

我正在尝试应用期望最大化算法来估计缺失的计数数据，但 R 中的所有包（例如 missMethods）都假设为多元高斯分布。假设泊松分布，我将如何应用期望最大化算法来估计缺失的计数数据？

假设我们有如下所示的数据：

x <- c(100,  96,  79, 109, 111,  NA,  93,  95, 119,  90, 121,  96,  NA,  
       NA,  85,  95, 110,  97,  87, 104, 101,  87,  87,  NA,  89,  NA, 
       113,  NA,  95,  NA, 119, 115,  NA, 105,  NA,  80,  90, 108,  90,  
       99, 111,  93,  99,  NA,  87,  89,  87, 126, 101, 106)

使用missMethods (missMethods::impute_EM(x, stochastic = FALSE)) 应用 impute_EM 会给出答案，但数据不是连续的而是离散的。

我知道像这样的问题需要一个最少的、可重复的例子，但老实说我不知道从哪里开始。甚至建议阅读以将我指向正确的方向也会有所帮助。

【问题讨论】：

您不希望计数数据的插补结果是“离散的”吗？毕竟，在许多情况下，你会将该输出提交给函数，如果你有非整数输入，那么这些函数至少会给你警告，最坏的情况是只会出错。
我很抱歉，但我对你的评论有点困惑。基于平均值或其他一些统计数据的插补与期望最大化不同。在将缺失值提交给另一个分析之前，我想使用更严格的方法来估计缺失值。

标签： r expectation-maximization

【解决方案1】：

定义x0：

x0 <- x[!is.na(x)]

平均为lambda 的泊松分布的杰弗里斯/参考先验是1/sqrt(lambda)。根据观察到的值，这导致lambda 具有带有形状参数sum(x0) + 0.5 和速率参数1/length(x0) 的伽马参考后验。您可以通过以下方式获取n 的lambda 样本：

lambda <- rgamma(n, sum(x0) + 0.5, length(x0))

然后对 n 缺失值 (xm) 进行采样

xm <- rpois(n, lambda)

或者，由于 Gamma-Poisson 复合分布可以表示为负二项式（在积分出 lambda 之后）：

xm <- rnbinom(n, sum(x0) + 0.5, length(x0)/(length(x0) + 1L))

作为一个函数：

MI_poisson <- function(x, n) {
  x0 <- x[!is.na(x)]
  rbind(matrix(x0, ncol = n, nrow = length(x0)),
        matrix(rnbinom(n*(length(x) - length(x0)), sum(x0) + 0.5, length(x0)/(length(x0) + 1L)), ncol = n))
}

这将返回一个包含n 列的矩阵，其中每列包含原始向量x，其中所有NA 值均已估算。每列可以单独用于进一步分析，然后可以汇总结果。

【讨论】：

您确实意识到这不是 EM 插补，但问题清楚地谈到了 EM 插补。
@Onyambu 不会产生非 NA 值的均值，因为它是单个样本空间，而 ML 是均值？
问题是这种方法没有捕捉到缺失数据中存在的固有随机性。这就是我想使用 EM 的原因。让 EM 算法适应非正态的指数分布有那么难吗？
您是否希望获得单个值来替换所有 NA 值（单一插补），还是希望每个 NA 独立替换为随机值（多重插补）？
我的理解是 EM 用于单次插补。请参阅修改后的答案以了解多重插补解决方案。