在R中模拟贝叶斯定理出现超界值，该如何处理？

问题背景 #

贝叶斯定理公式为：
$p(A|B) = \frac{p(A) \times p(B|A)}{p(B)}$

我们的目标是计算阴天时下雨的概率$p(rain|overcast)$，基于20天观测数据：5天下雨，其中4天是阴天；总共10天阴天。用R的rbeta函数对三个概率做1000次模拟时，发现有22次计算结果超出0-1的合理范围。现提出三种处理方案，需分析其有效性与缺陷，并找到正确处理方式：

• 将大于1的值截断为1
• 重新模拟直到分子≤分母
• 重新模拟直到分母>分子

模拟核心代码 #

library(tidyverse)
set.seed(4434)
sims = 1000
p_rain = rbeta(sims,5,15)
p_overcast_rain = rbeta(sims, 4,1)
p_overcast = rbeta(sims, 10,10)

data <- tibble(p_rain, p_overcast, p_overcast_rain) %>% 
  mutate(numerator = p_rain * p_overcast_rain,
         p_rain_overcast = numerator / p_overcast)

summary(data$p_rain_overcast)
table(data$p_rain_overcast > 1)

问题根源 #

超界的本质是独立模拟三个概率违反了概率的逻辑约束：根据全概率公式，$p(overcast)$必须满足：
$p(overcast) = p(overcast|rain) \times p(rain) + p(overcast|not rain) \times (1-p(rain))$

而代码中直接独立模拟$p(overcast) \sim Beta(10,10)$，完全忽略了它和$p(rain)$、$p(overcast|rain)$的依赖关系，导致可能出现$p(overcast) < p(rain) \times p(overcast|rain)$的情况，此时计算出的后验概率就会大于1。

三种方案的缺陷分析 #

1. 截断为1：人为压缩了分布的上尾，会低估后验分布的均值和方差，无法真实反映结果的不确定性，导致统计推断出现偏差。
2. 重新模拟直到分子≤分母：属于拒绝采样，会丢弃$p(overcast)$较小的样本，但这些样本本身是符合先验分布的合理值，丢弃后会引入选择偏差，改变原本的分布特性，得到的结果不再是真实的后验分布。
3. 重新模拟直到分母>分子：和第二种方案本质相同，同样存在选择偏差问题，无法得到无偏的后验分布。

正确处理方式 #

核心是让三个概率满足全概率公式的约束，以下是两种可行方法：

方法一：联合模拟符合约束的概率 #

根据观测数据，非雨天共15天，其中阴天6天，因此$p(overcast|not rain)$的先验分布为$Beta(6,9)$。我们先模拟$p(rain)$和两个条件概率，再通过全概率公式计算$p(overcast)$，确保其满足逻辑约束：

library(tidyverse)
set.seed(4434)
sims = 1000
p_rain = rbeta(sims,5,15)
p_overcast_rain = rbeta(sims,4,1)
p_overcast_norain = rbeta(sims,6,9) # 非雨天15天，6天阴天
p_overcast = p_rain * p_overcast_rain + (1 - p_rain) * p_overcast_norain

data <- tibble(p_rain, p_overcast, p_overcast_rain) %>% 
  mutate(numerator = p_rain * p_overcast_rain,
         p_rain_overcast = numerator / p_overcast)

summary(data$p_rain_overcast)
table(data$p_rain_overcast > 1) # 此时无超界值

方法二：直接模拟后验分布 #

根据贝叶斯共轭先验的性质，我们可以直接推导$p(rain|overcast)$的后验分布。结合观测数据（雨天阴天4次、雨天非阴天1次、非雨天阴天6次、非雨天非阴天9次），后验分布为$Beta(4+5, 6+15)$，直接模拟该分布即可得到无超界的结果：

library(tidyverse)
set.seed(4434)
sims = 1000
p_rain_overcast = rbeta(sims, 4+5, 6+15)

summary(p_rain_overcast)
table(p_rain_overcast > 1) # 无超界值

内容的提问来源于stack exchange，提问作者Derek DeBellis