# 教材位置：第10章，第10.7节；图 10-5。
# 用户画像的教学模拟：120 个用户、月订单数、月消费金额、月浏览时长。
# 数据由明确给出的潜在变量和噪声生成，不代表真实用户或任何经营结论。
# 三列单位依次为次、元、分钟；X 形状 (120,3)，每行一个模拟用户。
# 先按样本均值和样本标准差标准化，再用 SVD 求描述性主成分。
# Z 形状 (120,3)，C 为 (3,3)，scores 为 (120,3)，ratio 为 (3,)。
# 固定种子 2026；输出列尺度、解释率、主成分方向与二维重建误差。
# 执行 python e41_customer_pca.py，图片写入当前目录。

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

rng = np.random.default_rng(2026)
n = 120
latent = rng.normal(size=n)
# 潜在活跃度影响订单强度，订单为整数；后两列仍有各自独立噪声。
orders = rng.poisson(np.exp(2.2 + 0.35 * latent))
amount = 100 + 45 * orders + 30 * latent + rng.normal(0, 30, n)
minutes = 40 + 8 * latent + 2 * orders + rng.normal(0, 4, n)
X = np.column_stack([orders, amount, minutes])
# axis=0 跨用户计算每列统计量；不能沿 axis=1 把每个用户的三种单位混合。
mean = X.mean(axis=0)
sd = X.std(axis=0, ddof=1)
assert np.all(sd > 0)
Z = (X - mean) / sd
assert np.allclose(Z.mean(axis=0), 0, atol=1e-12)
assert np.allclose(Z.std(axis=0, ddof=1), 1)
# C 是标准化样本的协方差，也就是原三特征的样本相关矩阵。
C = Z.T @ Z / (n - 1)
U, s, Vh = np.linalg.svd(Z, full_matrices=False)
directions = Vh.T
scores = Z @ directions
ratio = s**2 / np.sum(s**2)
# 用独立的对称特征值算法核对，eigh 默认升序，须同时重排特征向量。
eigenvalues, eigenvectors = np.linalg.eigh(C)
order = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[order]
eigenvectors = eigenvectors[:, order]
assert np.allclose(eigenvalues, s**2 / (n - 1))
assert np.allclose(ratio, eigenvalues / eigenvalues.sum())
# 向量可能整体异号；外积给出同一一维子空间的正交投影矩阵。
for j in range(3):
    assert np.allclose(np.outer(directions[:, j], directions[:, j]),
                       np.outer(eigenvectors[:, j], eigenvectors[:, j]))
assert np.allclose(scores, U * s)
assert np.allclose(scores.T @ scores / (n - 1), np.diag(eigenvalues))
# 保留前两成分重建标准化数据，再还原单位；误差先在无量纲 Z 上度量。
reconstructed_z = scores[:, :2] @ directions[:, :2].T
reconstructed_x = reconstructed_z * sd + mean
relative_error = np.linalg.norm(Z - reconstructed_z) / np.linalg.norm(Z)
assert np.isclose(relative_error, np.sqrt(ratio[2]))
print("三列样本均值（次、元、分钟）：", mean)
print("三列样本标准差（次、元、分钟）：", sd)
print("标准化样本协方差：\n", C)
print("主成分解释率：", ratio)
print("主成分方向（每列一个方向）：\n", directions)
print("前两成分的标准化重建相对误差：", relative_error)
# 主成分得分不相关不等于独立；这里不进行检验、置信区间或因果推断。
# 未标准化时金额的大数值会主导平方损失，这是另一种分析问题。

configure_chinese_plot()
fig, axes = plt.subplots(1, 3, figsize=(6.3, 2.9), constrained_layout=True)
axes[0].scatter(orders, amount, alpha=0.6, s=12)
axes[0].set(title="明确模拟的原始画像", xlabel="月订单数（次）",
            ylabel="月消费金额（元）")
axes[1].scatter(scores[:, 0], scores[:, 1], alpha=0.6, s=12)
axes[1].set(title="标准化后的二维得分", xlabel="第一主成分得分",
            ylabel="第二主成分得分")
axes[2].bar([1, 2, 3], ratio)
for index, value in enumerate(ratio, 1):
    axes[2].text(index, value + 0.025, f"{value:.1%}", ha="center", fontsize=7)
axes[2].set(title="标准化方差的解释率", xlabel="主成分序号",
            ylabel="解释率", ylim=(0, 1.1))
axes[2].set_xticks([1, 2, 3])
for ax in axes:
    ax.grid(alpha=0.2)
fig.savefig("customer_pca.png", dpi=180, bbox_inches="tight")
plt.close(fig)