使用条件概率表进行抽样答案

【问题标题】：Sampling using conditional probability table使用条件概率表进行抽样
【发布时间】：2018-09-12 02:08:45
【问题描述】：

我正在尝试模拟某些描述“世界真实状态”（例如“红色”、“绿色”或“蓝色”）的离散变量及其指标，但描述起来有些不完美。

r_names <- c("real_R", "real_G", "real_B")

假设我对“现实”变量的分布有一些先验信念，我将用它来采样它。

r_probs <- c(0.3, 0.5, 0.2)
set.seed(100)
reality <- sample(seq_along(r_names), 10000, prob=r_probs, replace = TRUE)

现在，假设我有条件概率表，它规定了给定每个“现实”的指标值

ri_matrix <- matrix(c(0.7, 0.3, 0, 
                      0.2, 0.6, 0.2, 
                      0.05,0.15,0.8), byrow=TRUE,nrow = 3)
dimnames(ri_matrix) <- list(paste("real", r_names, sep="_"),
                        paste("ind", r_names, sep="_"))

ri_matrix

>#            ind_R ind_G ind_B
># real_Red    0.70  0.30   0.0
># real_Green  0.20  0.60   0.2
># real_Blue   0.05  0.15   0.8

由于base::sample() 没有针对prob 参数进行矢量化，我必须：

sample_cond <- function(r, rim){
  unlist(lapply(r, function(x) 
    sample(seq_len(ncol(rim)), 1, prob = rim[x,], replace = TRUE)))
 }

现在我可以使用条件概率矩阵对我的“指标”变量进行采样

set.seed(200)
indicator <- sample_cond(reality, ri_matrix)

只是为了确保分发结果符合预期：

prop.table(table(reality, indicator), margin = 1)

 #>        indicator
 #> reality          1          2          3
 #>       1 0.70043610 0.29956390 0.00000000
 #>       2 0.19976124 0.59331476 0.20692400
 #>       3 0.04365278 0.14400401 0.81234320

是否有更好（即更惯用和/或更有效）的方法来对以另一个离散随机变量为条件的离散变量进行采样？

更新：

正如@Mr.Flick 所建议的，这至少快了 50 倍，因为它重用了概率向量，而不是条件概率矩阵的重复子集。

sample_cond_group <- function(r, rim){
il <- mapply(function(x,y){sample(seq(ncol(rim)), length(x), prob = y, replace = TRUE)}, 
       x=split(r, r),
       y=split(rim, seq(nrow(rim))))
unsplit(il, r)
}

【问题讨论】：

有趣的是，在原始函数中使用 for 循环的运行速度平均比 unlist 包装的 lapply 快 28-29%。所以像这样：循环中的mm[i] <- sample(1:3, 1, prob=rim[r[i],]) 更快。 Whodathunkit？

标签： r simulation probability sampling

【解决方案1】：

您可以通过使用拆分/组合类型策略在每组中抽取所有随机样本来提高效率。这可能看起来像这样

simFun <- function(N, r_probs, ri_matrix) {
  stopifnot(length(r_probs) == nrow(ri_matrix))
  ind <- sample.int(length(r_probs), N, prob = r_probs, replace=TRUE)
  grp <- split(data.frame(ind), ind)
  unsplit(Map(function(data, r) {
    draw <-sample.int(ncol(ri_matrix), nrow(data), replace=TRUE, prob=ri_matrix[r, ])
    data.frame(data, draw)
    }, grp, as.numeric(names(grp))), ind)
}

你可以打电话给

simFun(10000, r_probs, ri_matrix)

【讨论】：