「入門はじめての多変量解析」をPythonで写経 Vol.9 ～ ３章「はじめての主成分分析」①主成分分析をはじめる前に

３章「はじめての主成分分析」

書籍の著者　石村貞夫 先生、石村光資郎 先生

書籍「入門はじめての多変量解析」３章「はじめての主成分分析」の Python写経活動記録 です。 

多変量解析の入門を Python と一緒に学ぶ写経シリーズです。

この記事は、書籍の最初のデータで ざっくり主成分分析 に取り組みます。
書籍の計算とPython ライブラリの結果を比べて、次回以降の主成分分析深堀りに備えます。

とても長い旅になりそうです。
ひとまず 「データの偏差平方和積和行列」の主成分分析を最初のマイルストーンにして 頑張りましょう。

ChatGPT 活用型学習で進めてまいります！
では書籍を開いて多変量解析の旅に出かけましょう🚀

はじめに

このブログシリーズは、書籍「入門はじめての多変量解析」（東京図書、「テキスト」と呼びます）の Python 写経を通じて得た「多変量解析の楽しさ」をご紹介します。

書籍の紹介と引用表記はリンク先の記事に掲載しています。

３章 はじめての主成分分析

この記事は３章の以下のSectionを取り扱います。

記事に用いるデータは、テキストに掲載されたデータそのものを引用しています。
データ件数の少ないものはコード上でデータを登録し、データ件数の多いものはCSVファイル化してデータを読み込みしています。

この記事で用いるライブラリをインポートします。

### インポート

# 数値計算
import string        # 大文字アルファベットを取得
import numpy as np
import pandas as pd

# PCA
from sklearn.decomposition import PCA
from statsmodels.multivariate.pca import PCA as sm_PCA  # sklearnと名前重複を回避

# 描画
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo'  # または import japanize_matplotlib

イントロダクション

この記事から「主成分分析」に取り組みます！
「主成分分析って何？」を ChatGPTにざっくり訊いてみました。

1️⃣ 主成分分析って何？どんなことができる？

• 本質イメージ
  複数の変数を１本—or 少数—の「要約指標」（総合指標とも）にまとめ、データの大まかな傾向をつかむ手法です。

• できること

  • 変数の次元を削減し、扱いやすい指標でデータを要約

  • データの全体像を可視化して俯瞰

  • 雑音を取り除き、特徴を抽出

◆

2️⃣ 主成分分析の利用シーン（例）

• アンケートの要約：
  20〜30項目の設問を、３〜５本の合成指標にまとめて傾向把握

• センサーデータの総合化：
  複数センサーの信号を１本の「状態スコア」に圧縮

• マーケティング指標づくり：
  売上・アクセス・滞在時間…複数数値を集約した「興味度スコア」の作成

◆

3️⃣ 数学的にはどんなことをしている？

1. 相関・共分散行列 を計算

2. 固有値・固有ベクトル を求める

   • 固有ベクトル＝要約軸（データのばらつきが大きい方向）

   • 固有値＝その軸の“情報量”（分散の大きさ）

3. 上位の固有ベクトルを「主成分」と呼び、
   元データをその軸に投影して得る値が 主成分得点（要約スコア）

4. 寄与率／累積寄与率 で「要約指標として何本取れば十分か」を判断

◆

4️⃣ 学習を進める上でぜひ押さえたいポイント

• 標準化の重要性：
  異なる単位の変数を公平に要約するため

• スクリープロット（寄与率グラフ）：
  情報を十分残しつつ次元削減する本数の目安

• 主成分負荷量（loadings）の解釈：
  各元変数が“要約軸”にどれだけ寄与しているか

• 要約の落とし穴：

  • 次元削減で失われる情報もある

  • 抽出した要約軸は「元の変数名そのまま」ではないので、解釈は「総合的な意味合い」で行う

ChatGPT さん、ありがとうございます！
想像できたでしょうか？

改めますと…

💡「要約」「総合化」「次元削減」「教師なし学習」がポイント

• 主成分分析はデータを「要約」して「少数の新しい変数」をつくります。

  • 新しい変数はデータの特徴を最大限残して「総合化」するものです。

  • 「次元削減」の手法とも呼ばれ、新しい変数間の相関係数は０です。

• 重回帰分析との違いは目的変数がないこと。
  機械学習の文脈では「教師なし学習」に該当します。

◆

📖 今回記事でやることの概要

① 主成分分析の体感
テキストの 表 3.1.1 のデータを用いて、主成分分析を概観します。
主成分分析のイメージをざっくり感じてみましょう。

② 計算の体感
テキストの計算方法とPythonライブラリの計算結果の違いをざっくり眺めます。

◆

📖 今回記事でやらないことの概要

• テキストの「情報損失量の最小化」には取り組みません。
  代わりに「分散の最大化」に着目して進めます。
  なおテキストでは、分散の最大化は Section 3.9 以降で実践します。

• 「データの標準化」を行いません。
  つまり「相関行列」ベースの主成分分析には触れていません。
  テキストの前半部がそうなっているからです。

◆

⚠️ 注意書き：用語の定義 ⚠️
主成分分析のさまざまな書籍・Webサイトの情報を横断的に眺めると、用語の意味が異なる場面にによく遭遇します。
例えば、上記ChatGPTの主成分（固有ベクトルと同義）は、この記事の主成分とは異なっています。
特に主成分負荷量（Loadings）の多様さには難儀します。
記事の用語遣いがみなさんの理解と異なる場合には、すみませんm(_ _)m

では分析データの確認からはじめていきます！

分析データ

テキスト p.87 表 3.1.1「介護施設と医療施設（その１）」のデータをお借りします。

### 介護・医療施設データ p.87 表3.1.1

# データの登録
N = 10
data1 = pd.DataFrame(
    {'地域名': list(string.ascii_uppercase[:N]),
     '介護施設': [22, 22, 18, 18, 15, 19, 19, 24, 21, 25],
     '医療施設': [12, 8, 6, 15, 7, 9, 7, 17, 14, 11]},
     index=range(1, N+1))
data1.index.name = 'No.'
data1

【実行結果】
10 の地域における 65 歳以上人口 1 万人あたりの介護施設数、人口 1 万人あたりの医療施設数です。
２つの施設数は単位が異なるのでご注意ください。

要約統計量と相関係数を見ましょう。

# 要約統計量
data1.describe().round(1).T

【実行結果】
特段ありません。。。

# 相関係数
data1.corr(numeric_only=True).round(3)

【実行結果】
２つの変数の相関係数 $${0.5}$$ は中程度の正の相関です。

各変数のヒストグラムと２変数の散布図を見ましょう。

# ヒストグラムの描画
fig, axes = plt.subplots(1, 2, figsize=(8, 4), sharey=True, tight_layout=True)
for col, ax in zip(data1.columns[1:], axes.flat):
    sns.histplot(data=data1[col], kde=True, ec='white', ax=ax)

【実行結果】
介護施設は山型になっています。
医療施設は左側（小さい値）に偏り気味です。

### 散布図の描画 p.88 図3.2.1

# 描画領域の設定
fig, ax = plt.subplots(figsize=(6, 6))
# 散布図の描画
sns.scatterplot(data=data1, x='介護施設', y='医療施設', s=70, alpha=0.7, ax=ax)
# 施設名のテキストの表示
for (s, x, y) in data1.values:
    ax.text(x=x, y=y+0.2, s=s, fontsize=12)
# 修飾
ax.set(xlim=(13, 27), ylim=(4, 19), aspect='equal')
ax.grid(lw=0.5, alpha=0.5);

