# 教材位置：第3章，第3.6节；仅终端输出。
# 单位、变量与形状：sample 为 (50000,)，frequency 为 (21,)；事件频率和标准化差异无量纲。
# 学习目标：用重复模拟核对二项分布的点概率，并用标准误衡量模拟差异的通常尺度。
# 核心关系：独立重复生成 X~二项分布(20,0.3)，令 I=1{X=6}，则 I~Bernoulli(q)，
# q=P(X=6)。count 次重复的样本比例标准误为 sqrt(q(1-q)/count)。
# 数据与输出：生成 50000 个成功次数，打印理论/观察概率、概率估计的标准误与标准化差异，
# 以及理论/观察均值；频数数组用于检查计数方式和归一化是否正确。
# 运行方式：安装 NumPy、SciPy 后执行 python e12_simulation_check.py，结果输出到终端。

import numpy as np
from scipy import stats

# 固定种子，让同一环境下的模拟结果可复现；不同种子产生不同的模拟误差。
rng = np.random.default_rng(11)
# trials=20 是一次二项实验内部的伯努利试验数，count=50000 是整个二项实验的重复次数。
# 概率估计的标准误使用 count；不能误用 trials，因为估计是对 count 个二项观测求比例。
trials, p, count = 20, 0.3, 50000
# sample 的形状为 (50000,)，每个元素是一次独立二项实验中的成功次数，取值为 0,...,20。
sample = rng.binomial(trials, p, size=count)
# bincount 按整数值计数：位置 k 保存 sample 中等于 k 的观测个数。
# minlength=21 保留可能未出现的高端类别；除以 count 后得到经验点概率，和为 1。
# 对离散类别而言这是相对频数，不是连续直方图中的概率密度。
frequency = np.bincount(sample, minlength=trials + 1) / count
# q 由理论 概率质量函数 直接给出，表示“一轮实验恰好有 6 次成功”的概率，不是单次成功概率
# p。
q = stats.binom.pmf(6, trials, p)
# sample==6 形成布尔指示数组；对它求均值就是相应事件的观察比例。
estimate = np.mean(sample == 6)
# 理论 q 在本例已知，因此这里使用的是比例估计量的理论标准误，而不是用 estimate 代入的估
# 计值。
# 它是多次重新生成整个 50000 观测数据集时，estimate 的标准差，不是原始成功次数的标准差。
se = np.sqrt(q * (1 - q) / count)
print("理论概率、模拟频率：", q, estimate)
# (estimate-q)/se 说明本次有符号差异相当于多少个标准误；其有限样本分布并非严格正态。
# 因指示变量独立同分布且方差有限，大样本下可用 中心极限定理 将这个量近似看作标准正态变量
# 。
print("理论标准误、标准化差异：", se, (estimate - q) / se)
# 原始成功次数的理论均值为 trials*p；sample.mean() 是对它的另一项模拟核对。
print("理论均值、样本均值：", trials * p, sample.mean())
# 验证从 0 到 trials 的每个可能成功次数都有一个频率槽位，包括零频数类别。
assert frequency.shape == (trials + 1,)
# 所有观测都恰好归入一个类别，因此经验概率和应为 1，只需允许浮点舍入误差。
assert np.isclose(frequency.sum(), 1)
# 两种统计同一事件的方法应一致；这检查统计实现，没有要求经验概率精确等于理论 q。
# 上述断言通过也不等同于完整的分布拟合检验，仍应结合模拟误差尺度解释理论与观察的差异。
assert np.isclose(frequency[6], estimate)