R simulation

我们前面主要介绍了simulation在样本量计算中的应用：非劣效、生存分析、有效事件数，实际上，我们在学习概率论时候，解的那些数学题，同样也可以用simulation的方式来逼近一个正确的答案，从而避开了复杂的数学计算与推导。

一旦涉及到概率，就会涉及到比如组合数学，比如：

从40个球中抽取5个球，至少有一个红球、一个蓝球、一个绿球的概率是多少？

这个问题用数学公式计算起来非常复杂，但是用simulation的方式就非常简单了。

数学概率计算

上面这个问题也是比较典型的题目，我们首先用数学的方式来计算一下：用总体概率减去不符合条件的概率。

$$ \begin{aligned} \text{Pr}(\text{at least one red, blue, and green}) &= 1 - \frac{\text{Ways to get no red or no blue or no green}}{\text{Ways to get 5 balls from 40}} \[10pt] &= 1 - \dfrac{ \begin{array}{@{}c@{}} (\text{Ways to get no red}) + (\text{Ways to get no blue}) + (\text{Ways to get no green}) - \ (\text{Ways to get no red or blue}) - (\text{Ways to get no red or green}) - (\text{Ways to get no blue or green}) \end{array} }{\text{Ways to get 5 balls from 40}} \[10pt] &= 1 - \dfrac{ \dbinom{30}{5} + \dbinom{30}{5} + \dbinom{20}{5} - \dbinom{20}{5} - \dbinom{10}{5} - \dbinom{10}{5} }{ \dbinom{40}{5} } \end{aligned} $$

用R来解：

#| label: full-probability-math-r
no_red <- choose(30, 5)
no_blue <- choose(30, 5)
no_green <- choose(20, 5)

no_red_blue <- choose(20, 5)
no_red_green <- choose(10, 5)
no_blue_green <- choose(10, 5)

total_ways <- choose(40, 5)

prob_real <- 1 - (no_red + no_blue + no_green - no_red_blue - no_red_green - no_blue_green) / total_ways
prob_real

这样我们能算出概率，但是这样太麻烦了，而且很容易出错。

暴力模拟

简单粗暴地解这个题，就是用计算机模拟抽球的过程，然后重复很多次，最后计算符合条件的概率。这其实也是现代贝叶斯统计的核心思想，计算复杂的积分来找到后验分布有时候太难了，所以我们可以在后验分布的可能空间中用马尔科夫链蒙特卡洛（MCMC）的方法，直到其收敛。

#| label: simulated-urn
urn= c(rep('red', 10), rep('blue', 10), rep('green', 20))
simulate= function(n, num_simulations = 2000000) {
  results <- map_lgl(1:num_simulations, ~{
    draw <- sample(urn, n)
    all(c('red', 'blue', 'green') %in% draw)
  })
  mean(results)
}
prob_simulated= withr::with_seed(12345, {
  simulate(5)
})
prob_simulated

可以拿到一样的概率。

这样我们就得到了一个答案，而且不需要做任何数学计算。通过重新构造问题的场景数据，也提供了一种解决问题的思路。

代码已经放进了星球里。

Did you find this page helpful? Consider sharing it 🙌