第09回主成分分析（PCA）¶

科目：ビッグデータ応用Ⅱ〔特徴量エンジニアリング〕 担当：小野原彩香 対象：3年生（後期）

学習目標¶

この回を終えると、以下のことができるようになります。

「次元の呪い」が引き起こす問題を具体例で説明できる
主成分分析（PCA）の原理（固有値・固有ベクトル）を理解する
寄与率・累積寄与率をもとに主成分数を選択できる
scikit-learnを使ってPCAを実装し、スクリープロットとバイプロットを描ける
空間特徴量における標準化の重要性を説明できる
発展：第1主成分が「射影の分散が最大になる向き」であることと、寄与率が固有値の割合になることを、自分で導ける

1. 次元削減の必要性¶

1.1 次元の呪い（Curse of Dimensionality）¶

地理空間LLMのログデータから抽出した特徴量を想像してください。

密度系：クエリ密度、POI密度、道路密度
距離系：最寄り駅距離、市街地からの距離
面積系：緑地面積率、建物被覆率
時間系：ピーク時間帯、利用頻度の変動係数
テキスト系：エラーメッセージの長さ、地名のあいまい度スコア

これらを組み合わせると、30〜50次元以上になることがあります。

次元が増えると何が困るのか？

直感的な説明：

1次元（線）の上に100点のデータを置くと、点と点の平均距離はそれなりに密です。これを2次元（正方形）に拡張すると、同じ100点で同じ密度を保つには 10倍 のデータが必要です。 10次元になると、10億倍 のデータが必要になります。

つまり、次元が増えるほどデータが空間を埋めるのに必要なサンプル数が爆発的に増えます。手元のデータは実質的に「スカスカ」になり、機械学習モデルが不安定になります。

問題のまとめ：

問題	内容
過学習リスク	次元が増えると、偶然のパターンを学習しやすくなる
距離指標の劣化	高次元では「ユークリッド距離」が意味をなさなくなる
計算コスト増大	学習・推論の時間が指数的に増える
可視化の困難	3次元を超えると人間が直接見えない

1.2 解決策としての次元削減¶

次元削減には大きく2種類あります。

特徴選択（Feature Selection）：元の特徴量の中から重要なものを選ぶ
特徴抽出（Feature Extraction）：元の特徴量を組み合わせて新しい軸を作る

PCAは「特徴抽出」の代表手法です。データの分散が最大になる方向を新しい軸として定義します。

2. 主成分分析（PCA）の原理¶

2.1 直感的な理解¶

PCAは「データが一番ばらついている方向」を順番に見つけます。

たとえば、縦長に分布した2次元データがあるとします。 - 第1主成分：データが最もばらつく方向（長軸） - 第2主成分：第1主成分に直交し、残りのばらつきを最大化する方向（短軸）

これを多次元に拡張したものがPCAです。

2.2 数学的な手順¶

ステップ1：標準化

各特徴量を平均0、分散1にします。これをしないと、スケールが大きな特徴量がPCAを支配します。

\[z_i = \frac{x_i - \bar{x}}{s}\]

ステップ2：共分散行列の計算

標準化後のデータ行列 \(Z\)（\(n \times p\)）から共分散行列 \(\Sigma\) を計算します。

\[\Sigma = \frac{1}{n-1} Z^T Z\]

ステップ3：固有値分解

共分散行列を固有値分解します。

\[\Sigma v = \lambda v\]

\(\lambda\)（固有値）：その主成分が「どれだけのばらつきを説明するか」
\(v\)（固有ベクトル）：主成分の「向き」（元の特徴量の線形結合）

固有値が大きい順に固有ベクトルを並べたものが主成分です。

導出①：第1主成分は「射影したときの分散が最大になる向き」である

「分散が最大になる方向」という言葉を式にします。標準化済みデータ \(Z\) を、長さ1のベクトル \(w\) の向きに射影したスコアは

\[t = Zw\]

です。このスコアの分散は、共分散行列 \(\Sigma\) を使って

