7.4 應用：PCA 降維#

import numpy as np

rng = np.random.default_rng(0)
X = rng.multivariate_normal([0, 0], [[3, 2.4], [2.4, 3]], size=300)
print(X.shape)          # (300, 2)：兩個高度相關的特徵
print(np.corrcoef(X.T)) # 相關係數約 0.8——一個座標幾乎夠用

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(0)
X = rng.multivariate_normal([0, 0], [[3, 2.4], [2.4, 3]], size=300)
Xc = X - X.mean(axis=0)                       # 中心化
w, V = np.linalg.eigh(Xc.T @ Xc / (len(X) - 1))   # 對稱矩陣用 eigh

plt.figure(figsize=(5, 5))
plt.scatter(Xc[:, 0], Xc[:, 1], s=10, alpha=0.5)
for i, lam in enumerate(w[::-1]):             # 由大到小畫主成分
    v = V[:, ::-1][:, i] * 2.5 * np.sqrt(lam)
    plt.arrow(0, 0, v[0], v[1], width=0.03, color=f"C{i}",
              label=f"$PC_{i+1}$ ($\\lambda$={lam:.2f})")
plt.axis("equal"); plt.legend(); plt.title("PCA：找到變異數最大的方向")
plt.show()

import numpy as np

rng = np.random.default_rng(1)
X = rng.random((200, 5)) * np.array([10, 5, 1, 0.5, 0.1])
Xc = X - X.mean(axis=0)

# 路線一：特徵分解協方差矩陣
w, V = np.linalg.eigh(Xc.T @ Xc / (len(X) - 1))
w, V = w[::-1], V[:, ::-1]                    # eigh 回傳由小到大，翻轉

# 路線二：SVD
U, s, Vt = np.linalg.svd(Xc, full_matrices=False)
w2 = s**2 / (len(X) - 1)

print(np.allclose(w, w2))                     # True：兩條路殊途同歸
print(np.allclose(np.abs(V), np.abs(Vt.T)))   # 方向可能差一個正負號

import numpy as np
import matplotlib.pyplot as plt

def pca(X, n_components=2, standardize=True):
    """標準 PCA 流程：標準化 → 中心化 → SVD → 投影"""
    if standardize:
        X = (X - X.mean(axis=0)) / X.std(axis=0)   # 標準化：每個特徵均值 0、標準差 1
    Xc = X - X.mean(axis=0)                        # 中心化
    U, s, Vt = np.linalg.svd(Xc, full_matrices=False)
    Z = Xc @ Vt[:n_components].T                   # 投影到前 n 個主成分
    ratio = s**2 / (s**2).sum()                    # 變異數解釋比例
    return Z, Vt[:n_components], ratio

# 延續上一格的 X（5 個特徵）
Z, components, ratio = pca(X)
print(Z.shape)          # (200, 2)：5 維壓到 2 維
print(ratio[:2])        # 前兩個主成分解釋了多少變異

plt.scatter(Z[:, 0], Z[:, 1], s=15, alpha=0.7)
plt.xlabel("$PC_1$"); plt.ylabel("$PC_2$")
plt.title("投影到前兩個主成分")
plt.show()

import numpy as np
import matplotlib.pyplot as plt

def scree_plot(X, standardize=True):
    if standardize:
        X = (X - X.mean(axis=0)) / X.std(axis=0)
    Xc = X - X.mean(axis=0)
    s = np.linalg.svd(Xc, compute_uv=False)
    ratio = s**2 / (s**2).sum()

    fig, ax = plt.subplots(1, 2, figsize=(10, 4))
    ax[0].plot(np.arange(1, len(s)+1), ratio, "o-")
    ax[0].set_xlabel("主成分 $i$"); ax[0].set_ylabel("解釋比例")
    ax[0].set_title("Scree plot")
    ax[1].plot(np.arange(1, len(s)+1), np.cumsum(ratio), "o-", color="tomato")
    ax[1].axhline(0.95, ls="--", c="gray")
    ax[1].set_xlabel("主成分數 $k$"); ax[1].set_ylabel("累積解釋比例")
    ax[1].set_title("累積解釋比例（虛線：95%）")
    plt.tight_layout(); plt.show()

scree_plot(X)

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris

iris = load_iris()
X, y = iris.data, iris.target

Z, components, ratio = pca(X)          # 沿用前面的 pca 函數
print("前兩個主成分解釋比例：", ratio[:2].round(3))     # 約 [0.73, 0.23]

plt.figure(figsize=(7, 5))
for label, name, color in zip([0, 1, 2], iris.target_names, ["C0", "C1", "C2"]):
    plt.scatter(Z[y == label, 0], Z[y == label, 1],
                label=name, s=25, alpha=0.8, color=color)
plt.xlabel(f"$PC_1$（解釋 {ratio[0]:.0%}）")
plt.ylabel(f"$PC_2$（解釋 {ratio[1]:.0%}）")
plt.legend(); plt.title("鳶尾花：4 維資料的 2 維 PCA 視覺化")
plt.show()

print("各主成分的組成（欄 = 原特徵權重）：")
print(np.round(components, 2))         # 看第一主成分主要由哪些特徵構成

scree_plot(X)      # 鳶尾花的累積曲線：前 2 個主成分就超過 95%