【実行結果】
散布図の各点には地域名を添えています。
散布図から読み取れる介護・医療分野への取り組みは：
・地域 H が力を入れている感じ
・地域 C, E あたりが力を入れていない感じ

単回帰分析では散布図に回帰直線を引きました。
主成分分析も散布図に線を引きます！
ただし、回帰直線と異なる線になる予定です。

■ これから確認していくものたちの紹介
本記事の献立てをご紹介します。
ツールで計算する方法（計算過程はブラックボックス）と、計算過程を明示する方法に大別できます。
計算に使う２種類の行列との組合せを踏まえて、次の順番で進めます。

1. 偏差平方和積和行列から計算
   計算過程を明示します。
   偏差平方和積和行列ベースの計算です。
   テキストが用いるツール SPSS と同等と想定しています（ChatGPT談）。
   テキストに近いという理由で一番最初に見ておきます。

2. scikit-learn で主成分分析
   機械学習の文脈でよく利用される Python のライブラリです。
   分散共分散行列ベースの計算です。
   偏差平方和積和行列ベースの計算との違いも確認します。

3. statsmodels で主成分分析
   統計解析の文脈でよく利用される Python のライブラリです。
   偏差平方和積和行列ベースの計算です。
   １の結果と若干違うことも確認します。

4. 分散共分散行列から計算
   計算過程を明示します。
   分散共分散行列ベースの計算です。

ベースにする行列とツールの関係です。

「データの偏差平方和積和行列」で主成分分析

ChatGPT によると、統計解析ソフト SPSS はデフォルトでデータの偏差平方和積和行列を用いるそうです。
テキストで使われるツールですので、まずは偏差平方和積和行列のケースから見ていきます。
特に注目したいことは：

データの偏差平方和積和行列の「固有値・固有ベクトル」

固有値・固有ベクトルは本シリーズ Vol.1 で肩慣らししたアレです！
こちらの記事でご確認ください。

■ 準備
主成分分析に利用する変数 X と、データ作成時に用いる列名等を準備します。

### データ等の準備

# 変数Xの設定 ※データから地域列を除いたもの
X = data1.iloc[:, 1:]

# 主成分の列名の設定　※PC1：第１主成分、PC2：第２主成分
columns = ['PC1', 'PC2']

# 地域名の取得
area = data1['地域名']

【実行結果】なし

■ 偏差平方和積和行列の算出
ひとまず Python で計算してみましょう。

## データのSSCP行列の作成 p.100 平方和積和行列

# 中心化した行列：偏差行列
Xc = X - X.mean(axis=0)

# SSCP行列の算出
S_sscp = Xc.T @ Xc

# 結果表示
print('SSCP行列:')
display(S_sscp)

【実行結果】
とてもシンプルなデータに変わりました！

【計算過程の確認】
偏差、平方和、積和に分解して個別に計算します。

1️⃣ 偏差
データからそのデータの平均を差し引いた結果を偏差と呼んでいます。
「データの中心化」とも呼ばれます。

Python で偏差を算出しましょう。

# 偏差の算出
deviation = X - X.mean(axis=0)
deviation

【実行結果】
２つの変数の偏差です。
平均値からの距離に相当します。バラツキのことです！

偏差を足すと０になります。
これではバラツキ度合いがわからないのです。

# 偏差の合計は０
deviation.sum(axis=0).rename('合計').to_frame().T.round(10)

【実行結果】

2️⃣ 平方和
個々の偏差を二乗（＝平方）して合計（＝和）します。
二乗しないで足すと０になりますが、二乗するとすべて正の値になるので、０になりません。
これでバラツキ度合いが分かります！

# 偏差平方和の算出

# 個々の偏差の二乗＝偏差平方を算出
deviation_squared = deviation**2
display(deviation_squared)

# 偏差平方を合計して偏差平方和
deviation_squared.sum().rename('偏差平方和').to_frame().T

【実行結果】
上の表が「偏差平方」です。
下の表が「偏差平方和」です。

下の表に注目しましょう。
「偏差平方和積和行列」の対角成分（左上から右下にかけての対角セル）になっています！

3️⃣ 積和
個々の偏差どうしを行単位で掛けて（＝積）して合計（＝和）します。
変数間の同調度合いです！

# 偏差積和の算出

# 個々の偏差どうしを掛ける＝偏差積を算出
deviation_product = deviation.prod(axis=1).rename('偏差積').to_frame()
display(deviation_product)

# 偏差積を合計して偏差積和
deviation_product.sum().rename('偏差積和').to_frame().T

【実行結果】上の表が「偏差積」です。
下の表が「偏差積和」です。

下の表に注目しましょう。
「偏差平方和積和行列」の対角以外の成分になっています！

偏差平方和積和行列は、テキストでは「平方和積和行列」と呼ばれ、英語呼称「Sum of Squares and Cross Products Matrix」の頭文字をとって「SSCP行列」とも呼ばれます。

これからは SSCP行列と呼んでいきます！

■ SSCP 行列の固有値・固有ベクトルの算出
SSCP行列の固有値・固有ベクトルは本日のメインディッシュ級です！
ギュッと圧縮した結果のSSCP行列から「本質」（みたいなもの）を取り出します！

numpy.linalg の eig() で固有値・固有ベクトルを簡単に算出できます！
固有値の大きな順に並び替えしています。

## SSCP行列の固有値・固有ベクトルの算出 p.100

# 固有値・固有ベクトルの算出
eig_val1_sscp, eig_vec1_sscp = np.linalg.eig(S_sscp)

# 固有値の降順ソートでindex取得
sort_idx = eig_val1_sscp.argsort()[::-1]
eig_val1_sscp = eig_val1_sscp[sort_idx]
eig_vec1_sscp = eig_vec1_sscp[:, sort_idx]

【実行結果】なし

固有値と固有ベクトルを表示しましょう。
pandas データフレームにして見やすくします。
PC1、PC2 は第１主成分、第２主成分の略称です。

# 固有値の表示
pd.DataFrame(eig_val1_sscp, index=columns, columns=['固有値']).T

【実行結果】
こちらは固有値です。

# 固有ベクトルの表示 ※テキストと比べて正負が逆転
pd.DataFrame(eig_vec1_sscp, columns=columns, index=X.columns)

【実行結果】
こちらは固有ベクトルです。
PC1 の固有値と固有ベクトル、PC2 の固有値と固有ベクトル、のように、固有値と固有ベクトルには対応関係があります。

固有値は主成分の重要度、固有ベクトルは主成分の方向を示すそうです。

■ 主成分とは！？
ここで「主成分」の正体を「式」で明かします！

$${z_1}$$ は第１主成分（略称 PC1）、$${z_2}$$ は第２主成分（略称 PC2）です。
「固有値の大きい順」で第１、第２の順序をとります。
主成分の数は最大で変数の数です。

$${z_1}$$ は第１主成分です。
固有値の大きい方の固有ベクトルを係数にしています。
続いて、$${z_2}$$ は第２主成分であり、固有値の小さい方の固有ベクトルを係数にしています。

$${z_1, z_2}$$ をデータの散布図上で示してみましょう。

## 主成分z1, z2の描画 p.110 図 3.8.1 の改造

## z1, z2の傾き・切片の算出 ※横軸x1, 縦軸x2のグラフ描画用
# x1, x2の平均値の算出
x1_mean, x2_mean = X.mean(axis=0).values
# z1,z2の傾きの算出
slope_z1 = eig_vec1_sscp[1, 0] / eig_vec1_sscp[0, 0]
slope_z2 = eig_vec1_sscp[1, 1] / eig_vec1_sscp[0, 1]
# z1,z2の切片の算出 ※z1,z2はx1,x2の平均を通る
intercept_z1 = x2_mean - slope_z1 * x1_mean
intercept_z2 = x2_mean - slope_z2 * x1_mean

