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

Question

我正在尝试应用期望最大化算法来估计缺失的计数数据，但 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 给出了答案，但数据不是连续的而是离散的。

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

Answer 1

正在定义 x0:

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

具有均值 lambda 的泊松分布的 Jeffreys/reference 先验是 1/sqrt(lambda)。根据观察值，这导致 lambda 具有具有形状参数 sum(x0) + 0.5 和速率参数 1/length(x0) 的伽马参考后验。您可以使用 lambda 的 n 个样本：

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))
}

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

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

Expectation Maximization using a Poisson likelihood function

r

expectation-maximization