首先我们先通过贝叶斯定理计算出理论值，这样你可以对比模拟结果是否正确：

理论计算结果 #

1. P(D|T)（阳性时患病的概率） #

根据贝叶斯公式：

P(D|T) = [P(T|D) * P(D)] / [P(T|D)*P(D) + P(T|Dc)*P(Dc)]

代入数值：

• P(D) = 0.0005（男性患病率）
• P(T|D) = 0.9（患病者检出率）
• P(T|Dc) = 0.01（未患病者假阳性率）
• P(Dc) = 1 - 0.0005 = 0.9995

计算得：

分子 = 0.9 * 0.0005 = 0.00045
分母 = 0.00045 + 0.01 * 0.9995 = 0.010445
P(D|T) ≈ 0.00045 / 0.010445 ≈ 0.0431（约4.31%）

2. P(D|Tc)（阴性时患病的概率） #

同样用贝叶斯公式：

P(D|Tc) = [P(Tc|D)*P(D)] / [P(Tc|D)*P(D) + P(Tc|Dc)*P(Dc)]

代入数值：

• P(Tc|D) = 1 - 0.9 = 0.1（患病者未检出率）
• P(Tc|Dc) = 1 - 0.01 = 0.99（未患病者真阴性率）

计算得：

分子 = 0.1 * 0.0005 = 0.00005
分母 = 0.00005 + 0.99 * 0.9995 = 0.989555
P(D|Tc) ≈ 0.00005 / 0.989555 ≈ 0.0000505（约0.00505%）

你的代码中的问题 #

现在来看你的R模拟代码，有几个关键错误导致结果不正确：

• 未处理未患病样本的测试结果：你只模拟了患病（D=1）的情况，完全忽略了未患病（D=0）时的假阳性（T=1）和真阴性（Tc=1）情况，这会导致T和Tc的样本量严重失真。
• 患病时的测试逻辑错误：你用了两个独立的随机抽样来生成T和Tc，这会出现T和Tc同时为1或同时为0的不合理情况——测试结果只能是阳性或阴性，二者互斥。正确的做法是用一次抽样决定测试结果：患病时90%概率阳性，10%概率阴性。
• 样本量太小：患病率只有0.05%，1000个样本里可能只有0或1个患病者，统计结果的误差会极大，需要增大样本量（比如100万）才能得到稳定的模拟结果。

修正后的R模拟代码 #

set.seed(110)
sims = 1e6  # 增大样本量到100万
D = numeric(sims)
T_result = numeric(sims)  # 避免用T做变量名，T在R里是逻辑真的缩写

# 第一步：模拟是否患病
D = rbinom(sims, size = 1, prob = 0.0005)

# 第二步：根据患病情况模拟测试结果
for(i in 1:sims){
  if(D[i] == 1){
    # 患病：90%阳性，10%阴性
    T_result[i] = rbinom(1, size = 1, prob = 0.9)
  } else {
    # 未患病：1%阳性，99%阴性
    T_result[i] = rbinom(1, size = 1, prob = 0.01)
  }
}

# 计算P(D|T)：测试阳性时患病的概率
positive_cases = which(T_result == 1)
P_D_given_T = mean(D[positive_cases])
cat("P(D|T)模拟值：", P_D_given_T, "\n")

# 计算P(D|Tc)：测试阴性时患病的概率
negative_cases = which(T_result == 0)
P_D_given_Tc = mean(D[negative_cases])
cat("P(D|Tc)模拟值：", P_D_given_Tc, "\n")

运行这段代码后，你会得到接近理论值的结果：

• P(D|T)大约在0.043左右
• P(D|Tc)大约在0.00005左右

另外，也可以不用循环，用向量化操作让代码更高效：

set.seed(110)
sims = 1e6
D = rbinom(sims, 1, 0.0005)
# 生成测试结果：根据D的值选择对应的概率
T_result = rbinom(sims, 1, ifelse(D == 1, 0.9, 0.01))

# 计算概率
P_D_given_T = mean(D[T_result == 1])
P_D_given_Tc = mean(D[T_result == 0])

cat("P(D|T)：", round(P_D_given_T, 4), "\n")
cat("P(D|Tc)：", round(P_D_given_Tc, 6), "\n")

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