## z1, z2の直線の算出
# z1のx軸:x1,y軸:x2の値
x1_z1_lines = np.linspace(16.5, 25, 2)
x2_z1_lines = slope_z1 * x1_z1_lines + intercept_z1
# z2のx軸:x1,y軸:x2の値
x1_z2_lines = np.linspace(16, 24.3, 2)
x2_z2_lines = slope_z2 * x1_z2_lines + intercept_z2

## 描画
# 描画領域の設定
fig, ax = plt.subplots(figsize=(6, 6))

# z1、z2の直線の描画
ax.plot(x1_z1_lines, x2_z1_lines, color='tab:green', label='第１主成分 $z_1$')
ax.plot(x1_z2_lines, x2_z2_lines, color='tab:red', label='第２主成分 $z_2$')

# 散布図の描画
sns.scatterplot(data=data1, x='介護施設', y='医療施設', s=70, alpha=0.7, ax=ax)
# 施設名のテキストの表示
for (s, x, y) in data1.values:
    ax.text(x=x, y=y+0.2, s=s, fontsize=12)

# x1,x2の平均点の描画
ax.plot([x1_mean], [x2_mean], '*', color='tab:orange', ms=15)
ax.text(x=x1_mean+0.5, y=x2_mean, va='center', s=(x1_mean, x2_mean))

# 修飾
ax.set(xlim=(13, 27), ylim=(4, 19), aspect='equal')
ax.grid(lw=0.5, alpha=0.5);

【実行結果】
十字形の直線が入りました。
緑色の直線が第１主成分 $${z_1}$$、赤色の直線が第２主成分 $${z_2}$$ です。
それぞれ右上向き・右下向きの「方向」がありますね！
２つの直線は垂直に交わっています。「直交」と呼びます。

【データ点と主成分の直線の関係 ～ 分散最大化～】
データ点を第１主成分の直線に「垂直に」移動させてみましょう。

## 主成分z1とデータ点の直交射影の描画

## z1の傾き・切片の算出 ※横軸x1, 縦軸x2のグラフ描画用
# x1, x2の平均値の算出
X_mean = X.mean(axis=0).values
# 傾きの算出
slope_z1 = eig_vec1_sscp[1, 0] / eig_vec1_sscp[0, 0]
# 切片の算出 ※z1,z2はx1,x2の平均を通る
intercept_z1 = X_mean[1] - slope_z1 * X_mean[0]

## z1の直線の算出
# z1のx軸:x1,y軸:x2の値
x1_z1_lines = np.linspace(16.5, 25, 2)
x2_z1_lines = slope_z1 * x1_z1_lines + intercept_z1

## データ点をz1へ直交射影
# 平坦な配列を列ベクトルに変換する関数の定義
vec = lambda x: x.reshape(-1, 1)
# 第１主成分への射影行列の算出
P1 = vec(eig_vec1_sscp[:, 0]) @ vec(eig_vec1_sscp[:, 0]).T
# データXを第１主成分z1へ直交射影
proj1 = vec(X_mean) + P1 @ (X.T - vec(X_mean))


## 描画
# 描画領域の設定
fig, ax = plt.subplots(figsize=(6, 6))

# z1の直線の描画
ax.plot(x1_z1_lines, x2_z1_lines, color='tab:green', label='第１主成分 $z_1$')

# 散布図の描画
sns.scatterplot(data=data1, x='介護施設', y='医療施設', s=70, alpha=0.7, ax=ax)
# 施設名のテキストの表示
for (s, x, y) in data1.values:
    ax.text(x=x, y=y+0.2, s=s, fontsize=12)

# 直交射影点の描画
ax.plot(*proj1.values, 'o', color='green', ms=3.5)

# データXから直交射影点への直線の描画
for i in range(N):
    ax.plot([X.iloc[i, 0], proj1.iloc[0, i]], [X.iloc[i, 1], proj1.iloc[1, i]],
            color='gray', lw=1, ls='--', zorder=0)

# x1,x2の平均点の描画
ax.plot(X_mean[0], X_mean[1], '*', color='tab:orange', ms=15)
ax.text(x=X_mean[0]+0.5, y=X_mean[1], va='center', s=(X_mean[0], X_mean[1]))

# 修飾
ax.set(xlim=(13, 27), ylim=(4, 19), aspect='equal')
ax.grid(lw=0.5, alpha=0.5);

【実行結果】
図の点線に沿って第１主成分の直線に「垂直に」移動するイメージです。
緑の点が変換後のデータ点です。

変換は「２次元のデータが１次元に変換される」ことを示しています。
第１主成分の１次元の軸で表現してみましょう。

## 主成分z1・１次元に要約された様子を描画

## 描画
# 描画領域の設定
fig, ax = plt.subplots(figsize=(6, 1))

# z1の直線の描画
zero_vec1 = np.zeros_like(x1_z1_lines)
ax.plot(x1_z1_lines, zero_vec1, color='tab:green', label='第１主成分 $z_1$')

# 直交射影点の描画
z1_flat = proj1.values[0]
zero_vec2 = np.zeros_like(proj1.values[0])
ax.plot(z1_flat, zero_vec2, 'o', color='green', ms=3.5)

# 施設名のテキストの表示
for s, x, y in zip(data1.values[:, 0], z1_flat, zero_vec2):
    ax.text(x=x, y=y+0.001, s=s, ha='center', fontsize=10)

# x1,x2の平均点の描画
ax.plot(X_mean[0], 0, '*', color='tab:orange', ms=10, label='平均')

# 修飾
ax.set(xlim=(13, 27), ylim=(-0.01, 0.01), yticks=[])
ax.grid(lw=0.5, alpha=0.5)
ax.legend(bbox_to_anchor=(1, 1));

【実行結果】

第１主成分に着目して、２次元平面から１次元直線の点に「要約」されました！
第１主成分の直線は、「移動後の各データ点とデータの中心・平均（★印）の距離」の合計が最も大きくなるように引かれています。

さらに１つの点に拡大して、細部を見ていきましょう！

「移動後の各データ点とデータの中心（平均）の距離」は、データのバラツキを意味して「分散」なのです。
分散最大化とは「移動後の各データ点とデータの中心との距離」を最大にする軸＝$${z_1}$$ を見つけることなのです。

続いて第２主成分です。
第２主成分の直線は、第１主成分の直線と直交する直線を検討します。
直交する直線にはたくさんの候補がありますが、やはり分散最大化の観点で「直行する直線への移動後の各データ点とデータの中心との距離」の合計が最も大きいものを選びます。

そして「固有値」は、主成分がデータのバラツキをどれだけ説明しているかを表す指標です。
固有値は、主成分分析の文脈で「説明された分散」（explained variance）と呼ばれることもあるようです。

固有値に着目して進みます。

■ 寄与率・累積寄与率
固有値の合計に対する「各主成分の固有値」の比率を寄与率と呼びます。
寄与率は主成分がデータを説明できている割合を示します。
「主成分の数」を絞り込むに当たり、例えば全体の７割の説明ができるようにするには、寄与率の累計である「累積寄与率」が７割に達するまでの主成分を１から順に選べばよいのです。

寄与率、累積寄与率を可視化しましょう。

## 寄与率、累積寄与率の可視化 p.104～

