# 教材位置：第9章，第9.9节；图 9-5。
# 单位、变量与形状：X、Xc、scores 为 (300,2)；特征单位统一，得分沿用该单位，方差比例无量
# 纲。
# 案例 31：从中心化、样本协方差和 奇异值分解 三个角度理解 主成分分析。
# 学习目标：提取二维相关数据中方差最大的正交方向，理解主成分得分与解释方差比例。
# 核心关系：Xc=X-列均值，C=XcᵀXc/(N-1)；若 Xc=USVᵀ，
# 则 C=V diag(s²/(N-1)) Vᵀ，scores=XcV=US，不同得分列之间的样本协方差为零。
# 数据与形状：N=300 个样本、2 个特征；X/Xc/scores 均为 (300,2)，C 为 (2,2)。
# 本例只中心化，没有按标准差标准化；主成分分析 方向因此依赖原特征的量纲与尺度。
# 输出：打印两个主成分的方差比例，验证等价关系，生成 e31_pca.png。
# 运行方式：安装依赖后在本目录执行 python e31_pca.py，图片保存到当前工作目录。

# NumPy 用于模拟相关变量、样本协方差、对称特征分解与 奇异值分解。
import numpy as np
# 左图展示中心化数据和主轴，中图展示旋转到主轴坐标系后的得分，右图展示解释率。
import matplotlib.pyplot as plt
# 在创建第一幅图之前配置可跨系统选择的中文字体。
from plot_style import configure_chinese_plot

configure_chinese_plot()
# 固定种子 5，便于复现相同的点云和解释方差比例。
rng = np.random.default_rng(5)
# 生成 300 个共同潜在值，均值 0.5、标准差 0.18；两个特征共用它来制造相关性。
base = rng.normal(0.5, 0.18, 300)
# 第一特征为 base 加标准差 0.05 的独立噪声，第二特征为 0.8*base+0.25 加噪声。
# column_stack 把两个长度 300 的数组按列拼成 (300,2)：行是样本，列是特征。
X = np.column_stack([base + rng.normal(0, 0.05, 300),
                     0.8 * base + 0.25 + rng.normal(0, 0.06, 300)])
# axis=0 沿样本行求各列均值，形状 (2,) 自动广播到 300 行。
# 先移除均值才能让 主成分分析 描述围绕样本中心的变化；未中心化的 奇异值分解 会混入均值位
# 置。
Xc = X - X.mean(axis=0)
N = len(X)
# (2,300)@(300,2) 得到 2×2 散布矩阵；除以 N-1 采用样本协方差口径。
# 对角线是各特征样本方差，非对角线是两特征样本协方差。
C = Xc.T @ Xc / (N - 1)
# C 为实对称矩阵，eigh 返回升序特征值，directions 的各列为相应的标准正交方向。
values, directions = np.linalg.eigh(C)
# 同时反转特征值和特征向量列，使第一列对应最大方差方向；保持配对关系不能只反转一方。
values = values[::-1]
directions = directions[:, ::-1]
# 中心化矩阵的精简 奇异值分解 形状为 U(300,2)、s(2,)、Vh(2,2)，奇异值已按降序排列。
U, s, Vh = np.linalg.svd(Xc, full_matrices=False)
# V 的列为特征空间里的主轴；Xc@V 把每个样本投影到这些轴上，得到主成分得分。
# 每轴方差为 sⱼ²/(N-1)，求比例时公因子 N-1 抵消，因此 ratio=s²/∑s²。
V = Vh.T
scores = Xc @ V
ratio = s ** 2 / (s @ s)
# 两个比例之和为 1，表示保留各主成分能解释多少中心化数据的总样本方差。
print("主成分方差比例 =", np.round(ratio, 6))
# 核对协方差特征值与中心化数据奇异值平方之间的关系，分母保持 N-1 一致。
assert np.allclose(values, s ** 2 / (N - 1))
# 得分均值仍为零，因此该式正是得分样本协方差矩阵。
# 非对角元接近零说明主成分线性不相关；这通常不意味着统计独立。
assert np.allclose(scores.T @ scores / (N - 1), np.diag(values))
# 本例存在明显的最大方差主方向；特征值不等也使该方向除正负号外唯一。
assert values[0] > values[1]
# 特征分解与 奇异值分解 可能为同一条主轴返回相反方向，直接比较向量会误判。
# 比较 vvᵀ 投影矩阵可以消除符号差异，因为 (-v)(-v)ᵀ=vvᵀ。
assert np.allclose(np.outer(directions[:, 0], directions[:, 0]),
                   np.outer(V[:, 0], V[:, 0]))
# XcV=US；U*s 用广播把每列 U 乘以对应奇异值，形状仍为 (300,2)。
assert np.allclose(scores, U * s)
# 创建中心化原坐标、主成分得分与解释率三幅图；全都使用上面同一批300个样本。
fig, axes = plt.subplots(1, 3, figsize=(6.3, 3.0))
# Xc.T 的两行被展开为散点横纵坐标；小点和半透明可减轻重叠遮挡。
axes[0].scatter(*Xc.T, s=8, alpha=0.4)
# 分别绘制 PC1 与 PC2；两条方向线应正交。
for j in range(2):
    # 特征值是该轴方差，其平方根为标准差；取两倍标准差作为每一侧示意轴长度。
    v = 2 * np.sqrt(values[j]) * V[:, j]
    # 从 -v 画到 +v，线段经过中心化后的原点；它只是主轴示意，并非置信区间。
    axes[0].plot([-v[0], v[0]], [-v[1], v[1]], label=f"主成分 {j + 1}")
# 左图横纵单位等长，保证图上的主轴角度忠实反映几何正交关系。
axes[0].set_aspect("equal")
axes[0].legend(fontsize=7.8)
# 左图数据中心在原点，展示中心化后的原特征坐标。
axes[0].set(xlabel="中心化特征 1", ylabel="中心化特征 2", title="中心化数据")
# 中图以 PC1、PC2 得分为横纵坐标，点云的主要变化沿 PC1 方向展开。
axes[1].scatter(*scores.T, s=8, alpha=0.4)
axes[1].set(xlabel="主成分 1 得分", ylabel="主成分 2 得分", title="主成分得分")
# 得分图也采用等比例坐标，避免纵横自动缩放改变读者对主次变化的观感。
axes[1].set_aspect("equal")
# 右图直接读取已算出的ratio；不重新生成点云或用另一套特征值。
axes[2].bar([1, 2], ratio, width=0.55, color=["tab:blue", "tab:orange"])
axes[2].set_xticks([1, 2], ["主成分 1", "主成分 2"])
axes[2].set(ylabel="方差比例", ylim=(0, 1.1), title="解释方差比例")
for index, value in enumerate(ratio, 1):
    axes[2].text(index, value + 0.035, f"{value:.2%}", ha="center", fontsize=8)
# 为各幅图添加淡网格，帮助比较点云位置、主轴方向与方差份额。
for ax in axes:
    ax.grid(alpha=0.3)
    ax.tick_params(labelsize=7)
# 调整布局、保存图片并关闭图对象，以支持脚本批量执行。
fig.tight_layout()
fig.savefig("e31_pca.png", dpi=180, bbox_inches="tight")
plt.close(fig)