# 教材位置：第4章，第4.6节；图 4-4。
# 单位、变量与形状：gaps、arrival_times 为 (5000,)；间隔和时刻按时间计，rate 为每时间单
# 位的事件数。
# 学习目标：从指数分布生成独立事件间隔，累加得到到达时刻，并验证无记忆性。
# 核心数学关系：间隔 T 的速率为 lambda 时，E[T]=1/lambda，
# P(T>t)=exp(-lambda*t)，且 P(T>s+t | T>s)=P(T>t)。
# 数据与输出：模拟 5000 个速率为 0.5 的间隔，打印前五个到达时刻、
# 理论/样本平均间隔和条件生存概率；数值自检检查性质，并保存 e15_exponential.png。
# 运行方式：安装 NumPy、SciPy、Matplotlib 后，在终端运行 python e15_exponential.py。

# NumPy 负责随机模拟和累加；SciPy 用于计算理论分布的性质。
import numpy as np
from scipy import stats

# 创建独立随机数生成器；固定种子便于在相同软件环境下复现教学结果。
rng = np.random.default_rng(15)
# rate 的单位是“每单位时间的平均事件数”；这里平均每单位时间 0.5 次。
rate = 0.5
# SciPy 和 NumPy 的指数分布参数都是 scale，即平均等待时间 1/rate=2。
# rate 与 scale 互为倒数；直接把 rate 传给 scale 会得到另一个分布。
rv = stats.expon(scale=1 / rate)
# size=5000 返回形状 (5000,) 的一维数组，每个元素是一段独立等待时间。
gaps = rng.exponential(scale=1 / rate, size=5000)
# 从时刻 0 开始累计间隔：第 i 个元素是第 i+1 次事件发生的绝对时刻。
# 区分“间隔 gaps”与“时刻 arrival_times”，才能正确解释后续统计量。
arrival_times = np.cumsum(gaps)
# [:5] 切片取前五个元素；样本到达时刻受随机间隔影响，不会等距。
print("前五个到达时刻：", arrival_times[:5])
# rv.mean() 是总体理论均值 2；gaps.mean() 是有限次模拟的样本均值，
# 两者接近但无需完全相等，差异来自蒙特卡洛抽样误差。
print("理论平均间隔、样本平均间隔：", rv.mean(), gaps.mean())
# 已等待 s 个时间单位，考察还需等待超过 t 的条件概率。
s, t = 2.0, 1.0
# sf(u)=P(T>u) 是生存函数。因为 {T>s+t} 包含于 {T>s}，
# 条件概率等于两者概率之比；sf 也避免了尾部用 1-cdf 的相减精度损失。
conditional = rv.sf(s + t) / rv.sf(s)
print("条件生存概率：", conditional)
# 无记忆性说明“已等待多久”不会改变剩余等待时间的分布。
# 这是指数分布的性质，并不适用于所有连续等待时间模型。
assert np.isclose(conditional, rv.sf(t))
# diff 计算相邻到达时刻之差，长度从 5000 减为 4999。
# 指数间隔理论上几乎必然为正；本断言检查本次样本的到达时刻严格递增。
assert np.all(np.diff(arrival_times) > 0)

# 教材图 4-4使用本例已经抽出的gaps和arrival_times，绝不重新设种子或额外抽样。
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)
# 左图只显示5000次到达中的前8次；在最前面补0仅用于画第一个等待间隔。
shown_count = 8
shown_times = np.r_[0.0, arrival_times[:shown_count]]
axes[0].hlines(1, 0, shown_times[-1] + 0.5, color='0.5', linewidth=1)
axes[0].vlines(shown_times[1:], 0.9, 1.1, color='C0', linewidth=1.4)
# 前四个带端点的水平线段展示“间隔”与“累计时刻”的区别，端点完全来自已有数组。
for gap_index in range(4):
    # 很短的相邻间隔可能挤在一起，错开线段和文字的高度即可避免标签重叠。
    arrow_height = 0.75 - 0.17 * (gap_index % 3)
    # 采用线段和端点短线，极短间隔也不会因双箭头尺寸超过间隔而交叉。
    axes[0].hlines(arrow_height, shown_times[gap_index], shown_times[gap_index + 1],
                   color=f'C{gap_index % 3}', linewidth=1.3)
    axes[0].vlines(shown_times[gap_index:gap_index + 2],
                   arrow_height - 0.025, arrow_height + 0.025,
                   color=f'C{gap_index % 3}', linewidth=1)
    axes[0].text((shown_times[gap_index + 1] + shown_times[gap_index]) / 2,
                 arrow_height - 0.11, rf'$T_{{{gap_index + 1}}}$', ha='center',
                     fontsize=7.8)
axes[0].set(title=f'{gaps.size:,} 次到达中的前 {shown_count} 次',
            xlabel='到达时刻', yticks=[], ylim=(0.15, 1.35),
            xlim=(0, shown_times[-1] + 0.5))
# 右图使用全部5000个原间隔，与同一scale=1/rate=2的理论指数密度比较。
# bins=35只控制显示分组；density=True把柱形总面积归一化为1。
axes[1].hist(gaps, bins=35, density=True, alpha=0.5, label='全部模拟间隔')
exponential_t = np.linspace(0, gaps.max(), 500)
axes[1].plot(exponential_t, rv.pdf(exponential_t), color='C1',
             label=f'指数：速率={rate:g}，尺度={1 / rate:g}')
axes[1].set(title='相同的等待时间模型', xlabel='间隔长度 t', ylabel='概率密度')
axes[1].legend(fontsize=7.8)
fig.savefig('e15_exponential.png', dpi=180, bbox_inches='tight')
plt.close(fig)