# 描画領域の設定
fig, ax = plt.subplots()
# 寄与率の棒グラフの描画
exp_var_ratio1_sscp = eig_val1_sscp / np.sum(eig_val1_sscp)
g = ax.bar(columns, exp_var_ratio1_sscp, width=0.4, alpha=0.5)
# 寄与率の値の表示
ax.bar_label(g, label_type='center', fmt='%.3f', fontsize=14)
# 累積寄与率の折れ線グラフの描画
cum_exp_var_ratio1_sscp = np.cumsum(eig_val1_sscp) / np.sum(eig_val1_sscp)
ax.plot(columns, cum_exp_var_ratio1_sscp, '-o')
# 累積寄与率 0.8 の水平点線の描画
ax.axhline(0.8, color='tab:red', ls='--')
# 修飾
ax.set_title('寄与率・累積寄与率');

【実行結果】

第１主成分だけでデータ全体の $${76.6\%}$$ を説明できます。
第２主成分を加えると $${100\%}$$ になります。
データを要約するために、説明割合のしきい値を $${70\%}$$ と置いたらなら、第１主成分を選択すればよいことになります。
これが「次元削減」であり、１つの新しい変数で「総合化」されたことを指します。

次は各データ点の要約後の指標に進みます。
第１主成分だけに絞り込まず、２つの主成分を見ます。

■ 主成分得点
テキストによると要約後の「新しい情報量＝主成分得点」です。
主成分得点はデータの中心（平均）からの距離です。
要約・次元削減後の新しいデータなのです！

要約・次元削減後の新しいデータである主成分得点を算出しましょう。
一般論では

「データ － データの平均」と「固有ベクトル」の「行列積」

で計算します。

SSCP行列に基づく場合には、上の計算結果をさらに「固有値 ÷ (N-1)」の正の平方根で割ります（ChatGPTによる）。

## 主成分得点の算出（想定上の情報損失量最小化からの導出結果） p.106～
# 中心化した変数のベクトルと固有ベクトルのドット積を、N-1でスケールした固有値で割る

score1_sscp = (
    (X - X.mean()).values @ eig_vec1_sscp.T / np.sqrt(eig_val1_sscp / (N-1)))
score1_sscp = pd.DataFrame(score1_sscp, columns=columns, index=area)
score1_sscp

【実行結果】
地域ごとの新しい情報量です。

２つの変数を１つに要約したい場合は、第１主成分（PC1）の主成分得点を使います。

ここで主成分分析を「要約」の場面で活用する例を考えてみます。

💡 要約の使い道１ 💡
変数が数百・数千あるようなデータを思い浮かべます。
表にすると列の数が数百・数千です。
さて、データの概観を可視化で確認しよう、と思った矢先…

「簡単に可視化できない！」

平面（２次元）、立体（３次元）が可視化によるデータ把握の限界でしょう。
変数が多い場合に２個・３個の変数で表現するには…
そうです！主成分分析で要約するのです！

数百・数千の変数を主成分分析にかけて、３個の新しい変数を得ます。
３次元でまとめて可視化するもよし、次の３つの平面もよし、です！
・第１主成分と第２主成分の平面チャート
・第１主成分と第３主成分の平面チャート
・第２主成分と第３主成分の平面チャート

では次の要約例の前振りに進みます。
主成分得点どうしの相関係数は０になります！（直交の効果）

# 主成分得点の各列の相関係数は０
score1_sscp.corr().round(10)

【実行結果】

💡 要約の使い道２ 💡
重回帰分析のときに問題視された「多重共線性」。
変数どうしが強い相関関係にある場合に引き起こされる現象でした。
「主成分得点」の相関係数が０になる性質を利用すると、「主成分得点」を重回帰分析の変数にすることで、多重共線性を回避できます。

要約の２話はこれにて終了です。

話は変わります。
先ほど、一般論では「データ － データの平均」と「固有ベクトル」の「行列積」で計算します、と書きました。
この方法で計算した主成分得点を「生の主成分得点」と呼びます（ChatGPTの受け売り）。
生の主成分得点を計算します。

## 生の主成分得点の検証　scikit-learnの主成分得点と一致
# ※中心化したデータ @ 固有ベクトル
 
score1_row_sscp = (X - X.mean()).values @ eig_vec1_sscp.T
score1_row_sscp = pd.DataFrame(score1_row_sscp, columns=columns, index=area)
score1_row_sscp

【実行結果】
生の主成分得点は scikit-learn の主成分得点と一致します。
（ただし、正負の符号が反転する場合があります）

■ 主成分負荷量
もう一つ、主成分分析の指標を検討します。
主成分負荷量です。loadings と呼ばれたりもします。
元の変数が主成分に寄与する程度を示す指標です。

主成分負荷量の定義にはいくつかのパターンが存在します。
（そして私はいつも混乱します）

ここでは、

固有ベクトル × 固有値の正の平方根 ÷ ばらつき調整量

とします。
SCCP行列を用いる場合、ばらつき調整量は「SSCP行列の対角成分」です。

## 主成分負荷量の算出（ChatGPTによる。テキスト未掲載のため、計算の妥当性は未確認）
# 主成分負荷量ij = 固有ベクトルij * √固有値i / √diag(SSCP行列jj), i=第i主成分,j=変数

loadings1_sscp = eig_vec1_sscp * (np.c_[eig_val1_sscp] / np.diag(S_sscp))**(1/2)
loadings1_sscp = pd.DataFrame(loadings1_sscp, columns=columns, index=X.columns)
loadings1_sscp

【実行結果】

■ バイプロット
「個々のデータ」の主成分の位置＝主成分得点と、主成分に対する「変数」の寄与度合い＝主成分負荷量を可視化するものが「バイプロット」です。
地域ごとの主成分得点、施設種類ごとの主成分負荷量を１つのチャートで表現しましょう。

## バイプロットの描画

# 描画領域の設定 ax:主成分得点の軸、twin:固有ベクトルの軸
fig, ax = plt.subplots(figsize=(6, 6))
twin = ax.twinx().twiny()

# 主成分得点の散布図の描画、施設名のテキストの表示
ax.scatter(score1_sscp['PC1'], score1_sscp['PC2'])
for s, (x, y) in zip(data1['地域名'].values, score1_sscp.values):
    ax.text(x=x+0.05, y=y, s=s)

# 主成分負荷量の散布図の描画、施設種類のテキストの表示
twin.scatter(loadings1_sscp['PC1'], loadings1_sscp['PC2'], color='tab:red')
for s, (x, y) in zip(X.columns, loadings1_sscp.values):
    twin.text(x=x, y=y+0.03, s=s, ha='right', va='bottom', color='tab:red',
              fontsize=12)
    twin.arrow(0, 0, x, y, color='tab:red')

# x=0, y=0の垂直線、水平線の描画
ax.axvline(0, color='black', lw=0.5, ls='--')
ax.axhline(0, color='black', lw=0.5, ls='--')

# 修飾
ax.set_xlabel('第１主成分', fontsize=12)
ax.set_ylabel('第２主成分', fontsize=12)
ax.set_title('バイプロット')
ax.set(xlim=(-2, 2), ylim=(-2, 2), xticks=range(-2, 3), yticks=range(-2, 3))
twin.set(xlim=(-1, 1), ylim=(-1, 1), xticks=range(-1, 2), yticks=range(-1, 2));

【実行結果】
青い点が地域ごとの主成分得点、赤い点＋線が施設種類ごとの主成分負荷量です。

ここで「第１主成分」と「第２主成分」の意味＝何を「総合化」したのか、を考えてみます。
$${x=0, y=0}$$ の線で区切られた４つの象限と赤い主成分負荷量の向き・長さに注目します。

◆

第１主成分に関しては、$${x=0}$$ の左側領域と右側領域に分けて考えます。
左側領域には、介護施設と医療施設の両方の主成分負荷量が含まれています。
右側領域には、主成分負荷量が含まれていません。
介護施設・医療施設は人口あたり施設数に関連する変数でした。

