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

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)(应用估算EM可以给出答案，但数据不是连续的而是离散的。

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

Answer 1

定义x0:

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

具有平均lambda的泊松分布的Jeffreys/参考先验是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)

或者，由于伽玛-泊松复合分布可以公式化为负二项式(在积分出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值都已估算。每一列可以在进一步的分析中单独使用，然后可以汇总结果。