# 教材位置：第5章，第5.6节；图 5-1。
# 单位、变量与形状：样本按同一一般测量单位计；均值区间保留该单位，t 统计量与 p 值无量纲
# 。
# 学习目标：手算单样本双侧 t 检验和均值置信区间，并与 SciPy 的实现核对。
# 核心数学关系：样本标准差为 s 时，SE=s/sqrt(n)，
# t=(样本均值-mu0)/SE；对正态总体的独立同分布样本，
# 该统计量在原假设下服从 t_(n-1)。
# 前提：总体标准差未知，观测应独立；n=9 很小，精确 t 推断依赖总体正态性，
# 数据若严重偏斜或有异常值，不能仅因使用 t 分布就认为推断可靠。
# 数据与输出：九个假想测量值，检验总体均值是否为 170；打印均值、t 值、
# 双侧 p 值和 95% 置信区间，并保存 e18_t_inference.png。
# 运行方式：安装 NumPy、SciPy、Matplotlib 后，在终端运行 python e18_t_inference.py。

import numpy as np
from scipy import stats

# 一维数组中每个元素代表一个观测，dtype=float 便于后续连续量计算。
x = np.array([168, 172, 175, 169, 171, 173, 170, 174, 167], dtype=float)
# mu0 是原假设 H0: mu=170 的均值；alpha=0.05 是预先选定的显著性水平。
# 此处备择假设是 mu!=170，因此需要考虑正负两个方向的偏离。
mu0, alpha = 170.0, 0.05
# size 返回观测总数；此处 x 是一维，所以也等于数组长度。
n = x.size
# ddof=1 使样本方差分母为 n-1；std 再开平方得到样本标准差 s。
# SE 描述样本均值的抽样波动，和描述个体离散程度的 s 是不同的量。
se = x.std(ddof=1) / np.sqrt(n)
# 标准化均值与原假设的差异；正 t 值表示样本均值大于 mu0。
# 用样本估计标准差后有 n-1=8 个自由度，所以采用 t 分布。
t_value = (x.mean() - mu0) / se
# t 分布关于 0 对称；先算 |t| 右侧尾概率，再乘 2 得到双侧 p 值。
# p 值是在 H0 和模型前提下，得到至少这么极端统计量的概率，
# 并不是“原假设为真”的概率；sf 用于直接计算右尾概率。
p_value = 2 * stats.t.sf(abs(t_value), df=n - 1)
# 双侧中心 95% 对应每尾 alpha/2；ppf 返回右侧正临界值。
critical = stats.t.ppf(1 - alpha / 2, df=n - 1)
# [-1,1] 与标量相乘，一次生成“均值-误差界、均值+误差界”两个端点。
# 95% 置信度是重复抽样时该构造方法覆盖固定总体均值的长期比例，
# 不能把它直接解释成已得到的固定区间对固定参数的后验概率。
interval = x.mean() + np.array([-1, 1]) * critical * se
# SciPy 默认做双侧单样本 t 检验；popmean 指定原假设的总体均值。
result = stats.ttest_1samp(x, popmean=mu0)
print("样本均值、t 统计量、p 值：", x.mean(), t_value, p_value)
print("95% 置信区间：", interval)
# 分别比较统计量、双侧 p 值和区间端点，防止手工推导出现自由度或单双侧错误。
# 浮点运算存在舍入误差，用 isclose/allclose 进行近似比较。
assert np.isclose(t_value, result.statistic)
assert np.isclose(p_value, result.pvalue)
# confidence_interval() 默认置信水平为 0.95，恰好与本例的 1-alpha 相同。
assert np.allclose(interval, result.confidence_interval())

# 教材图 5-1使用原九个观测及本次计算的置信区间；总体均值mu和总体标准差均未知。
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)
# 左图只画本例检验实际采用的n-1自由度t密度，并用标准正态密度作尾部参照。
t_grid = np.linspace(-5, 5, 500)
axes[0].plot(t_grid, stats.t.pdf(t_grid, df=n - 1), label=f't 分布：自由度={n - 1}')
axes[0].plot(t_grid, stats.norm.pdf(t_grid), '--', label='标准正态')
axes[0].axvline(t_value, color='C2', linestyle=':', label=f'观察 t={t_value:.3f}')
axes[0].set(title='原假设下的参考分布', xlabel='t / z', ylabel='概率密度')
axes[0].legend(fontsize=7.8)
# 右图误差条的端点直接取interval，不另写自由度8或分母3，避免数据变化后两处失配。
mean_x = x.mean()
ci_errors = np.array([[mean_x - interval[0]], [interval[1] - mean_x]])
axes[1].errorbar(mean_x, 0.6, xerr=ci_errors, fmt='o', capsize=5)
axes[1].axvline(mu0, color='C1', linestyle='--', label=rf'$H_0: \mu={mu0:g}$')
axes[1].text(mean_x, 0.82, f'样本均值 = {mean_x:.2f}', ha='center', fontsize=8)
axes[1].text(mean_x, 0.34,
             f'{1 - alpha:.0%} 置信区间 = [{interval[0]:.2f}, {interval[1]:.2f}]',
             ha='center', fontsize=7.8)
axes[1].set(title=f'模拟样本：n={n}', xlabel='总体均值尺度',
            yticks=[], ylim=(0.15, 1.05),
            xlim=(min(interval[0], mu0) - 0.7, max(interval[1], mu0) + 0.7))
axes[1].legend(fontsize=7.8, loc='upper right')
fig.savefig('e18_t_inference.png', dpi=180, bbox_inches='tight')
plt.close(fig)