# 教材位置：第3章，第3.3节；图 3-2。
# 单位、变量与形状：k、target、approx 为 (30,)；横轴为计数，概率质量和近似误差无量纲。
# 学习目标：观察二项分布在“试验次数增大、单次概率减小、均值固定”时的泊松极限。
# 核心关系：固定 λ，令 p=λ/n，当 n→∞ 时 二项分布(n,p) 在每个固定整数 k 上的
# 概率趋向 Poisson(λ) 的概率；两者均值均为 λ，但二项分布方差是 λ(1-λ/n)。
# 数据与输出：取 λ=3，比较 n=10、50、500 时 k=0,...,29 上的最大概率差，
# 再输出没有纳入该网格的泊松右尾概率；本例直接计算理论分布。
# 配套图件：运行后另保存 e08_poisson_limit.png，教材直接采用该文件。
# 运行方式：安装 NumPy、SciPy、Matplotlib 后执行 python e08_poisson_limit.py，结果输出到
# 终端。

import numpy as np
from scipy import stats

# lam 表示期望计数 λ；使用 lam 作为变量名，因为 lambda 是 Python 的保留字。
lam = 3.0
# k 是形状 (30,) 的整数数组，比较范围包括 0 和 29；泊松分布的支持集实际上没有有限上界。
k = np.arange(30)
# target[j] 对应 P(Poisson(3)=k[j])；这里得到概率质量，不是概率密度或累积概率。
target = stats.poisson.pmf(k, mu=lam)
# 每轮改变 n，同时让 p 随 n 改变，不能保持 p=0.3 不变后仍声称泊松参数固定为 3。
# 这些 n 均大于 λ，因此 p=λ/n 落在合法的概率区间内。
for n in [10, 50, 500]:
    # 在相同 k 网格计算二项分布；当 k>n 时，其概率为 0，如 n=10 时的 k=11,...,29。
    approx = stats.binom.pmf(k, n=n, p=lam / n)
    # 先逐项取概率差的绝对值，再取最大值；delta 只衡量当前有限网格上的最大点概率差。
    # 它不是全部尾概率的误差，也不是总变差距离，更不是每个 k 上误差的总和。
    delta = np.max(np.abs(approx - target))
    # 格式符只控制对齐和保留小数位，不改变计算所用的 p 或 delta。
    print(f"n={n:3d}, p={lam/n:.4f}，k=0…29 的最大绝对差={delta:.6f}")
    # 验证这三个二项分布的均值都精确对应 n×(λ/n)=λ，仅允许浮点舍入误差。
    assert np.isclose(stats.binom.mean(n, lam / n), lam)
# k[-1]=29；离散分布的 sf(29)=P(X>29)=P(X>=30)，恰好是网格之外的泊松概率质量。
# 输出它是为了交代截断范围；网格内概率之和不足 1 并不代表分布没有正确归一化。
print("网格外遗漏的泊松尾概率：", stats.poisson.sf(k[-1], mu=lam))

# 教材图 3-2复用原诊断网格k=0,...,29，图上最大差与上面打印值使用完全相同的范围。
import matplotlib.pyplot as plt
# 在创建第一幅图之前配置可跨系统选择的中文字体。
from plot_style import configure_chinese_plot

configure_chinese_plot()

fig, axes = plt.subplots(1, 3, figsize=(6.3, 2.8), constrained_layout=True,
                         sharex=True, sharey=True)
# 每个n对应p=lam/n；理论概率质量函数的重复计算不抽样，不改变任何随机状态或终端结果。
for ax, plot_n in zip(axes, [10, 50, 500]):
    plot_p = lam / plot_n
    plot_binomial = stats.binom.pmf(k, n=plot_n, p=plot_p)
    plot_delta = np.max(np.abs(plot_binomial - target))
    ax.bar(k, plot_binomial, alpha=0.55, label='二项分布')
    ax.plot(k, target, 'o-', color='C1', markersize=2.6, linewidth=1, label='泊松分布')
    # 只把主要概率区域放大到0..12；最大差的计算范围仍明确写为0..29。
    ax.set(title=f'n={plot_n}, p={plot_p:g}', xlabel='计数 k',
           xlim=(-0.5, 12.5), ylim=(0, 0.29), xticks=[0, 3, 6, 9, 12])
    ax.text(0.98, 0.95, f'最大差（k=0…29）\n{plot_delta:.6f}',
            transform=ax.transAxes, ha='right', va='top', fontsize=7.8)
axes[0].set_ylabel('概率质量')
axes[0].legend(fontsize=7.8, loc='upper right', bbox_to_anchor=(1, 0.75))
# 文件名明确对应本例；图注需说明显示到12与完整计算到29的区别。
fig.savefig('e08_poisson_limit.png', dpi=180, bbox_inches='tight')
plt.close(fig)