ですので、両方を含む左側領域は「介護・医療の充実度が高い地域」を示すと考えられます。
例えば地域 H は介護施設が 24、医療施設が 17 と両方とも多い地域です。
一方で右側領域は「介護・医療の充実度が低い地域」でしょう。
例えば地域 E は介護施設が 15、医療施設が 7 と両方とも多い地域です。

第１主成分は「地域の介護・医療の充実度」を総合的に示す指標と言えそうです。

◆

第２主成分に関しては、$${y=0}$$ の上側領域と下側領域に分けて考えます。
上側領域には医療施設、下側には介護施設が含まれています。
おそらく第２主成分は「地域の介護・医療の相対的な取り組み度」でしょう。

上側には医療施設が相対的に高い地域 H, D が含まれています。
下側には介護施設が相対的に高い地域 J, B が含まれています。
地域 H は介護施設も多いですが、相対的には医療施設も多い、という感じです。

まとめの可視化です。

■ ざっくりまとめ
主成分分析の手続きを見てきました！

• データから SSCP 行列を作る

• SSCP行列の固有値・固有ベクトルを求める

• 固有値・固有ベクトルに基づく「主成分」「寄与率」で主成分を理解する

• 固有値・固有ベクトルに基づく「主成分得点」「主成分負荷量」で分析する

お疲れ様でした！
一休みしましょう。

scikit-learn で主成分分析

scikit-learn は機械学習ライブラリです。
主成分分析用のクラス PCA を使っていきます！
こちらは「データの分散共分散行列」による主成分分析です。

■ 主成分分析の実行
２行で主成分分析の主要な指標を得ることができます！
PCA(・) でインスタンスを生成し、PCA.fit_transform(データ) でPCAの実行と主成分得点を取得します。

### sklearnのPCAで主成分分析

# PCAのインスタンスの生成
pca1 = PCA(n_components=2)

# PCAを実行して主成分得点を取得 shape=(標本サイズ, 主成分数)
score1 = pca1.fit_transform(X)

【実行結果】なし

■ 固有値・固有ベクトルの表示
PCAのインスタンス pca1 から固有値・固有ベクトルを取り出して、見やすい表に加工します。
固有値は「 .explained_variance_ 」で取り出せます。

# 固有値の表示
eig_val1 = pca1.explained_variance_
pd.DataFrame(eig_val1, index=columns, columns=['固有値']).T

【実行結果】
分散共分散行列を元にしていますので、SSCP行列ベースの値とは異なります。

固有ベクトルは「 .components_ 」で取り出せます。

# 固有ベクトルの表示
eig_vec1 = pca1.components_
pd.DataFrame(eig_vec1, index=columns, columns=X.columns).T

【実行結果】
分散共分散行列を元にしていますが、固有ベクトルは「符号が逆転している点を除けば」SSCP行列ベースの値と一致します！

■ 主成分
主成分 $${z_1, z_2}$$ は次のとおりです。

■ 寄与率・累積寄与率
SSCP行列ベースと固有値は異なりますが、率に変換した寄与率・累積寄与率は一致しています。
寄与率は 「 .explained_variance_ratio_」で取り出せます。

## 寄与率、累積寄与率の可視化

# 描画領域の設定
fig, ax = plt.subplots()
# 寄与率の棒グラフの描画
exp_var_ratio1 = pca1.explained_variance_ratio_
g = ax.bar(columns, exp_var_ratio1, width=0.4, alpha=0.5)
# 寄与率の値の表示
ax.bar_label(g, label_type='center', fmt='%.3f', fontsize=14)
# 累積寄与率の折れ線グラフの描画
cum_exp_var_ratio1 = np.cumsum(exp_var_ratio1)
ax.plot(columns, cum_exp_var_ratio1, '-o')
# 累積寄与率 0.8 の水平点線の描画
ax.axhline(0.8, color='tab:red', ls='--')
# 修飾
ax.set_title('寄与率・累積寄与率');

【実行結果】

■ 主成分得点
主成分得点は score1 に格納済みです。

「データ － データの平均」と「固有ベクトル」の「行列積」

で計算されています。

## 主成分得点の表示
pd.DataFrame(score1, columns=['PC1', 'PC2'], index=data1.index)

【実行結果】

「データ － データの平均」と「固有ベクトル」の「行列積」で主成分得点を計算してみましょう。

# 主成分得点の検証：中心化した変数のベクトルと固有ベクトルのドット積で計算する p.106
(X - X.mean()).values @ eig_vec1.T

【実行結果】
scikit-learn で計算された score と一致することが確認できました。

scilit-learn の主成分得点を SCCP行列ベースの主成分得点に変換しましょう。
scilit-learn の主成分得点を $${\sqrt{\text{固有値}}}$$ で割ります。

# 主成分得点の検証：SSCP行列の結果に合わせる ※√(固有値)で割る
score1 / np.sqrt(eig_val1)

【実行結果】
「符号が逆転している点を除けば」SSCP行列ベースの値と一致します！

■ 主成分負荷量
主成分負荷量の取得には計算が必要です。

固有ベクトル × 固有値の正の平方根 ÷ 標準偏差

で計算します。

## 主成分負荷量の算出
# 主成分負荷量ij = 固有値iの平方根 * 固有ベクトルij / 標準偏差j,  i=第i主成分,j=変数

loadings1 = np.c_[eig_val1**(1/2)] * eig_vec1 / X.std(ddof=1).values
pd.DataFrame(loadings1, columns=columns, index=X.columns)

【実行結果】
「符号が逆転している点を除けば」SSCP行列ベースの値と一致します！

主成分負荷量は、主成分 $${z}$$ と元の変数 $${X}$$ との相関係数です。
計算して、scikit-learn の主成分負荷量と比べてみましょう。

# 主成分負荷量の検証：主成分zと元の変数との相関係数で計算 shape=(主成分数, 変数数)

# 主成分zの算出
z = pca1.components_ @ X.values.T

# 主成分zと変数の相関係数を算出
loadings1 = np.array([np.corrcoef(z[i], X.iloc[:, j])[0, 1]
                      for i in range(2) for j in range(2)]
                    ).reshape((2, 2))

# 結果の表示
print('主成分負荷量の検証:')
print(loadings1)

【実行結果】
確かに一致しました！

■ バイプロット
scikit-learn 主成分得点と主成分負荷量を用いてバイプロットを描きます。

## バイプロットの描画

# 描画領域の設定 ax:主成分得点の軸、twin:固有ベクトルの軸
fig, ax = plt.subplots(figsize=(6, 6))
twin = ax.twinx().twiny()

# 主成分得点の散布図の描画、施設名のテキストの表示
ax.scatter(score1[:, 0], score1[:, 1])
for s, (x, y) in zip(data1['地域名'].values, score1):
    ax.text(x=x+0.2, y=y, s=s)

# 主成分負荷量の散布図の描画、施設種類のテキストの表示
twin.scatter(loadings1[:, 0], loadings1[:, 1], color='tab:red')
for s, (x, y) in zip(X.columns, loadings1):
    twin.text(x=x, y=y-0.1, s=s, ha='left', color='tab:red', fontsize=12)
    twin.arrow(0, 0, x, y, color='tab:red')

# x=0, y=0の垂直線、水平線の描画
ax.axvline(0, color='black', lw=0.5, ls='--')
ax.axhline(0, color='black', lw=0.5, ls='--')

# 修飾
ax.set_xlabel('第１主成分', fontsize=12)
ax.set_ylabel('第２主成分', fontsize=12)
ax.set_title('バイプロット')
ax.set(xlim=(-8, 8), ylim=(-8, 8),
       xticks=range(-8, 9, 2), yticks=range(-8, 9, 2))