\[\mathrm{Var}[t] = w^{T}\Sigma w\]

と書けます。PCAがやっているのは「\(\|w\| = 1\) のもとで \(w^{T}\Sigma w\) を最大にする \(w\) を探す」ことです。その答えが、最大固有値に対応する固有ベクトルになります。

2変数なら向きを角度 \(\theta\) ひとつで書けるので、実際に計算して確かめられます。

数値で確かめる：5点のデータから共分散行列が

\[\Sigma = \begin{pmatrix} 2.5 & 1.5 \\ 1.5 & 1.0 \end{pmatrix}\]

になったとします。\(w = (\cos\theta,\ \sin\theta)\) として \(w^{T}\Sigma w\) を計算すると、

\(\theta\)	0°	15°	30°	45°	60°	90°	120°
分散	2.500000	3.149519	3.424038	3.250000	2.674038	1.000000	0.075962

30°付近が最大です。一方この \(\Sigma\) の固有値は \(\lambda_1 = 3.427051\)、\(\lambda_2 = 0.072949\) で、第1固有ベクトルの向きは \(\theta = 31.7175°\)。そこでの分散はちょうど \(\lambda_1 = 3.427051\) になります。「分散が最大の向き」と「第1固有ベクトル」が同じものだと確かめられました。

同時に、\(\theta = 120°\) 付近（第1主成分と直交する向き）で分散が最小の \(\lambda_2\) に近づいていることも読めます。

一般の \(p\) 次元でこれを証明するには制約つき最大化（ラグランジュの未定乗数法）が必要で、この授業で扱う範囲外です。2変数で「そうなっている」ことを確かめるところまでにします。

使った数学：ベクトルの内積と射影。三角比。行列とベクトルの積。

2.3 寄与率と累積寄与率¶

指標	定義	読み方
寄与率	\(\lambda_k / \sum \lambda\)	その主成分が全分散の何%を説明するか
累積寄与率	\(\sum_{i=1}^{k} \lambda_i / \sum \lambda\)	k個の主成分で全分散の何%を説明するか

一般的な選択基準： - 累積寄与率が 70〜80% を超える主成分数を選ぶ - スクリープロット（後述）の「肘」の位置を選ぶ

導出②：寄与率が「固有値の割合」になる理由

寄与率の分母が \(\sum \lambda\) になっているのはなぜか。鍵は、分散の総和が座標の取り方によらないという性質です。

標準化したので元の各変数の分散はすべて1です。変数が \(p\) 個なら分散の総和は \(p\)。これは共分散行列の対角成分の和、つまり \(\mathrm{trace}(\Sigma)\) です。一方、固有値分解すると

\[\mathrm{trace}(\Sigma) = \sum_{k=1}^{p} \lambda_k\]

が成り立ちます。主成分は元の軸を回転させたものですが、回転しても分散の総和は変わりません。総和が保存されるからこそ、「全体のうち何割か」という言い方ができます。

\[\text{第}k\text{主成分の寄与率} = \frac{\lambda_k}{\sum_{j} \lambda_j} = \frac{\text{その軸の分散}}{\text{分散の総和}}\]

数値で確かめる：導出①の \(\Sigma\) では \(\mathrm{trace}(\Sigma) = 2.5 + 1.0 = 3.500000\)、固有値の和も \(3.427051 + 0.072949 = 3.500000\) で一致します。第1主成分の寄与率は \(3.427051 / 3.500000 = 0.979157\)、つまり97.9%です。

使った数学：対角成分の和（トレース）。割合の定義。

3. 空間特徴量へのPCA適用¶

3.1 空間特徴量のスケール問題¶

地理空間データの特徴量は単位がバラバラです。

特徴量	典型的な値	単位
クエリ密度	0.01〜5.0	件/km²
最寄り駅距離	100〜5000	m
建物被覆率	0.0〜1.0	比率
POI数	0〜500	件

これをそのままPCAにかけると、「最寄り駅距離（数百〜数千）」がほぼすべての第1主成分を占領してしまいます。

