基于Python的贝叶斯估计：部件已失败k次后全败概率及αβ参数求解

问题背景 #

生产的部件需经过1-3轮质检测试：单次测试失败概率p未知（给定p时测试相互独立），只要通过一次即判定合格，三次全败则判定不合格。需利用过往测试历史数据估计p的分布，计算当一个部件已失败k次时，最终三次全败的概率。

现有尝试代码 #

# 计算每个softbin的先验通过率
softbins_total_counts_df = softbins_df.withColumn('total', F.col('total_pass') + F.col('total_fail'))
prior_pass_prob = F.col('total_pass') / F.col('total')
prior_fail_prob = F.col('total_fail') / F.col('total')

# 基于经验先验的贝叶斯估计
rho = 0.3
prior_window = Window.partitionBy('softbin_first_test').orderBy(F.col('time_window')).rowsBetween(Window.unboundedPreceding, Window.currentRow)
prior_pass = F.mean(prior_pass_prob).over(prior_window)
prior_fail = F.mean(prior_fail_prob).over(prior_window)
alpha = prior_pass * (1 - rho) / rho
beta = prior_fail * (1 - rho) / rho
pass_prob = alpha / (alpha + beta)

当前问题 #

运行后发现prior_pass与pass_prob数值完全相同，原因是：
prior_pass + prior_fail = 1，因此alpha + beta = (prior_pass + prior_fail) * (1 - rho)/rho = (1 - rho)/rho，最终pass_prob = alpha/(alpha+beta) = prior_pass，完全无法体现贝叶斯估计对先验的修正效果。

调整方案 #

核心问题拆解 #

1. 原代码仅计算了先验Beta分布的参数，未结合当前部件的观测数据（已失败k次），因此后验概率等于先验概率；
2. 先验参数的计算逻辑导致Beta分布的伪计数强度（α+β）仅由rho决定，未赋予先验合理的信息量权重。

修正步骤 #

1. 合理设置先验Beta参数

Beta分布Beta(α, β)的均值对应先验通过率prior_pass，需引入先验伪计数强度n0（代表先验信息等价于n0个历史样本的置信度，比如n0=100，可根据业务经验调整）：

# 替换原alpha、beta计算逻辑
n0 = 100  # 先验伪计数强度，值越大先验权重越高
alpha = prior_pass * n0
beta = prior_fail * n0

此时α + β = n0，远大于1，先验具备足够的信息量权重。

2. 结合当前观测更新后验分布

当部件已失败k次时，根据Beta分布的共轭特性，后验参数更新为：

# 已失败k次，无成功记录，更新后验参数
alpha_post = alpha  # 成功计数不变
beta_post = beta + k  # 失败计数增加k次

3. 计算三次全败的概率

三次全败要求部件在已失败k次的基础上，再连续失败3 - k次。对于后验分布Beta(α_post, β_post)，p的m阶矩（即E[p^m]）可直接计算，这里m=3−k：

from pyspark.sql import functions as F

# 定义计算E[p^m]的函数
def calculate_p_moment(alpha_col, beta_col, m):
    numerator = F.product([alpha_col + i for i in range(m)])
    denominator = F.product([alpha_col + beta_col + i for i in range(m)])
    return numerator / denominator

调整后完整代码 #

from pyspark.sql import functions as F
from pyspark.sql.window import Window

# 计算每个softbin的先验通过率
softbins_total_counts_df = softbins_df.withColumn('total', F.col('total_pass') + F.col('total_fail'))
prior_pass_prob = F.col('total_pass') / F.col('total')
prior_fail_prob = F.col('total_fail') / F.col('total')

# 定义窗口计算滚动先验均值
prior_window = Window.partitionBy('softbin_first_test').orderBy(F.col('time_window')).rowsBetween(Window.unboundedPreceding, Window.currentRow)
prior_pass = F.mean(prior_pass_prob).over(prior_window)
prior_fail = F.mean(prior_fail_prob).over(prior_window)

# 设置先验伪计数强度，生成Beta先验参数
n0 = 100
alpha = prior_pass * n0
beta = prior_fail * n0

# 假设当前部件已失败k次（示例k=1，可替换为实际值）
k = 1
alpha_post = alpha
beta_post = beta + k

# 计算三次全败概率（需再失败3-k次）
m = 3 - k
full_fail_prob = calculate_p_moment(F.col('alpha_post'), F.col('beta_post'), m)

# 应用到DataFrame
result_df = softbins_total_counts_df.withColumn('full_fail_prob', full_fail_prob)

关键说明 #

• 先验强度n0：n0越大，先验信息权重越高，后验结果越贴近历史平均；n0越小，当前观测的影响越大，后验结果更反映单个部件的异常。
• 共轭先验优势：Beta分布作为伯努利试验的共轭先验，无需复杂积分即可直接更新后验参数，计算效率极高。

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