# 教材位置：第4章，第4.7节；图 4-5。
# 单位、变量与形状：原间隔数组为 (30000,4)，求和后 waits 为 (30000,)；theta 按时间计。
# 学习目标：理解整数形状参数的 伽马分布为何是多个指数等待时间之和，
# 并用 Poisson 计数从另一角度计算第 r 次事件发生的概率。
# 核心数学关系：若 r 个独立间隔均为 Exponential(scale=theta)，
# 则总等待时间 T_r ~ 伽马分布(shape=r, scale=theta)，
# E[T_r]=r*theta、Var(T_r)=r*theta^2，P(T_r<=t)=P(N(t)>=r)。
# 数据与输出：模拟 30000 次“等到第 4 次事件”的过程，比较理论/样本矩，
# 并打印四个时刻前已经发生第 4 次事件的概率，保存 e16_gamma.png。
# 运行方式：安装 NumPy、SciPy、Matplotlib 后，在终端运行 python e16_gamma.py。

import numpy as np
from scipy import stats

# 固定种子使当前环境中的样本可复现；每行模拟过程使用不同随机抽样。
rng = np.random.default_rng(21)
# r 是所等待的事件次数，也是 伽马分布 的形状参数；theta 是每段间隔的
# 平均时间（尺度），对应事件速率 lambda=1/theta=2，不能把尺度当作速率。
r, theta = 4, 0.5
# 先产生形状 (30000,4) 的数组：每行是一轮过程中的四个独立间隔。
# sum(axis=1) 沿列求和，把每行变成总等待时间，结果形状为 (30000,)。
# 形状参数为正整数的 伽马分布也称 Erlang 分布。
waits = rng.exponential(scale=theta, size=(30000, r)).sum(axis=1)
# SciPy 用 a 表示形状参数，用 scale 表示尺度；此处默认 loc=0。
rv = stats.gamma(a=r, scale=theta)
# 在四个非负时刻同时计算概率，数组形状为 (4,)。
t = np.array([0.5, 1.0, 2.0, 3.0])
# 对速率 1/theta 的齐次 Poisson 过程，N(t) 的均值为 t/theta。
# Poisson.sf(r-1)=P(N>r-1)=P(N>=r)，正好对应“第 r 次事件已到达”。
# 离散分布中必须注意这个 r-1；sf(r) 会漏掉 N=r 的情形。
count_formula = stats.poisson.sf(r - 1, mu=t / theta)
# allclose 逐元素比较后汇总为一个布尔值，验证四个时刻的等待/计数关系。
assert np.allclose(rv.cdf(t), count_formula)
# 独立间隔之和的均值、方差分别相加，得到 r*theta 和 r*theta**2。
assert np.allclose([rv.mean(), rv.var()], [r * theta, r * theta ** 2])
# 理论矩由分布参数确定，不受当前模拟样本影响。
print("理论均值、方差：", rv.mean(), rv.var())
# ddof=1 使用分母 30000-1，计算总体方差的无偏样本估计；
# 样本均值和方差存在随机误差，因此不要求它们与理论值完全相等。
print("样本均值、样本方差：", waits.mean(), waits.var(ddof=1))
# cdf(t) 对数组逐元素计算 P(T_r<=t)，输出顺序与 t 中时刻顺序一致。
print("第 r 次事件在 t 前到达的概率：", rv.cdf(t))

# 教材图 4-5采用本例r=4、theta=0.5，与前例rate=0.5的指数模型属于不同参数情景。
# 本图复用已有30000个waits及原t数组；新增理论曲线计算不进行任何随机抽样。
import matplotlib.pyplot as plt
# 在创建第一幅图之前配置可跨系统选择的中文字体。
from plot_style import configure_chinese_plot

configure_chinese_plot()

fig, axes = plt.subplots(1, 2, figsize=(6.3, 2.9), constrained_layout=True)
# 左图把r段指数间隔之和的样本密度与对应伽马分布理论密度比较。
wait_grid = np.linspace(0, max(waits.max(), rv.ppf(0.999)), 600)
axes[0].hist(waits, bins=45, density=True, alpha=0.5, label=f'{waits.size:,} 组间隔和')
axes[0].plot(wait_grid, rv.pdf(wait_grid), color='C1',
             label=f'伽马：形状={r}，尺度={theta:g}')
axes[0].set(title='指数间隔之和', xlabel='总等待时间', ylabel='概率密度')
axes[0].legend(fontsize=7.8)
# 右图两条理论曲线重合：第r次到达已发生，等价于时刻t前计数至少r次。
# count_curve中的Poisson参数为时间/theta；SF传r-1才能包含计数恰好等于r的情况。
cdf_grid = np.linspace(0, rv.ppf(0.999), 400)
count_curve = stats.poisson.sf(r - 1, mu=cdf_grid / theta)
axes[1].plot(cdf_grid, rv.cdf(cdf_grid), linewidth=2.5, alpha=0.6,
    label='伽马累积分布函数')
axes[1].plot(cdf_grid, count_curve, '--', color='C1', linewidth=1.2,
             label=f'泊松尾概率：次数 ≥ {r}')
# 四个实心点就是终端打印的t数组所对应的count_formula，连接图形和数值核验。
axes[1].scatter(t, count_formula, color='C2', s=20, zorder=3, label='核对的时刻 t')
axes[1].set(title=r'$P(T_r\leq t)=P(N(t)\geq r)$',
            xlabel='时间 t', ylabel='累积概率', ylim=(-0.03, 1.05))
axes[1].legend(fontsize=7.8)
fig.savefig('e16_gamma.png', dpi=180, bbox_inches='tight')
plt.close(fig)