# 教材位置：第4章，第4.5节；图 4-3。
# 泊松过程的独立增量与等价间隔法，模拟客服呼入的理想模型。
# 前提是恒定呼入率、独立增量；真实数据的时段变化需另行检验。
# 时间按小时计，lam=3 次/小时，dt=1 小时，每条路径覆盖 24 小时。
# increments 形状 (2000,24)：行是独立路径，列是不重叠的时间窗。
# paths 形状 (2000,25)，包含第 0 时刻；均值、方差数组形状 (24,)。
# 输出理论/模拟矩、相邻增量经验协方差及另一条间隔法路径的事件数。
# 固定种子 2026；执行 python e37_poisson_process.py 保存中文四面板图。

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
from plot_style import configure_chinese_plot

rng = np.random.default_rng(2026)
lam, dt, windows, repetitions = 3.0, 1.0, 24, 2000
# ΔN_j~Poisson(lam*dt)，不同窗及不同路径的抽样相互独立。
increments = rng.poisson(lam * dt, size=(repetitions, windows))
# 沿 axis=1 累积各行的时间窗，最左补一个 N(0)=0。
paths = np.column_stack([np.zeros(repetitions), increments.cumsum(axis=1)])
times = np.arange(windows + 1) * dt
# 沿 axis=0 跨 2000 条独立路径求矩，保留 24 个时间窗。
means = increments.mean(axis=0)
variances = increments.var(axis=0, ddof=1)
adjacent_cov = np.cov(increments[:, 0], increments[:, 1], ddof=1)[0, 1]
# 理论独立增量的协方差为 0；有限模拟的经验协方差通常不恰为 0。
# 不对经验协方差断言等于 0，也不从一次接近 0 推断真实独立性。
assert paths.shape == (repetitions, windows + 1)
assert np.array_equal(paths[:, 1:] - paths[:, :-1], increments)
assert np.all(np.diff(paths, axis=1) >= 0)
print("每窗理论均值、方差：", lam * dt, lam * dt)
print("首窗模拟均值、样本方差：", means[0], variances[0])
print("前两窗增量的经验协方差（理论为 0）：", adjacent_cov)
print("24 小时计数的理论均值、方差：", lam * windows, lam * windows)
print("24 小时计数的模拟均值、样本方差：",
      paths[:, -1].mean(), paths[:, -1].var(ddof=1))

# 同一过程还可用独立 Exp(lam) 间隔生成；尺度参数是 1/lam 小时。
# 逐次累计直至首次越过 24 小时，避免预设个数不足而漏掉事件。
arrival_times = []
clock = 0.0
while True:
    clock += rng.exponential(scale=1 / lam)
    if clock > windows * dt:
        break
    arrival_times.append(clock)
arrival_times = np.asarray(arrival_times)
interval_path = np.searchsorted(arrival_times, times, side="right")
print("间隔法另一条路径的事件总数：", arrival_times.size)
# 两种构造同分布，但使用不同随机抽样，不要求逐点相同。
assert np.all(np.diff(arrival_times) > 0)
assert interval_path[-1] == arrival_times.size

configure_chinese_plot()
fig, axes = plt.subplots(2, 2, figsize=(6.3, 5.1), constrained_layout=True)
for index in range(3):
    axes[0, 0].step(times, paths[index], where="post", label=f"路径 {index + 1}")
axes[0, 0].plot(times, lam * times, "k--", label="理论均值 3t")
axes[0, 0].set(title="独立增量构造", xlabel="时刻（小时）", ylabel="累计呼入次数")
axes[0, 0].legend(fontsize=7)
window_index = np.arange(1, windows + 1)
axes[0, 1].plot(window_index, means, label="逐窗样本均值")
axes[0, 1].plot(window_index, variances, label="逐窗样本方差")
axes[0, 1].axhline(lam * dt, color="gray", ls="--", label="理论值 3")
axes[0, 1].set(title="跨路径估计每窗矩", xlabel="小时窗序号", ylabel="计数矩")
axes[0, 1].legend(fontsize=7)
k = np.arange(12)
frequency = np.bincount(increments[:, 0], minlength=12)[:12] / repetitions
axes[1, 0].bar(k, frequency, alpha=0.5, label="首窗模拟频率")
axes[1, 0].plot(k, stats.poisson.pmf(k, lam * dt), "o-", label="泊松理论概率")
axes[1, 0].set(title="一小时窗的计数分布", xlabel="呼入次数", ylabel="概率质量")
axes[1, 0].legend(fontsize=7)
event_grid = np.r_[0.0, arrival_times, windows * dt]
event_counts = np.r_[0, np.arange(1, arrival_times.size + 1), arrival_times.size]
axes[1, 1].step(event_grid, event_counts, where="post", label="指数间隔法")
axes[1, 1].plot(times, lam * times, "k--", label="理论均值 3t")
axes[1, 1].set(title="另一条连续时间路径", xlabel="时刻（小时）", ylabel="累计呼入次数")
axes[1, 1].legend(fontsize=7)
for ax in axes.flat:
    ax.grid(alpha=0.2)
fig.savefig("poisson_process.png", dpi=180, bbox_inches="tight")
plt.close(fig)