主成分分析のシミュレーション

pca.py

#!/usr/bin/env python3
# coding: utf-8
# プログラミングのための確率統計 p.280 例題8.2 の numpy での実装
import numpy as np


# データ群
x1 = np.array([0, 5])
x2 = np.array([0, -5])
x3 = np.array([4, 3])
x4 = np.array([-4, -3])
X = np.array([x1, x2, x3, x4]).T

# 共分散行列を求める（どれで計算しても最後の固有値には影響しない）
# X_cov = np.cov(X)  # numpy
# 期待値を求める
# X_mean = np.mean(X)
# X_cov = np.dot(X - X_mean, (X - X_mean).T)  # まとめて
X_cov = np.dot(x1.reshape(2, 1), x1.reshape(1, 2)) \
    + np.dot(x2.reshape(2, 1), x2.reshape(1, 2)) \
    + np.dot(x3.reshape(2, 1), x3.reshape(1, 2)) \
    + np.dot(x4.reshape(2, 1), x4.reshape(1, 2))  # 1つ1つ

# 固有値計算
X_vals, X_vecs = np.linalg.eig(X_cov)

# 大きいほうの固有値に対応する固有ベクトルを取得
max_vec = X_vecs[np.argmax(X_vals)]
# 大きさを1に正規化（もともとなってはいる）
q = max_vec / np.linalg.norm(max_vec)

# 各データの第1主成分計算
x1_pca = np.dot(q, x1)
x2_pca = np.dot(q, x2)
x3_pca = np.dot(q, x3)
x4_pca = np.dot(q, x4)
print(x1_pca, x2_pca, x3_pca, x4_pca)