OpenBUGS后验样本生成失败，贝叶斯预测概率计算结果异常求助

问题背景 #

设定的计算假设：

• 两组标准差均为1.44（方差为(1.44^2=2.0736)）
• 组1均值=2.06，组2均值=1.34，预期均值差为0.72
• 组1总样本量(n=80)；组2总样本量(n=40)
• 从两组各借用(n=30)的样本作为先验信息，实际需收集组1 (n=50)、组2 (n=10)的样本

目标是计算贝叶斯预测概率 P(P(Diff>0)>0.95)

当前OpenBUGS代码 #

model{
# likelihoods 
for (i in 1:n){
    x[i]  ~ dnorm(mua,taua)
    y[i]  ~ dnorm(mup,taup)
    diff[i] <- x[i] - y[i]
    difffl[i] <- step(diff[i])
  }
  # normal priors
  mua ~ dnorm(mu1, inv.sigma)
  mup ~ dnorm(mu2, inv.sigma)
  inv.sigma <- n3/sigma.squared
  taua <- n1/sigma.squared
  taup <- n2/sigma.squared
  # calculate the probability P(diff > delta)
  # delta will be 0 here; # eta will be 0.95 here
  # calculate the predictive probability P(P(diff >delta)>0.95)
  sdx <- sd(x[])
  sdy <- sd(y[])
  mdiff <- mean(diff[])
  mdifffl <- mean(difffl[])
  powerfl <- step(mdifffl - 0.95)
}
list(mu1 = 2.06, mu2 = 1.34, n = 1000, n1 = 50, n2 = 10, n3 = 30, 
sigma.squared = 2.0736)

运行结果 #

mean    sd  MC_error    val2.5pc    median  val97.5pc   start   sample
mdiff   0.7152  0.3724  0.003383    -0.01549    0.7143  1.46    1   10000
mdifffl 0.8749  0.1376  0.001318    0.489   0.924   0.999   1   10000
powerfl 0.3932  0.4885  0.005197    0.0 0.0 1.0 1   10000
sdx     0.2036  0.004437    4.722E-5    0.195   0.2036  0.2124  1   10000
sdy     0.4552  0.01011 1.079E-4    0.4357  0.4551  0.4752  1   10000

异常现象 #

• 组1均值的标准差sdx计算异常：理论值应为1.44/sqrt(80)=0.161，但结果为0.2036（对应1.44/sqrt(50)）
• powerfl值仅为0.3932，远低于预期的约93%

错误排查与修正 #

1. 似然部分样本量定义错误 #

代码用同一个循环变量遍历两组样本，且循环上限设为n=1000，完全不符合实际样本量设定：

• 组1实际收集50个样本，组2收集10个，样本量不同，需分开循环
• 循环上限应分别对应n1=50和n2=10

修正后的似然部分：

# likelihoods 
for (i in 1:n1){
    x[i] ~ dnorm(mua, taua)
}
for (j in 1:n2){
    y[j] ~ dnorm(mup, taup)
}

2. 预测概率计算逻辑错误 #

当前代码计算的是模拟样本中diff>0的比例，而非目标的“后验分布中，未来样本Diff>0的概率大于0.95的概率”。

正确逻辑：

1. 从后验分布抽取mua和mup的样本
2. 对每组(mua, mup)，计算未来/总体中Diff>0的概率
3. 统计该概率大于0.95的比例

修正后的计算部分：

# 计算总体均值差
mean_diff <- mua - mup
# 计算总体中Diff>0的概率（用总样本量的均值差标准误）
p_diff_gt0 <- 1 - phi( -mean_diff / sqrt( sigma.squared*(1/(n1+n3) + 1/(n2+n3)) ) )
# 目标概率：P(P(Diff>0)>0.95)
powerfl <- step(p_diff_gt0 - 0.95)

注：phi()是OpenBUGS内置的标准正态累积分布函数

3. 标准差计算错误 #

当前sdx <- sd(x[])计算的是模拟样本的标准差，而非均值的标准误。需修正为：

# 组1均值的标准误（总样本量80）
se_mua <- sqrt( sigma.squared/(n1+n3) )
# 组2均值的标准误（总样本量40）
se_mup <- sqrt( sigma.squared/(n2+n3) )

完整修正代码示例 #

model{
# 似然部分：分别处理两组不同样本量
for (i in 1:n1){
    x[i] ~ dnorm(mua, taua)
}
for (j in 1:n2){
    y[j] ~ dnorm(mup, taup)
}

# 先验分布：基于30个历史样本的信息
mua ~ dnorm(mu1, inv_sigma)
mup ~ dnorm(mu2, inv_sigma)
inv_sigma <- n3 / sigma_squared
taua <- n1 / sigma_squared  # 当前收集的50个样本的精度
taup <- n2 / sigma_squared  # 当前收集的10个样本的精度

# 计算总体均值差
mean_diff <- mua - mup

# 计算总体中Diff>0的概率
p_diff_gt0 <- 1 - phi( -mean_diff / sqrt( sigma.squared*(1/(n1+n3) + 1/(n2+n3)) ) )

# 目标概率：P(P(Diff>0)>0.95)
powerfl <- step(p_diff_gt0 - 0.95)

# 计算均值的标准误（用于验证）
se_mua <- sqrt( sigma.squared/(n1+n3) )
se_mup <- sqrt( sigma.squared/(n2+n3) )
}

list(mu1 = 2.06, mu2 = 1.34, n1 = 50, n2 = 10, n3 = 30, 
sigma_squared = 2.0736)

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