twin.set(xlim=(-1, 1), ylim=(-1, 1),
         xticks=np.arange(-1, 1.1, 0.5), yticks=np.arange(-1, 1.1, 0.5));

【実行結果】
SSCP行列ベースと比べて指標の $${\pm}$$ 符号が逆転しているので、バイプロットも左右逆店・上下逆転しています。
４象限の位置は違いますが、意味合いはSSCP行列と同様の結果を示しています。

駆け足になりました🏃🏃‍♀️🏃‍♂️
最後の【「データの分散共分散行列」で主成分分析】で計算手順を見ていきますので、そのときにイメージアップしましょう！

statsmodels で主成分分析

statsmodels は統計解析ライブラリです。
主成分分析用のクラス PCA を使っていきます！
scikit-learn と名前 PCA がかぶるので、statmodels の方は sm_PCA とします。
こちらは「データの SSCP 行列」による主成分分析です。

■ 主成分分析の実行
１行で主成分分析の主要な指標を得ることができます！

### statsmodelsのpcaで主成分分析

# PCAの実行
pca1_sm = sm_PCA(X, standardize=False)

【実行結果】なし

■ 固有値・固有ベクトルの表示
PCAのインスタンス pca1_sm から固有値・固有ベクトルを取り出して、見やすい表に加工します。
固有値は「 .eigenvals 」で取り出せます。

# 固有値の表示
eig_val1_sm = pca1_sm.eigenvals.values
pd.DataFrame(eig_val1_sm, index=columns, columns=['固有値']).T

【実行結果】
SSCP 行列ベースの固有値と一致します。

scikit-learn の固有値に変換してみましょう。
statsmodels の固有値を $${N-1}$$ で割ります。

# 固有値の検証：scikit-learnに合わせる
N = X.shape[0]  # 標本サイズ
print('固有値:')
print(eig_val1_sm / (N-1))

【実行結果】

固有ベクトルは「 .eigenvecs 」で取り出せます。

# 固有ベクトルの表示
eig_vec1_sm = pca1_sm.eigenvecs
eig_vec1_sm.columns = columns
eig_vec1_sm.index = X.columns
eig_vec1_sm

【実行結果】
「符号が逆転している点を除けば」SSCP行列ベースの値と一致します！
scikit-learn の固有ベクトルとは同じ符号・同じ値です。

■ 主成分
主成分 $${z_1, z_2}$$ は次のとおりです。

■ 寄与率・累積寄与率

## 寄与率、累積寄与率の可視化

# 描画領域の設定
fig, ax = plt.subplots()
# 寄与率の棒グラフの描画
exp_var_ratio1_sm = eig_val1_sm / np.sum(eig_val1_sm)
g = ax.bar(columns, exp_var_ratio1_sm, width=0.4, alpha=0.5)
# 寄与率の値の表示
ax.bar_label(g, label_type='center', fmt='%.3f', fontsize=14)
# 累積寄与率の折れ線グラフの描画
cum_exp_var_ratio1_sm = np.cumsum(eig_val1_sm) / np.sum(eig_val1_sm)
ax.plot(columns, cum_exp_var_ratio1_sm, '-o')
# 累積寄与率 0.8 の水平点線の描画
ax.axhline(0.8, color='tab:red', ls='--')
# 修飾
ax.set_title('寄与率・累積寄与率');

【実行結果】
SSCP 行列ベース、scikit-learn の両方と同じです。

ちなみに statsmodels は、寄与率・累積寄与率の可視化手段の１つ「スクリープロット」を簡単に描画できます。

# スクリープロットの描画
pca1_sm.plot_scree(log_scale=False);

【実行結果】
主成分が２つしか無いので寂しいプロットになりました。。。
ちなみにこちらは固有値そのものを描画しています。

■ 主成分得点
固有ベクトルは「 .scores 」で取り出せます。

「（データ － データの平均）と固有ベクトルの行列積」$${\div \ \sqrt{\text{固有値}}}$$

で計算されています。

## 主成分得点の表示
score1_sm = pca1_sm.scores
score1_sm.columns = columns
score1_sm.index = area
score1_sm

【実行結果】
SSCP 行列ベースとも、scikit-learn とも異なる結果になりました。

statsmodels の主成分得点に $${\times \sqrt{\text{固有値}}}$$ で、scikit-learn の主成分得点に変換できます。

# 主成分得点の検証：sckit-learnの結果に合わせる ※√(固有値)を掛ける
score1_sm * np.sqrt(eig_val1_sm)

【実行結果】

statsmodels の主成分得点に $${\times \sqrt{N-1}}$$ で、SSCP 行列ベースの主成分得点に変換できます。

# 主成分得点の検証：SSCP行列ベースの結果に合わせる ※√(N-1)を掛ける
score1_sm * np.sqrt(N - 1)

【実行結果】

■ 主成分負荷量
主成分負荷量の取得には計算が必要です。

固有ベクトル × 固有値の正の平方根 ÷ SSCP行列の対角成分の正の平方根

で計算します。

## 主成分負荷量の算出：ChatGPTによる。テキスト未掲載のため、この計算の妥当性は未確認
# 主成分負荷量ij = 固有ベクトルij * √固有値i / √diag(SSCP行列jj), i=第i主成分,j=変数

loadings1_sm = eig_vec1_sm * (np.c_[eig_val1_sm] / np.diag(S_sscp))**(1/2)
loadings1_sm

【実行結果】
scikit-learn の主成分負荷量と一致します。
「符号が逆転している点を除けば」SSCP行列ベースの値と一致します！

■ バイプロット

## バイプロットの描画

# 描画領域の設定 ax:主成分得点の軸、twin:固有ベクトルの軸
fig, ax = plt.subplots(figsize=(6, 6))
twin = ax.twinx().twiny()

# 主成分得点の散布図の描画、施設名のテキストの表示
ax.scatter(score1_sm.iloc[:, 0], score1_sm.iloc[:, 1])
for s, (x, y) in zip(data1['地域名'].values, score1_sm.values):
    ax.text(x=x+0.02, y=y, s=s)

# 主成分負荷量の散布図の描画、施設種類のテキストの表示
twin.scatter(loadings1_sm.iloc[:, 0], loadings1_sm.iloc[:, 1],
             color='tab:red')
for s, (x, y) in zip(X.columns, loadings1_sm.values):
    twin.text(x=x*1.05, y=y, s=s, ha='left', color='tab:red', fontsize=12)
    twin.arrow(0, 0, x, y, color='tab:red')

# x=0, y=0の垂直線、水平線の描画
ax.axvline(0, color='black', lw=0.5, ls='--')
ax.axhline(0, color='black', lw=0.5, ls='--')

# 修飾
ax.set_xlabel('第１主成分', fontsize=12)
ax.set_ylabel('第２主成分', fontsize=12)
ax.set_title('バイプロット')
ax.set(xlim=(-0.7, 0.7), ylim=(-0.7, 0.7))
twin.set(xlim=(-1.1, 1.1), ylim=(-1.1, 1.1),
         xticks=np.arange(-1, 1.5, 0.5), yticks=np.arange(-1, 1.5, 0.5));

【実行結果】
介護施設・医療施設の第１主成分の主成分負荷量が正の値であり、かつ、各地域の主成分得点の配置バランスが良いので、私はこの statsmodels のバイプロットが一番見やすいです。

statsmodels の PCA の特徴を眺めることができました。
お疲れ様です。

「データの分散共分散行列」で主成分分析

データの分散共分散行列を用いた主成分分析について、「計算手順」を追って確認しましょう。

