# 教材位置：第3章，第3.3节；仅终端输出。
# 单位、变量与形状：sample 为 (50000,)；观测为单位暴露量中的事件次数，概率无量纲。
# 学习目标：核对泊松分布的均值、方差和二阶阶乘矩，理解无限求和的截断与模拟检验。
# 核心关系：X~Poisson(λ) 时 E[X]=Var(X)=λ，E[X(X-1)]=λ²；
# 理论期望本应对 k=0,1,... 的全部可能值加权求和，代码使用高分位数进行有限截断。
# 数据与输出：取 λ=5，输出精确理论矩、截断求和的矩、遗漏尾概率，以及 50000 次模拟的样本
# 矩。
# 运行方式：安装 NumPy、SciPy 后执行 python e09_poisson.py；结果输出到终端。

import numpy as np
from scipy import stats

# mu 是 SciPy 泊松分布的均值参数；冻结分布对象保存参数，供后续计算和抽样复用。
rv = stats.poisson(mu=5)
# ppf 是分位数函数：寻找使 P(X<=upper) 至少达到给定概率的整数界限。
# 使用接近 1 的分位数是为了把遗漏概率压得很小；int 转成构造整数网格所需的 Python 整数。
upper = int(rv.ppf(1 - 1e-12))
# arange 不包含上界，因此加 1 使网格覆盖 0,...,upper；这是有限数组，长度为 upper+1。
k = np.arange(upper + 1)
# prob 与 k 位置对应，保存各整数的理论点概率；截断后 prob.sum() 会略小于 1。
prob = rv.pmf(k)
# 这是 ∑_{k=0}^{upper} kP(X=k)，近似 E[X]；没有对截断概率重新归一化。
mean_sum = np.sum(k * prob)
# 这里计算 X(X-1) 的期望，而不是 X² 的期望；泊松分布的理论值为 λ²=25。
# 因为 X²=X(X-1)+X，也可以由这两个矩重建二阶原点矩与方差。
factorial_moment = np.sum(k * (k - 1) * prob)
# 固定随机种子，以复现下面的样本；理论矩与截断计算不使用随机数。
rng = np.random.default_rng(42)
# rvs 生成 (50000,) 形状的泊松整数观测；通过 random_state 显式传入生成器。
sample = rv.rvs(size=50000, random_state=rng)
print("理论均值、方差：", rv.mean(), rv.var())
# 有限和会遗漏尾部的加权贡献；尾概率很小并不自动给出同样大小的矩误差上界，
# 因为均值和阶乘矩的尾项还分别乘上了 k 和 k(k-1)。
print("有限求和的均值、二阶阶乘矩 E[X(X-1)]：", mean_sum, factorial_moment)
# sf(upper)=P(X>upper)，正好对应本次有限网格外的全部概率；这里不是 P(X>=upper)。
print("遗漏的尾概率：", rv.sf(upper))
# 模拟均值和方差会随机波动；ddof=1 用样本量减 1 为分母，估计总体方差。
print("样本均值、样本方差：", sample.mean(), sample.var(ddof=1))
# 首个断言核对理论均值和方差；并未要求有限样本的均值、方差正好等于 5。
assert np.allclose([rv.mean(), rv.var()], [5, 5])
# 网格内的概率加上严格大于 upper 的尾概率应为 1，检查截断端点没有漏算或重复计算。
assert np.isclose(prob.sum() + rv.sf(upper), 1)