# 教材位置：第7章，第7.10节；图 7-5。
# 从四个二维观测构造样本协方差，再核对线性组合的样本方差。
# X 形状 (4,2)：行是观测，列是特征；这里两列均用无量纲教学数值。
# mean 形状 (2,)，Xc 形状 (4,2)，C 形状 (2,2)，a 形状 (2,)。
# 核心关系：C=Xc.T@Xc/(n−1)，var(X@a,ddof=1)=a.T@C@a。
# 输入固定，不随机抽样；半正定性来自任何方向方差均非负。
# 输出中心、中心化平方和、样本协方差、特征值及两条方差计算路径。
# 执行 python e38_covariance_matrix.py，图片写入当前工作目录。

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

X = np.array([[1, 2], [2, 1], [3, 5], [4, 4]], dtype=float)
n = X.shape[0]
# axis=0 消去观测这一维，每列保留一个均值。
mean = X.mean(axis=0)
# 广播把两个列均值复制到四行，不对行向量逐行分别中心化。
Xc = X - mean
cross_products = Xc.T @ Xc
C = cross_products / (n - 1)
# 手算：列均值 (2.5,3)，平方和为 5、10，交叉乘积和为 5。
assert np.allclose(mean, [2.5, 3.0])
assert np.allclose(cross_products, [[5, 5], [5, 10]])
manual_cov = np.array([[5 / 3, 5 / 3], [5 / 3, 10 / 3]])
assert np.allclose(C, manual_cov)
# np.cov 默认把行视为变量；rowvar=False 才与本例“列为特征”一致。
assert np.allclose(C, np.cov(X, rowvar=False, ddof=1))

a = np.array([0.6, 0.4])
combined = X @ a
# combined=[1.4,1.6,3.8,4]；它的样本均值为 2.7。
direct_var = combined.var(ddof=1)
quadratic_var = a @ C @ a
manual_var = 0.6**2 * (5 / 3) + 2 * 0.6 * 0.4 * (5 / 3)
manual_var += 0.4**2 * (10 / 3)
assert np.allclose(combined, [1.4, 1.6, 3.8, 4.0])
assert np.isclose(direct_var, quadratic_var)
assert np.isclose(quadratic_var, manual_var)
# C 对称，可用 eigvalsh；浮点核验允许很小的负舍入误差。
eigenvalues = np.linalg.eigvalsh(C)
assert np.allclose(C, C.T)
assert eigenvalues.min() >= -1e-12
# 这里特征值均正，故四个观测在二维中并非完全落在同一条直线上。
print("列均值：", mean)
print("中心化后的矩阵：\n", Xc)
print("交叉乘积和矩阵：\n", cross_products)
print("样本协方差矩阵：\n", C)
print("协方差矩阵的特征值：", eigenvalues)
print("线性组合观测及均值：", combined, combined.mean())
print("直接样本方差、二次型方差、手算方差：",
      direct_var, quadratic_var, manual_var)

# 左图保留原始观测和中心，右图逐元素显示 C，帮助对应三个手算和。
configure_chinese_plot()
fig, axes = plt.subplots(1, 2, figsize=(6.3, 3.1), constrained_layout=True)
axes[0].scatter(X[:, 0], X[:, 1], label="四个观测")
axes[0].scatter(*mean, marker="x", s=70, color="red", label="样本中心")
for index, point in enumerate(X, 1):
    axes[0].annotate(str(index), point, xytext=(4, 4), textcoords="offset points")
axes[0].quiver(*mean, *a, angles="xy", scale_units="xy", scale=1, color="gray")
axes[0].set(title="观测与线性组合方向 a", xlabel="第一特征", ylabel="第二特征")
axes[0].set_aspect("equal")
axes[0].legend(fontsize=7.5)
axes[0].grid(alpha=0.2)
axes[1].imshow(C, cmap="Blues", vmin=0, vmax=C.max() * 1.3)
for i in range(2):
    for j in range(2):
        axes[1].text(j, i, f"{C[i, j]:.4f}", ha="center", va="center")
axes[1].set(xticks=[0, 1], yticks=[0, 1],
            xticklabels=["第一特征", "第二特征"],
            yticklabels=["第一特征", "第二特征"], title="样本协方差 C")
fig.savefig("covariance_matrix.png", dpi=180, bbox_inches="tight")
plt.close(fig)