■ データの分散共分散行列の算出
いきなりまとめです！

高濃度な情報から固有値・固有ベクトルを取り出します。

ところで変数行列 X は pandas データフレームです。
データフレームのメソッド「 .cov() 」で分散共分散行列 S_cov を算出しましょう。

## データの分散共分散行列の算出

S_cov = X.cov(ddof=1)
print('分散共分散行列:')
display(S_cov.round(3))

【実行結果】

ちなみに…
SSCP 行列を $${N-1}$$ で割ると、分散共分散行列になります！
計算してみましょう。

# SSCP行列をN-1で割ると分散共分散行列になる
(S_sscp / (N-1)).round(3)

【実行結果】
はい！このとおりです！
SSCP 行列の偏差平方和は分散の素、偏差積和は共分散の素だったのです！

■ 分散共分散行列の固有値・固有ベクトルの算出
分散共分散行列 S_cov について、numpy.linalg の eig() を用いて固有値・固有ベクトルを算出します。
あわせて固有値の大きい順に固有値・固有ベクトルを並び替えます。

## 分散共分散行列の固有値・固有ベクトルの算出 p.112～

# 固有値・固有ベクトルの算出
eig_val1_hand, eig_vec1_hand = np.linalg.eig(X.cov().values)

# 固有値の降順ソートで並べ替え
sort_idx = eig_val1_hand.argsort()[::-1]
eig_val1_hand = eig_val1_hand[sort_idx]
eig_vec1_hand = eig_vec1_hand[:, sort_idx]

【実行結果】なし

固有値を表示しましょう。

# 固有値の表示
pd.DataFrame(eig_val1_hand, index=columns, columns=['固有値']).T

【実行結果】
scikit-learn の固有値と一致します。

固有ベクトルを表示しましょう。

# 固有ベクトルの表示 ※SSCP行列ベースの値と比べて正負が逆転
pd.DataFrame(eig_vec1_hand, columns=columns, index=X.columns)

【実行結果】
scikit-learn の固有ベクトルとは、符号は逆転していますが、絶対値は一致しています。
しかも SSCP 行列ベースの固有ベクトルと一致しています。

■ 主成分
主成分 $${z_1, z_2}$$ は次のとおりです。

■ 寄与率・累積寄与率

## 寄与率、累積寄与率の可視化

# 描画領域の設定
fig, ax = plt.subplots()
# 寄与率の棒グラフの描画
exp_var_ratio1_hand = eig_val1_hand / np.sum(eig_val1_hand)
g = ax.bar(columns, exp_var_ratio1_hand, width=0.4, alpha=0.5)
# 寄与率の値の表示
ax.bar_label(g, label_type='center', fmt='%.3f', fontsize=14)
# 累積寄与率の折れ線グラフの描画
cum_exp_var_ratio1_hand = np.cumsum(eig_val1_hand) / np.sum(eig_val1_hand)
ax.plot(columns, cum_exp_var_ratio1_hand, '-o')
# 累積寄与率 0.8 の水平点線の描画
ax.axhline(0.8, color='tab:red', ls='--')
# 修飾
ax.set_title('寄与率・累積寄与率');

【実行結果】
すべての方法で寄与率・累積寄与率は一致しています。

■ 主成分得点
主成分得点は

「データ － データの平均」と「固有ベクトル」の「行列積」

で計算します。

## 主成分得点の算出：中心化した変数のベクトルと固有ベクトルのドット積で計算する p.106

score1_hand = ((X - X.mean()).values @ eig_vec1_hand.T)
pd.DataFrame(score1_hand, columns=columns, index=area)

【実行結果】
scikit-learn の主成分得点とは、符号は逆転していますが、絶対値は一致しています。

■ 主成分負荷量
主成分負荷量は

固有ベクトル × 固有値の正の平方根 ÷ 標準偏差

で計算します。

## 主成分負荷量の算出
# 主成分負荷量ij = 固有ベクトルij * 固有値iの平方根 / 標準偏差j,  i=第i主成分,j=変数

loadings1_hand = (
    eig_vec1_hand * np.c_[eig_val1_hand**(1/2)] / X.std(ddof=1).values
)
pd.DataFrame(loadings1_hand, columns=columns, index=X.columns)

【実行結果】
scikit-learn の主成分負荷量とは、符号は逆転していますが、絶対値は一致しています。

■ バイプロット

## バイプロットの描画

# 描画領域の設定 ax:主成分得点の軸、twin:固有ベクトルの軸
fig, ax = plt.subplots(figsize=(6, 6))
twin = ax.twinx().twiny()

# 主成分得点の散布図の描画、施設名のテキストの表示
ax.scatter(score1_hand[:, 0], score1_hand[:, 1])
for s, (x, y) in zip(data1['地域名'].values, score1_hand):
    ax.text(x=x+0.05, y=y, s=s)

# 主成分負荷量の散布図の描画、施設種類のテキストの表示
twin.scatter(loadings1_hand[:, 0], loadings1_hand[:, 1], color='tab:red')
for s, (x, y) in zip(X.columns, loadings1_hand):
    twin.text(x=x, y=y*1.1, s=s, ha='left', color='tab:red', fontsize=12)
    twin.arrow(0, 0, x, y, color='tab:red')

# x=0, y=0の垂直線、水平線の描画
ax.axvline(0, color='black', lw=0.5, ls='--')
ax.axhline(0, color='black', lw=0.5, ls='--')

# 修飾
ax.set_xlabel('第１主成分', fontsize=12)
ax.set_ylabel('第２主成分', fontsize=12)
ax.set_title('バイプロット')
ax.set(xlim=(-8, 8), ylim=(-8, 8),
       xticks=range(-8, 9, 2), yticks=range(-8, 9, 2))
twin.set(xlim=(-1.1, 1.1), ylim=(-1.1, 1.1),
         xticks=np.arange(-1, 2, 0.5), yticks=np.arange(-1, 2, 0.5));

【実行結果】
scikit-learn のバイプロットの左右反転、上下反転したものになりました。
固有ベクトルの符号がscikit-learn と逆転しているためです。
４象限の位置は異なりますが、解釈・意味合いは scikit-learn と同等です。

ここまでの長旅、お疲れ様でした！
SSCP行列と分散共分散行列の違いで主成分分析が影響を受けることや、Python ライブラリの計算過程が少し見えてきましたか？

おわりに～「情報損失量最小化」で主成分分析

ここまでは「分散の最大化」に基づく主成分分析を行ってきました。
ところで、テキスト Section 3.3 ～ 3.6 は「情報損失量最小化」に基づく主成分分析を取り扱っています。

情報損失量とは要約の過程で失う情報量のことであり、上の図の破線部分に相当します。
情報損失量の最小化を解いて固有値・固有ベクトルを求めます。
ちなみに「分散最大化」は上の図の「新しい情報量」の最大化を解いたものです。
最小化・最大化問題の対象が双方異なる点に留意しましょう。

ChatGPTによると、情報損失量の最小化で扱う行列は、SSCP 行列から導出できるようです。

情報損失量最小化では、行列 $${M}$$ の固有値・固有ベクトルを求めます。
$${2 \times 2}$$ の SSCP 行列の場合、SSCP 行列の対角成分以外の成分について正負を逆転すると行列 $${M}$$ になるようです。

■ 行列 $${M}$$ の算出

## テキストの行列の作成

# SSCP行列を変形してテキストの行列を作成
matrix_text = S_sscp * np.array([[1, -1], [-1, 1]])

# 結果の表示
print('テキストの行列:')
matrix_text

【実行結果】

■ 固有値・固有ベクトルの算出
上記の行列について、固有値・固有ベクトルを算出します。

