# 教材位置：第10章，第10.8节；图 10-6。
# NIST StRD Norris 的真实臭氧监测仪标定数据，36 对观测原始列顺序为 y,x。
# 来源：https://www.itl.nist.gov/div898/strd/lls/data/Norris.shtml
# 数据：https://www.itl.nist.gov/div898/strd/lls/data/LINKS/DATA/Norris.dat
# 原始资料未声明数值单位；这里不补设单位，也不指定两列为某具体仪器。
# data 形状 (36,2)，A 为 (36,2)，coef 为 (2,)，residual 为 (36,)。
# 设计列顺序 [1,x]，所以 coef[0] 是截距，coef[1] 是斜率。
# 输出手算/最小二乘系数、残差标准差、平方和及 x=500 处拟合值。
# 数据直接嵌入程序，离线可运行；执行 python e42_norris_calibration.py。

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

# 每行严格保留官方文件的 [y,x] 顺序；不能交换两列后仍沿用认证系数。
data = np.array([
    [0.1, 0.2], [338.8, 337.4], [118.1, 118.2],
    [888.0, 884.6], [9.2, 10.1], [228.1, 226.5],
    [668.5, 666.3], [998.5, 996.3], [449.1, 448.6],
    [778.9, 777.0], [559.2, 558.2], [0.3, 0.4],
    [0.1, 0.6], [778.1, 775.5], [668.8, 666.9],
    [339.3, 338.0], [448.9, 447.5], [10.8, 11.6],
    [557.7, 556.0], [228.3, 228.1], [998.0, 995.8],
    [888.8, 887.6], [119.6, 120.2], [0.3, 0.3],
    [0.6, 0.3], [557.6, 556.8], [339.3, 339.1],
    [888.0, 887.2], [998.5, 999.0], [778.9, 779.0],
    [10.2, 11.1], [117.6, 118.3], [228.9, 229.2],
    [668.4, 669.1], [449.2, 448.9], [0.2, 0.5],
], dtype=float)
y, x = data[:, 0], data[:, 1]
n = x.size
A = np.column_stack([np.ones(n), x])
coef, sums, rank, singular = np.linalg.lstsq(A, y, rcond=None)
yhat = A @ coef
residual = y - yhat
# 用中心化求和公式作独立手算核对，而非对 A.T@A 显式求逆。
xc, yc = x - x.mean(), y - y.mean()
manual_slope = (xc @ yc) / (xc @ xc)
manual_intercept = y.mean() - manual_slope * x.mean()
assert n == 36 and rank == 2
assert np.allclose(coef, [manual_intercept, manual_slope])
# 官方认证截距、斜率；只用作数值核验，不反向覆盖本次实际拟合结果。
certified_coef = np.array([-0.262323073774029, 1.00211681802045])
assert np.allclose(coef, certified_coef, rtol=1e-12, atol=1e-11)
sse = residual @ residual
residual_sd = np.sqrt(sse / (n - rank))
r_squared = 1 - sse / (yc @ yc)
# 残差标准差分母为 36−2=34；它不同于直接除以36的均方根误差。
assert np.isclose(sse, 26.6173985294224, rtol=1e-12)
assert np.isclose(residual_sd, 0.884796396144373, rtol=1e-12)
assert np.isclose(r_squared, 0.999993745883712, rtol=1e-12)
assert np.isclose(sums[0], sse)
assert np.allclose(A.T @ residual, 0, atol=1e-8)
# x=500 在已观测范围内；这是拟合均值，不是下一次仪器读数的保证。
x_new = 500.0
prediction = np.array([1.0, x_new]) @ coef
print("数据来源：NIST StRD Norris，36 对臭氧监测仪标定观测")
print("截距、斜率：", coef)
print("手算截距、斜率：", manual_intercept, manual_slope)
print("残差平方和、残差自由度：", sse, n - rank)
print("残差标准差、决定系数：", residual_sd, r_squared)
print("x=500 处的拟合均值：", prediction)
# 高决定系数还受 x 的覆盖范围影响，不能代替残差诊断或外部验证。
# 若要置信/预测区间，还需说明噪声模型和推断条件；本例只做拟合核验。

configure_chinese_plot()
fig, axes = plt.subplots(1, 2, figsize=(6.3, 3.0), constrained_layout=True)
order = np.argsort(x)
axes[0].scatter(x, y, s=14, alpha=0.7, label="官方观测")
axes[0].plot(x[order], yhat[order], color="tab:red", label="最小二乘拟合")
axes[0].set(title="臭氧监测仪标定数据", xlabel="官方数据列 x", ylabel="官方数据列 y")
axes[0].legend(fontsize=8)
axes[1].scatter(x, residual, s=14, alpha=0.7)
axes[1].axhline(0, color="gray", ls="--")
axes[1].set(title="残差需单独检查", xlabel="官方数据列 x", ylabel="观测减拟合值")
for ax in axes:
    ax.grid(alpha=0.2)
fig.savefig("norris_calibration.png", dpi=180, bbox_inches="tight")
plt.close(fig)