必ずStandardScalerで標準化してからPCAを適用します。

導出③：標準化しないと、単位を変えただけで主成分が変わる

「スケールが大きな特徴量がPCAを支配する」を、数値で確かめます。

導出①の5点データで、2列目の単位だけを1000倍にします。kmをmに直したような操作で、データの中身は何も変わっていません。

すると第1主成分の向きは

\[w_1 = (0.0015,\ 0.999999)\]

つまりほとんど第2変数そのものになり、寄与率は \(0.99999975\)（99.99997%）に跳ね上がります。導出①では第1主成分が2変数を混ぜた向き（\(\theta = 31.7175°\)）で寄与率97.9%だったのに、単位を変えただけで別物になりました。

理由は導出②から分かります。PCAは分散の総和のうち大きい取り分を探す手続きなので、ある変数の分散だけが \(1000^2\) 倍になれば、その変数が総和のほぼ全部を占めてしまいます。PCAは「どの単位で測ったか」を知りません。大きい数字が並んでいる列を重要だと判断するだけです。

標準化すればすべての変数の分散が1に揃い、この問題は起きません。この資料の特徴量は「m」「件/km²」「比率」が混在しているので、標準化は省略できる手順ではありません。

使った数学：定数倍したときの分散の変化（\(\mathrm{Var}[aX] = a^2\mathrm{Var}[X]\)）。導出②の寄与率の定義。

3.2 標準化の方法¶

from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)  # fit()とtransform()を同時実行

注意：fit()はあくまで訓練データで行い、テストデータにはtransform()だけ適用します。

4. 主成分の解釈¶

4.1 因子負荷量（ローディング）¶

各主成分が元の変数とどれだけ関係しているかを示す値です。

\[\text{ローディング} = v_k \times \sqrt{\lambda_k}\]

値が正の場合はその変数と正の相関、負の場合は負の相関があることを意味します。

4.2 スコア（主成分得点）¶

各データポイントが新しい主成分空間でどこに位置するかを示します。

\[T = Z \cdot V\]

\(Z\)：標準化済みデータ行列
\(V\)：固有ベクトル行列
\(T\)：主成分スコア行列

4.3 バイプロットの読み方¶

バイプロットは「スコア」と「ローディング」を1つのグラフに重ね表示します。

点（スコア）：各観測データの位置
矢印（ローディング）：元の変数の方向と強さ

読み方： - 矢印が長い → その変数は主成分の説明に貢献している - 矢印の向きが近い → それらの変数は正の相関がある - 矢印の向きが反対 → 負の相関がある

5. ハンズオン手順¶

5.1 環境準備¶

# 必要なライブラリのインストール
pip install scikit-learn matplotlib seaborn pandas numpy

# データの確認
python -c "import pandas as pd; df = pd.read_parquet('data/spatial_features.parquet'); print(df.shape)"

5.2 データ読み込みと前処理¶

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
import seaborn as sns
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA

# 日本語フォント設定（macOS）
plt.rcParams['font.family'] = 'Hiragino Sans'
plt.rcParams['axes.unicode_minus'] = False

# データ読み込み
df = pd.read_parquet('data/spatial_features.parquet')
print("データ形状:", df.shape)
print("\n特徴量一覧:")
print(df.dtypes)

# 数値列のみ選択
feature_cols = [
    'query_density',      # クエリ密度（件/km²）
    'poi_density',        # POI密度
    'road_density',       # 道路密度
    'nearest_station_dist',  # 最寄り駅距離（m）
    'green_area_ratio',   # 緑地面積率
    'building_coverage',  # 建物被覆率
    'elevation_mean',     # 平均標高
    'population_density', # 人口密度
    'error_rate',         # エラー発生率（後でターゲットとして使用）
]

X = df[feature_cols].dropna()
feature_names = feature_cols[:-1]  # error_rateを除く
X_feat = X[feature_names].values