## テキストの行列の固有値・固有ベクトルの算出 p.100

# 固有値・固有ベクトルの算出
eig_val1_text, eig_vec1_text = np.linalg.eig(matrix_text)

# 固有値の降順ソートでindex取得
sort_idx = eig_val1_text.argsort()
eig_val1_text = eig_val1_text[sort_idx]
eig_vec1_text = eig_vec1_text[:, sort_idx][::-1]  # 行の入れ替え

【留意点】

• 固有値の小さい順に並び替えました。
  分散最大化のときは固有値の大きい順に並び替えました。

• 固有ベクトルの行の順番を入れ替えました。
  分散最大化のときは行の順番は変えていません。

固有値を表示します。

# 固有値の表示
pd.DataFrame(eig_val1_text, index=columns, columns=['固有値']).T

【実行結果】
並び替え順が異なるものの、SSCP 行列と同じ値です。

「第１主成分の固有値が小さい」

このことが波乱の幕開けとなる予感です…

固有ベクトルを表示します。

# 固有ベクトルの表示 ※テキストと比べて正負が逆転
pd.DataFrame(eig_vec1_text, columns=columns, index=X.columns)

【実行結果】
こちらは SSCP 行列の固有ベクトルと順序・値が一致しました。

■ 主成分
主成分 $${z_1, z_2}$$ は次のとおりです。

$${z_1, z_2}$$ をデータの散布図上で示してみましょう。

## 主成分z1, z2の描画 p.110 図 3.8.1 の改造

## z1, z2の傾き・切片の算出 ※横軸x1, 縦軸x2のグラフ描画用
# x1, x2の平均値の算出
x1_mean, x2_mean = X.mean(axis=0).values
# z1,z2の傾きの算出
slope_z1 = eig_vec1_text[1, 0] / eig_vec1_text[0, 0]
slope_z2 = eig_vec1_text[1, 1] / eig_vec1_text[0, 1]
# z1,z2の切片の算出 ※z1,z2はx1,x2の平均を通る
intercept_z1 = x2_mean - slope_z1 * x1_mean
intercept_z2 = x2_mean - slope_z2 * x1_mean

## z1, z2の直線の算出
# z1のx軸:x1,y軸:x2の値
x1_z1_lines = np.linspace(16.5, 25, 2)
x2_z1_lines = slope_z1 * x1_z1_lines + intercept_z1
# z2のx軸:x1,y軸:x2の値
x1_z2_lines = np.linspace(16, 24.3, 2)
x2_z2_lines = slope_z2 * x1_z2_lines + intercept_z2

## 描画
# 描画領域の設定
fig, ax = plt.subplots(figsize=(6, 6))

# z1、z2の直線の描画
ax.plot(x1_z1_lines, x2_z1_lines, color='tab:green', label='第１主成分 $z_1$')
ax.plot(x1_z2_lines, x2_z2_lines, color='tab:red', label='第２主成分 $z_2$')

# 散布図の描画
sns.scatterplot(data=data1, x='介護施設', y='医療施設', s=70, alpha=0.7, ax=ax)
# 施設名のテキストの表示
for (s, x, y) in data1.values:
    ax.text(x=x, y=y+0.2, s=s, fontsize=12)

# x1,x2の平均点の描画
ax.plot([x1_mean], [x2_mean], '*', color='tab:orange', ms=15)
ax.text(x=x1_mean+0.5, y=x2_mean, va='center', s=(x1_mean, x2_mean))

# 修飾
ax.set(xlim=(13, 27), ylim=(4, 19), aspect='equal')
ax.grid(lw=0.5, alpha=0.5);

【実行結果】
SSCP 行列ベースのときと全く同じになりました。

■ これからのことを考える

さて「寄与率を求めよう」と考えたとき

• 情報損失量最小化ベースの固有値は「失った情報量」に基づいており、「新しい情報量」の説明度合いを示していない

• 主成分が２つのみの場合は、テキスト p.105 のように、固有値合計から固有値を引いた差分を「新しい情報量」とみなすことは可能である

• しかし、主成分が３つ以上ある場合に、新しい情報量の指標を固有値で示すことが難しい

• これらに加えて、固有ベクトルの行の順番を入れ替えたけど、主成分が３つ以上ある場合はどのように並び替えすればいいのだろう…？

みたいなことで頭がループ状態になりました。

なんか、情報損失量最小化ベースは使いづらいのかな…

こう思ってしまい、深堀りするのを諦めました。
しかも他の書籍やWeb情報は「分散最大化」を前提としており、「情報損失量最小化」ベースの説明を見つけられませんでした。

多数派（っぽい）分散最大化にのっかっておこう

これです！

ただし、次回・次々回はテキストに準拠して「情報損失量最小化」に取り組みます。

おあとがよろしいようで…

記事の最後はChatGPTが締めくくります。
今回は窓越しの風景に想いを寄せて。

📘 ChatGPTのひとこと：

今回の写経は以上です。

【追補】主成分負荷量（2026/03/27追記）

入力行列・計算式・計算結果の意味を Gemini といっしょに整理したので、こちらで共有いたします。

◆ 主成分負荷量の定義ざっくりマトリクス

$$
\begin{array}{lclc}
入力行列 & 計算式 & 主成分負荷量の意味 & 数値の範囲 \\
\hline
\\
相関行列 & v_{ij} \sqrt{\lambda_j} & 変数 x_i と 主成分得点 z_j の & -1 〜 1 \\
&&「相関係数」& \\
\\
共分散行列 & v_{ij} \sqrt{\lambda_j} & 変数 x_i と 主成分得点 z_j の & 制限なし \\
&&「共分散」& \\
\\
共分散行列 & \cfrac{v_{ij} \sqrt{\lambda_j}}{\sigma_i} & 変数 x_i と 主成分得点 z_j の  & -1 〜 1 \\
&&「相関係数」& \\
\\
\text{SSCP}行列 & v_{ij} \sqrt{\lambda_j} & 変数 x_i と 主成分得点 z_j の & 制限なし \\
&&「偏差積和」& \\
\end{array}
$$

（留意事項）掲載を省略したケース

• 「固有ベクトルを主成分負荷量だとする」ケース

◆ 重要ポイント

1. 「負荷量」は橋渡し役

主成分負荷量は、「元のデータの世界（変数）」 と 「新しいデータの世界（主成分得点）」 を繋ぐ、いわば「親密度」の指標です。

• 相関係数としての負荷量: 「パターンがどれくらい似ているか」を測る。

• 共分散としての負荷量: 「連動して動くエネルギーの勢い」を測る。

• 偏差積和としての負荷量: 「連動して動くエネルギーの総量」を測る。

2. なぜ「主成分」と略されるのか

多くの文脈で「変数と主成分の相関」と略されるのは、主成分得点がその主成分（軸）上での位置を完璧に代表しているからです。概念としては「軸」を指し、計算としては「得点」を指している、と解釈すると非常にスッキリします。

【追補】入力行列と固有値の関係（2026/03/27追記）

Gemini が整理してくれた表を共有いたします。

$$
\begin{array}{lll}
入力行列 & 固有値の合計 (\sum \lambda) \\
\hline
\\
相関行列 & 変数の数 (p)  \\
\\
共分散行列 & 全分散の和 (\sum \sigma^2)  \\
\\
\text{SSCP}行列 & 偏差平方和の総計  \\
\end{array}
$$

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！

１．のんびり統計

統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。

２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で

書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！

３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で

書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！

４．楽しい写経 ベイズ・Python等

ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀

５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で

書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。

６．データサイエンスっぽいことを綴る

統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

７．Python機械学習プログラミング実践記

書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。