# 教材位置：第1章，第1.6节；图 1-1。
# 单位、变量与形状：prior、posterior 为 (4,)；概率无量纲，图中乘100表示百分比，期望数量
# 按件计。
# 学习目标：理解阳性检测结果怎样更新次品概率，以及为什么基准次品率会影响后验概率。
# 核心关系：设 D 为次品、+ 为阳性，则 P(D|+)=P(+|D)P(D)/P(+)；
# P(+)=P(+|D)P(D)+P(+|非D)P(非D)，后验优势=先验优势×阳性似然比。
# 数据与输出：对四种假定次品率使用同一组灵敏度、特异度，逐行打印先验和阳性后的后验。
# 这些参数用于教学演示，结果是给定模型下的理论概率，不是从产品抽检数据估计出的产品质量结
# 论。
# 运行方式：安装 NumPy、Matplotlib 后执行 python e06_bayes.py；输出终端结果和教材图
# e06_bayes.png。

import numpy as np

# prior 形状为 (4,)，各元素是不同情景的 P(D)，并非四个互斥类别的概率。
# 因此这四个数不需要相加为 1；数组运算让四种情景一次完成相同的公式计算。
prior = np.array([0.001, 0.01, 0.05, 0.10])
# 灵敏度=P(+|D)，特异度=P(-|非D)；二者都不能直接解释成 P(D|+)。
sensitivity, specificity = 0.99, 0.95
# 乘法公式得到联合概率 P(D 且 +)，不是只针对次品群体内部的条件比例。
true_positive = sensitivity * prior
# 1-specificity 为假阳性率 P(+|非D)；再乘合格品比例，得到 P(非D 且 +)。
false_positive = (1 - specificity) * (1 - prior)
# 分母把所有阳性来源相加，得到 P(+)；分子占其中的比例才是阳性后的次品概率。
# 即使检测参数相同，先验很低时假阳性仍可能占阳性样本中的较大部分。
posterior = true_positive / (true_positive + false_positive)
# 优势=p/(1-p) 比较“发生”与“不发生”的概率，不是概率本身，允许大于 1。
odds_prior = prior / (1 - prior)
# 阳性似然比比较阳性结果在次品与合格品两种条件下的可能性，不是后验概率。
ratio = sensitivity / (1 - specificity)
# 贝叶斯公式的优势写法将证据更新化为乘法；标量 ratio 会广播到四个情景。
odds_post = odds_prior * ratio
# 用 p=优势/(1+优势) 转回概率，验证概率写法和优势写法逐元素一致。
# allclose 允许浮点舍入误差；这里只检验两种代数形式，没有进行随机模拟。
assert np.allclose(posterior, odds_post / (1 + odds_post))
# zip 按对应位置配对先验、后验；百分比格式只影响显示，不改变内部计算使用的小数。
for pre, post in zip(prior, posterior):
    print(f"先验={pre:.1%}，后验={post:.2%}")

# 教材图 1-1只采用上面已经计算的四种先验情景；不增加新的检测参数或随机样本。
import matplotlib.pyplot as plt
# 在创建第一幅图之前配置可跨系统选择的中文字体。
from plot_style import configure_chinese_plot

configure_chinese_plot()

# 左图比较同一检测规则下的先验和阳性后验；乘100只改变显示单位为百分数。
fig, axes = plt.subplots(1, 2, figsize=(6.3, 2.9), constrained_layout=True)
positions = np.arange(prior.size)
axes[0].bar(positions - 0.18, 100 * prior, width=0.36, label='先验 P(D)')
axes[0].bar(positions + 0.18, 100 * posterior, width=0.36, label='后验 P(D | +)')
axes[0].set(xticks=positions, xticklabels=[f'{value:.1%}' for value in prior],
            xlabel='先验次品率', ylabel='概率（%）',
            title='四种先验情景', ylim=(0, 85))
# 后验百分数直接来自posterior；避免在图中硬写与计算结果分离的数字。
for position, value in zip(positions, posterior):
    axes[0].text(position + 0.18, 100 * value + 1.5, f'{value:.2%}',
                 ha='center', fontsize=7.8)
axes[0].legend(fontsize=7.8, loc='upper left')
# 右图只把第一种先验情景换算为10万件产品的期望数量，未把概率改成实际观测频数。
population_size = 100000
positive_counts = population_size * np.array([true_positive[0], false_positive[0]])
axes[1].bar(['真阳性', '假阳性'], positive_counts,
            color=['C1', 'C0'], width=0.55)
for position, value in enumerate(positive_counts):
    axes[1].text(position, value + 100, f'{value:.0f}', ha='center', fontsize=8)
axes[1].set(ylabel='产品期望件数', ylim=(0, 6600),
            title=f'{population_size:,} 件；先验={prior[0]:.1%}')
# 图中的阳性件数分母与上面posterior[0]使用同一个全概率分母。
axes[1].text(0.5, 0.95,
             f'P(D | +) = {positive_counts[0]:.0f} / {positive_counts.sum():.0f}'
             f' = {posterior[0]:.2%}',
             transform=axes[1].transAxes, ha='center', fontsize=7.8)
fig.savefig('e06_bayes.png', dpi=180, bbox_inches='tight')
plt.close(fig)