print(f"\n使用する特徴量数: {len(feature_names)}")
print(f"サンプル数: {len(X_feat)}")

5.3 標準化¶

# 標準化
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X_feat)

# 標準化前後の比較
print("=== 標準化前 ===")
print(pd.DataFrame(X_feat, columns=feature_names).describe().round(2))

print("\n=== 標準化後 ===")
print(pd.DataFrame(X_scaled, columns=feature_names).describe().round(2))
# 平均≈0、標準偏差≈1になっていることを確認

5.4 PCAの実行¶

# PCA実行（全主成分）
pca = PCA()
pca.fit(X_scaled)

# 寄与率の確認
explained_variance_ratio = pca.explained_variance_ratio_
cumulative_variance_ratio = np.cumsum(explained_variance_ratio)

print("=== 主成分ごとの寄与率 ===")
for i, (ev, cv) in enumerate(zip(explained_variance_ratio, cumulative_variance_ratio)):
    print(f"PC{i+1}: 寄与率={ev:.3f} ({ev*100:.1f}%), 累積寄与率={cv:.3f} ({cv*100:.1f}%)")

5.5 スクリープロットの描画¶

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# スクリープロット（固有値）
eigenvalues = pca.explained_variance_
axes[0].plot(range(1, len(eigenvalues)+1), eigenvalues, 'bo-', markersize=8)
axes[0].axhline(y=1, color='r', linestyle='--', label='Kaiser基準（固有値=1）')
axes[0].set_xlabel('主成分番号')
axes[0].set_ylabel('固有値')
axes[0].set_title('スクリープロット')
axes[0].legend()
axes[0].xaxis.set_major_locator(ticker.MaxNLocator(integer=True))

# 累積寄与率プロット
axes[1].bar(range(1, len(explained_variance_ratio)+1),
            explained_variance_ratio * 100,
            alpha=0.6, label='寄与率')
axes[1].plot(range(1, len(cumulative_variance_ratio)+1),
             cumulative_variance_ratio * 100,
             'ro-', markersize=8, label='累積寄与率')
axes[1].axhline(y=80, color='g', linestyle='--', label='80%基準線')
axes[1].set_xlabel('主成分番号')
axes[1].set_ylabel('寄与率（%）')
axes[1].set_title('寄与率・累積寄与率')
axes[1].legend()
axes[1].xaxis.set_major_locator(ticker.MaxNLocator(integer=True))

plt.tight_layout()
plt.savefig('output/scree_plot.png', dpi=150, bbox_inches='tight')
plt.show()

# 累積寄与率80%に必要な主成分数
n_components_80 = np.argmax(cumulative_variance_ratio >= 0.80) + 1
print(f"\n累積寄与率80%に必要な主成分数: {n_components_80}")

5.6 主成分空間への変換¶

# 選択した主成分数でPCAを再実行
n_components = n_components_80  # または手動で設定

pca_final = PCA(n_components=n_components)
X_pca = pca_final.fit_transform(X_scaled)

print(f"変換後の形状: {X_pca.shape}")
print(f"（{X_feat.shape[1]}次元 → {n_components}次元に削減）")

# 主成分スコアをDataFrameに変換
pc_cols = [f'PC{i+1}' for i in range(n_components)]
df_pca = pd.DataFrame(X_pca, columns=pc_cols)
df_pca['error_rate'] = X['error_rate'].values

5.7 バイプロットの描画¶

