# 教材位置：第4章，第4.4节；图 4-2。
# 单位、变量与形状：μ、σ、threshold 采用同一一般测量单位；z 无量纲，原密度单位为测量单位
# 的倒数。
# 学习目标：理解正态变量的标准化、若干倍标准差的覆盖率及尾部分位点。
# 核心数学关系：X ~ N(mu,sigma^2) 时，Z=(X-mu)/sigma ~ N(0,1)，
# 因此 P(X<=x)=Phi((x-mu)/sigma)，分位数为 mu+sigma*Phi^{-1}(p)。
# 数据与输出：使用均值 170、标准差 8 的假想测量模型，打印阈值 178
# 的累计概率、中心覆盖率，以及单侧/双侧 95% 对应的临界值；保存 e14_normal.png。
# 运行方式：安装 NumPy、SciPy、Matplotlib 后，在终端运行 python e14_normal.py。

# NumPy 用于浮点自检；SciPy 提供正态分布的 累积分布函数 和逆 累积分布函数。
import numpy as np
from scipy import stats

# sigma 表示标准差，须为正；方差为 sigma**2=64。
# threshold 是要评估的原始量纲阈值，与 mu 使用相同单位。
mu, sigma, threshold = 170.0, 8.0, 178.0
# norm 的 loc 是均值，scale 是标准差，不能把方差 64 填入 scale。
rv = stats.norm(loc=mu, scale=sigma)
# 标准化后 z=1，表示阈值位于均值上方一个标准差处。
z = (threshold - mu) / sigma
# 累积分布函数 返回左侧累计概率 P(X<=178)，不是阈值处的密度。
print("标准分数 z、累积概率：", z, rv.cdf(threshold))
# stats.norm 不指定参数时是标准正态分布；验证两种计算途径相同。
assert np.isclose(rv.cdf(threshold), stats.norm.cdf(z))
# width 以标准差为单位；分别计算中心区间 mu±sigma、mu±2sigma、mu±3sigma。
for width in [1, 2, 3]:
    # 累积分布函数 相减给出 P(-width<=Z<=width)，即中心覆盖率，而非单侧概率。
    coverage = stats.norm.cdf(width) - stats.norm.cdf(-width)
    # :.6% 将 0 到 1 的概率按百分数输出，并保留六位小数。
    print(f"均值两侧 {width} 个标准差内：{coverage:.6%}")
# ppf 是 累积分布函数 的反函数；0.95 分位点右侧还剩 5% 的概率质量。
# 此处调用带参数的 rv，所以结果仍使用原始测量单位。
print("右尾概率 5% 的阈值：", rv.ppf(0.95))
# 单侧 95% 界把全部 5% 留在一侧，临界值为 Phi^{-1}(0.95)。
# 双侧中心 95% 区间把 5% 均分为两尾，右临界值是 Phi^{-1}(0.975)，
# 左临界值为其相反数；两种临界值都处于标准化后的 z 尺度。
print("单侧 95%、双侧 95% 临界值：",
      stats.norm.ppf(0.95), stats.norm.ppf(0.975))

# 教材图 4-2把原尺度模型与标准化后的模型并列；主参数仍是mu=170、sigma=8、threshold=178。
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)
# 原尺度密度网格覆盖均值左右四个标准差；网格只作显示，不改变总体参数。
normal_x = np.linspace(mu - 4 * sigma, mu + 4 * sigma, 600)
normal_pdf = rv.pdf(normal_x)
axes[0].plot(normal_x, normal_pdf)
# 阴影是阈值左侧的累计概率；阈值线与上面计算的z=1对应。
axes[0].fill_between(normal_x, normal_pdf, where=normal_x <= threshold,
                     alpha=0.25)
axes[0].axvline(threshold, color='C1', linestyle='--',
                label=f'P(X <= {threshold:g}) = {rv.cdf(threshold):.6f}')
axes[0].set(title=f'正态：均值={mu:g}，标准差={sigma:g}',
            xlabel='原尺度数值 x', ylabel='概率密度')
axes[0].legend(fontsize=7.8, loc='upper left')
# 右图的横轴是无量纲z；同一幅标准正态密度标出三个中心覆盖区间。
standard_z = np.linspace(-4, 4, 801)
standard_pdf = stats.norm.pdf(standard_z)
for band_width, band_color in [(3, '#C9DFEB'), (2, '#86B5CE'), (1, '#E7A067')]:
    band_coverage = stats.norm.cdf(band_width) - stats.norm.cdf(-band_width)
    axes[1].fill_between(standard_z, standard_pdf,
                         where=np.abs(standard_z) <= band_width,
                         color=band_color, alpha=0.9,
                         label=f'+/-{band_width} 个标准差内：{band_coverage:.2%}')
axes[1].plot(standard_z, standard_pdf, color='C0', linewidth=1)
axes[1].axvline(z, color='C1', linestyle='--', linewidth=1)
axes[1].set(title=r'标准化 $Z=(X-\mu)/\sigma$',
            xlabel='标准分数 z', ylabel='概率密度')
axes[1].legend(fontsize=7.8, loc='upper right')
fig.savefig('e14_normal.png', dpi=180, bbox_inches='tight')
plt.close(fig)