R用户求助：如何用OpenBUGS计算设备平均寿命>100小时的后验概率？

问题背景 #

已知设备寿命服从参数为theta的指数分布（OpenBUGS中dexp(theta)为率参数形式，均值为1/theta，单位：千小时），theta的先验分布为Gamma(alpha=45, beta=5)（OpenBUGS中dgamma(shape, rate)，对应shape=45，rate=5）。现有数据向量：
Y = c(0.052, 0.105, 0.234, 0.031, 0.061, 0.243, 0.118, 0.173, 0.042, 0.007)（单位：千小时），需要计算设备平均寿命大于100小时（即0.1千小时）的后验概率。

现有代码的核心错误 #

你当前代码中计算的是prob_theta_less_than_0.01，这是完全错误的：

• 平均寿命（千小时）= 1/theta，要满足平均寿命>0.1千小时（100小时），等价于1/theta > 0.1 → theta < 10，而非theta < 0.01。
• 你的先验Gamma(45,5)均值为45/5=9，数据样本均值为0.1066千小时，对应theta的样本估计约为9.38，远大于0.01，因此prob_theta_less_than_0.01的后验均值为0是合理的，但这不是你需要的概率。

修正后的OpenBUGS代码 #

model {
  # 似然函数：数据服从指数分布（率参数theta）
  for (i in 1:10) {
    Y[i] ~ dexp(theta)
  }
  
  # 先验分布：theta ~ Gamma(shape=45, rate=5)
  theta ~ dgamma(45, 5)
  
  # 计算平均寿命>0.1千小时（100小时）的后验概率
  # 等价于 theta < 10，用step函数指示：step(a)返回1当a>0，否则0
  prob_mean_gt_100h <- step(10 - theta)
}

# 数据列表
list(Y = c(0.052, 0.105, 0.234, 0.031, 0.061, 0.243, 0.118, 0.173, 0.042, 0.007))

OpenBUGS GUI操作步骤 #

• 打开OpenBUGS，点击左上角File → New，新建空白窗口，粘贴上述模型代码（从model {到}部分），保存为.bug文件（比如lifetime_model.bug）。
• 加载数据：点击File → Load Data，选择包含数据列表的文本文件（或直接在新建窗口粘贴数据列表部分，保存为.txt后加载）。
• 编译模型：点击Model → Compile，弹出窗口输入数据长度10，点击OK，底部状态栏显示model compiled即成功。
• 初始化参数：点击Model → Initialize，选择Generate inits，点击OK，状态栏显示model initialized即成功。
• 设置抽样参数：点击Inference → Samples，在弹出窗口输入要追踪的变量：theta和prob_mean_gt_100h，点击Set。
• 开始抽样：点击Model → Update，输入迭代次数（比如10000次，其中前2000次为燃烧期），点击OK。
• 查看结果：点击Inference → Summary，选择prob_mean_gt_100h，即可看到该变量的后验均值（即你需要的后验概率），同时可查看theta的后验分布验证合理性。

结果解释 #

prob_mean_gt_100h的后验均值就是平均寿命大于100小时的后验概率。基于你的先验和数据，这个概率应该接近1（因为theta的后验分布集中在9-10附近，大部分样本满足theta<10）。

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