def biplot(pca, feature_names, X_pca, figsize=(10, 8)):
    """
    バイプロットを描画する関数
    - 点: 各観測値のPC1-PC2スコア
    - 矢印: 元の変数のローディング
    """
    fig, ax = plt.subplots(figsize=figsize)

    # スコアの散布図（最初の200点のみ表示）
    n_plot = min(200, len(X_pca))
    ax.scatter(X_pca[:n_plot, 0], X_pca[:n_plot, 1],
               alpha=0.4, s=20, c='steelblue', label='観測値')

    # ローディングベクトルの矢印
    loadings = pca.components_.T  # shape: (n_features, n_components)
    scale = 3.0  # 矢印のスケール調整

    for i, (name, load) in enumerate(zip(feature_names, loadings)):
        ax.annotate('', xy=(load[0]*scale, load[1]*scale), xytext=(0, 0),
                    arrowprops=dict(arrowstyle='->', color='crimson', lw=2))
        ax.text(load[0]*scale*1.1, load[1]*scale*1.1, name,
                fontsize=9, color='darkred', ha='center')

    ax.axhline(y=0, color='gray', linestyle='-', alpha=0.3)
    ax.axvline(x=0, color='gray', linestyle='-', alpha=0.3)
    ax.set_xlabel(f'PC1 (寄与率: {pca.explained_variance_ratio_[0]*100:.1f}%)')
    ax.set_ylabel(f'PC2 (寄与率: {pca.explained_variance_ratio_[1]*100:.1f}%)')
    ax.set_title('バイプロット（PC1 vs PC2）')
    plt.tight_layout()
    return fig

fig = biplot(pca_final, feature_names, X_pca)
plt.savefig('output/biplot.png', dpi=150, bbox_inches='tight')
plt.show()

5.8 因子負荷量ヒートマップ¶

# 因子負荷量の計算
loadings_df = pd.DataFrame(
    pca_final.components_.T,
    index=feature_names,
    columns=pc_cols
)

# ヒートマップ
plt.figure(figsize=(10, 6))
sns.heatmap(loadings_df, annot=True, fmt='.2f', cmap='RdBu_r',
            center=0, vmin=-1, vmax=1,
            linewidths=0.5, cbar_kws={'label': '因子負荷量'})
plt.title('因子負荷量ヒートマップ')
plt.tight_layout()
plt.savefig('output/loadings_heatmap.png', dpi=150, bbox_inches='tight')
plt.show()

print("\n=== 各主成分の主要な寄与変数 ===")
for col in pc_cols:
    top_pos = loadings_df[col].nlargest(3)
    top_neg = loadings_df[col].nsmallest(2)
    print(f"\n{col}:")
    print(f"  正の方向: {', '.join([f'{v}({s:.2f})' for v, s in top_pos.items()])}")
    print(f"  負の方向: {', '.join([f'{v}({s:.2f})' for v, s in top_neg.items()])}")

6. 演習課題¶

課題1（基本）¶

data/spatial_features.parquetを読み込み、以下を実行してください。

全特徴量のスケールの違いをboxplotで可視化する
StandardScalerで標準化後、PCAを実行する
スクリープロットを描画し、累積寄与率80%に必要な主成分数を求める

課題2（応用）¶

PCAで3次元に削減したデータを3D散布図で可視化する（matplotlibのAxes3Dを使用）
バイプロットを描画し、PC1とPC2が何を表しているか、因子負荷量から解釈してください（100字以内のコメントをコードに記述すること）

課題3（発展）¶

PCA適用前と適用後で、ランダムフォレストによるエラー発生率予測精度（RMSE）を比較する
主成分数を2, 3, 4, 5と変えたときの予測精度の変化をグラフ化する

提出形式：Jupyterノートブック（.ipynb）、ファイル名：09_pca_[学籍番号].ipynb

まとめ¶

項目	ポイント
次元の呪い	次元増加で距離指標が劣化、データが「スカスカ」になる
PCAの基本	分散最大方向を主成分として抽出（固有値分解）
標準化の重要性	空間特徴量は単位が異なるため、必ずStandardScalerを適用
主成分数の選択	累積寄与率80%基準またはスクリープロットの肘を使う
バイプロット	スコアとローディングを重ね、特徴量の関係を可視化
次回への接続	削減した主成分スコアを第10回の時空間キューブ分析に活用する

次回予告： 第10回では、LLMログデータの「時間次元」を加えた時空間データキューブを構築します。

作成：小野原彩香 / ビッグデータ応用Ⅱ〔特徴量エンジニアリング〕

第09回 主成分分析（PCA）¶