機械学習のPython実現ガイド

機械学習の用語

基本概念

• 教師あり学習・教師なし学習：正解データ（教師信号、ラベル）の有無による学習方法の分類である。教師あり学習では入力と正解のペアからパターンを学習する。教師なし学習では正解データを用いずにデータの構造やパターンを発見する。
• 分類・回帰：予測問題の種類である。分類は離散的なカテゴリ（クラス）への割り当て、回帰は連続的な数値の予測である。
• 訓練データ：モデルのパラメータを学習するために使用するデータセットである。入力（特徴量）と出力のペアで構成される。
• 特徴量：データの特性や属性を数値化した変数である。画像認識ではピクセル値や色情報、テキスト分析では単語の出現頻度などが用いられる。PCAなどの手法により元の特徴量から新たな特徴量を生成することもある。
• 距離：データ点間の類似度を測る尺度である。ユークリッド距離、マンハッタン距離などの定義がある。距離が小さいほどデータ点は類似している。
• 確率分布：確率変数がとりうる値とその確率の関係を表す数学的モデルである。正規分布、ベルヌーイ分布、ポアソン分布などがある。
• 分散：データが平均値からどれだけ散らばっているかを表す統計量である。分散が大きいほどばらつきが大きい。
• 相関：2つの変数間の関連性の強さと方向を示す統計量である。相関係数は-1から1の値をとり、正の相関では一方が増えると他方も増える傾向を示す。
• 次元：データを表現するために必要な特徴量（変数）の数である。
• 最適化：目的関数（損失関数や評価関数）を最小化または最大化するパラメータを探索する手法である。勾配降下法などのアルゴリズムが用いられる。
• 陽性・陰性（正例・負例）：二値分類における2つのクラスの呼称である。検出対象の事象を陽性とすることが多いが、定義は問題設定による。

教師あり学習手法

• ロジスティック回帰：線形モデルとシグモイド関数を組み合わせた分類手法である。多クラス分類にはソフトマックス関数を適用する。最尤推定により重みパラメータを最適化する。
• サポートベクターマシン（SVM）：マージン最大化により分離超平面を求める手法である。カーネルトリック（高次元空間での内積を元の空間での計算により求める手法）により、データを高次元空間に写像した非線形分類を実現する。
• k近傍法（k-NN）：新規データに対し、訓練データからk個の最近傍点を距離計算で特定し、近傍点の値に基づいて予測する手法である。明示的なモデルを構築せず訓練データを保持するインスタンスベース学習（遅延学習：訓練時には学習を行わず、予測時に計算を行う方式）に分類される。
• 決定木：特徴量に基づく階層的な木構造により、データを再帰的に分割して判断ルールを生成する手法である。
• ランダムフォレスト：複数の決定木を組み合わせたアンサンブル学習手法である。アンサンブル学習とは、複数のモデルの予測を統合して単一モデルより予測性能の向上を図る手法である。
• 線形判別分析（LDA）：分類と次元削減の両機能を持つ教師あり学習手法である。クラス間分散を最大化し、クラス内分散を最小化する軸を求めて特徴抽出を行う。kクラスでは最大でk-1個の判別軸が得られる。

教師なし学習手法（クラスタリング）

• k-meansクラスタリング：クラスタ中心の初期配置後、データ点の割り当てと中心の再計算を収束まで反復する手法である。初期値に依存するため、n_init で指定した回数だけ異なる初期値で実行し、クラスタ内平方和が最小の結果を採用する。
• 階層的クラスタリング：データの階層構造を抽出する手法である。ボトムアップ（凝集型）は個別データ点から開始して近接クラスタを段階的に統合し、トップダウン（分割型）は単一クラスタから段階的に分割する。
• DBSCAN：密度ベースのクラスタリング手法である。データ点を、近傍に一定数以上のデータ点を持つコアポイント、コアポイントの近傍にある境界ポイント、どのクラスタにも属さないノイズポイントの3種類に分類する。クラスタ数の事前指定が不要で、外れ値を自動識別できる。結果はeps（近傍の半径）とmin_samples（コアポイント判定に必要な近傍点数）に依存する。

次元削減・特徴抽出手法

• PCA（主成分分析）：データの分散が最大となる軸を算出して次元削減を行う統計的手法である。直交する主成分を順次抽出し、各主成分の寄与率で重要度を評価できる。線形変換に基づくため、非線形な構造の抽出には適さない。
• 因子分析：観測変数の背後にある潜在的な共通因子を抽出する統計的手法である。PCAがデータの分散を最大化する軸を探索するのに対し、因子分析は観測変数間の相関構造を説明する潜在因子を仮定する。心理学や社会科学で用いられる。scikit-learnの実装は既定では因子回転を行わないため、解釈しやすい負荷量を得るには rotation='varimax' 等を指定する。
• t-SNE：確率分布の一致度を最適化する非線形次元削減手法である。近傍点間の関係（局所構造）を保持しながら、クラスタ構造の可視化に適した低次元表現を生成する。サンプル数の増加に伴って計算時間が増えるため、大規模データではPCAで事前に次元を削減してから適用する。乱数を用いるため、random_state を固定しない限り実行ごとに結果が異なる。

行列分解手法

• 特異値分解（SVD）：行列を3つの行列（U、Σ、V^T）の積に分解する手法である。特異値の大きさは各成分の重要度を示す。データ圧縮や潜在意味解析に応用される。
• QR分解：行列を直交性を持つQと上三角行列Rの積に分解する手法である。正方行列の完全QR分解ではQは直交行列であり、長方行列のreduced QR分解ではQの列が互いに直交する。最小二乗法や固有値計算で用いられる。

評価指標

• 精度（Accuracy）：全予測における正解の割合である。
• 適合率（Precision）：陽性予測における正解の割合である。
• 再現率（Recall）：実際の陽性サンプルのうち正しく検出された割合である。
• F1スコア：適合率と再現率の調和平均である。
• 適合率・再現率・F1スコアは二値分類の定義である。多クラス分類では各クラスについて算出し、平均（マクロ平均・加重平均など）をとる。classification_report は各クラスの値と平均を一括出力する。
• コフェネティック相関係数：階層的クラスタリングの評価指標である。元の距離とデンドログラム上の距離の相関を測定する。
• KL-divergence（Kullback-Leibler情報量）：2つの確率分布間の差異を測定する指標である。t-SNEでは高次元空間と低次元空間の確率分布の一致度評価に使用される。
• 寄与率：各主成分が全体の分散を説明する割合である。PCAにおける主成分の重要度を表す。
• 累積寄与率：複数の主成分による累積的な分散説明率である。次元削減後に保持される分散の割合を評価する。
• 固有値：各主成分が持つ分散の大きさである。PCAにおいて主成分の重要度を示す。
• 主成分負荷量：各主成分と元の変数との相関である。主成分の解釈に用いられる。
• 共通性：因子分析において、観測変数が共通因子によって説明される分散の割合である。
• 判別係数：線形判別分析（LDA）において、各クラスを判別するための線形結合の係数である。

実装のための技術

• 標準化：各特徴量の平均を0、標準偏差を1に変換する前処理である。特徴量間のスケールを統一し、距離、内積、分散、共分散、正則化の影響を受けるアルゴリズムで特定の特徴量の影響が過大になることを防ぐ。決定木やランダムフォレストはスケールに依存しにくいため標準化は不要である。
• 次元削減：高次元データを低次元空間に変換する手法である。計算コストの削減、可視化、過学習の抑制を目的とする。
• k分割交差検証：データセットをk個に分割し、1つをテストデータ、残りを訓練データとして使用する処理をk回繰り返す手法である。各分割で得られた評価指標の平均と標準偏差でモデルの汎化性能を評価する。
• グリッドサーチ：ハイパーパラメータの候補値の全組み合わせを試行し、交差検証で最適なパラメータセットを決定する手法である。
• データの可視化：データの分布、パターン、クラスタ構造、予測結果を図示する手法である。
• エルボー法：k-meansクラスタリングにおいて、クラスタ数kとクラスタ内平方和（inertia）の関係をプロットし、値の減少が緩やかになる転換点から最適なクラスタ数を判断する手法である。
• デンドログラム：階層的クラスタリングの結果を樹形図として可視化する手法である。データ点間の距離関係と各段階での統合距離を示す。
• 最尤推定：観測データが得られる確率を最大化するパラメータを推定する統計的手法である。ロジスティック回帰や因子分析で用いられる。

実装結果の検証・評価

• 汎化性能：モデルが未知データに対して正確な予測を行う能力である。交差検証で評価する。
• ハイパーパラメータ：学習前に設定する、学習アルゴリズムの動作を制御するパラメータである。k-NNのk値、SVMのカーネル種類、決定木の最大深度などが該当する。
• 混同行列（Confusion Matrix）：分類モデルの予測結果を真陽性、偽陽性、真陰性、偽陰性に分類して表形式で示す評価ツールである。精度、適合率、再現率、F1スコアの算出基礎となる。
• シルエット係数（Silhouette Coefficient）：クラスタリング結果の品質を評価する指標である。各データ点について、同一クラスタ内の平均距離と最近隣クラスタへの平均距離から算出し、-1から1の範囲で評価する。1に近いほどクラスタが明確に分離されている。

Pythonでの実現

演習

以下の演習では、Irisデータセット（4特徴量、3クラス、150サンプルから構成される分類用の標準データセット）を使用する。スケールの影響を受けるアルゴリズムでは標準化を行う。教師あり学習で標準化を含む交差検証やグリッドサーチを行う場合は Pipeline を用い、検証用データの情報が標準化の学習に混入しないようにする。分類問題の訓練・テスト分割では stratify=y を指定し、訓練データとテストデータのクラス比率を保つ。決定木およびランダムフォレストは特徴量のスケールに依存しにくいため標準化を行わない。

演習1．ロジスティック回帰

from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import classification_report, confusion_matrix

# Irisデータセットの読み込み
iris = load_iris()
X = iris.data
y = iris.target

# データの分割（分類問題なのでクラス比率を保つ）
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42, stratify=y
)

# 標準化とモデルをPipelineで一体化
model = Pipeline([
    ('scaler', StandardScaler()),
    ('logistic_regression', LogisticRegression(max_iter=1000, random_state=42))
])

# モデルの構築と学習
model.fit(X_train, y_train)

# テストデータでの予測と性能評価
y_pred = model.predict(X_test)
print(classification_report(y_test, y_pred, target_names=iris.target_names))

# 混同行列の表示
print("Confusion Matrix:")
print(confusion_matrix(y_test, y_pred))

# 5分割交差検証による性能評価
# Pipelineを用いるため、各訓練分割の中だけで標準化が学習される
cv_scores = cross_val_score(model, X, y, cv=5)
print(f"Cross-validation scores: {cv_scores}")
print(f"Average CV score: {cv_scores.mean():.3f} (+/- {cv_scores.std() * 2:.3f})")

線形モデルとシグモイド関数を組み合わせた分類手法である。多クラス分類では各クラスへの所属確率に基づいて予測する。LogisticRegression は既定で正則化を用いるため、特徴量のスケールを揃えることが重要である。本例では Pipeline に StandardScaler と LogisticRegression を入れることで、訓練・テスト分割および交差検証の各分割に対して正しい順序で標準化を行う。classification_report で精度、適合率、再現率、F1スコアを算出する。混同行列で予測結果を確認する。k分割交差検証で汎化性能を評価する。出力の「+/-」以降は標準偏差の2倍であり、スコアのばらつきを示す。

演習2．SVM

import pandas as pd
from sklearn.datasets import load_iris
from sklearn.svm import SVC
from sklearn.model_selection import GridSearchCV, train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.metrics import classification_report

# Irisデータセットの読み込み
iris = load_iris()
X = pd.DataFrame(iris.data, columns=iris.feature_names)
y = iris.target

# データの分割（分類問題なのでクラス比率を保つ）
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42, stratify=y
)

# 標準化とSVMをPipelineで一体化
svm_model = Pipeline([
    ('scaler', StandardScaler()),
    ('svc', SVC(kernel='rbf', random_state=42))
])

# パラメータグリッド
param_grid = {
    'svc__C': [0.1, 1, 10],
    'svc__gamma': ['scale', 'auto', 0.1, 1],
}

# グリッドサーチ
grid_search = GridSearchCV(svm_model, param_grid, cv=5, scoring='accuracy', n_jobs=-1)

# モデルの学習
grid_search.fit(X_train, y_train)

# 最適パラメータと評価
best_params = grid_search.best_params_
print(f"Best parameters: C={best_params['svc__C']}, gamma={best_params['svc__gamma']}")
y_pred = grid_search.predict(X_test)
print(classification_report(y_test, y_pred, target_names=iris.target_names))

マージン最大化により分離超平面を求める手法である。RBFカーネル（Radial Basis Function：データ間の類似度をガウス関数で測定するカーネル関数）を使用し、カーネルトリックで非線形分類を実現する。グリッドサーチでハイパーパラメータCとgammaの全組み合わせを探索し、k分割交差検証で最適値を決定する。Cは誤分類に対するペナルティの強さ、gammaはカーネルの影響範囲を制御するパラメータである。SVMは特徴量のスケールに強く影響されるため、標準化を行う。本例では Pipeline を用いるため、交差検証の各分割で標準化が訓練データだけから学習される。

演習3．k-NN

import numpy as np
import pandas as pd
from sklearn.datasets import load_iris
from sklearn.neighbors import KNeighborsClassifier
from sklearn.model_selection import cross_val_score, train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.metrics import classification_report

# Irisデータセットの読み込み
iris = load_iris()
X = pd.DataFrame(iris.data, columns=iris.feature_names)
y = iris.target

# データの分割（分類問題なのでクラス比率を保つ）
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42, stratify=y
)

# k値の範囲でモデル評価
# Pipelineを用いるため、交差検証の各分割内だけで標準化が学習される
k_values = list(range(1, 31))
k_scores = []
for k in k_values:
    knn_model = Pipeline([
        ('scaler', StandardScaler()),
        ('knn', KNeighborsClassifier(n_neighbors=k))
    ])
    scores = cross_val_score(knn_model, X_train, y_train, cv=5)
    k_scores.append(scores.mean())

# 最適なkの選択
best_k = k_values[int(np.argmax(k_scores))]
print(f"Best k: {best_k}")

# 最適モデルの構築と評価
knn = Pipeline([
    ('scaler', StandardScaler()),
    ('knn', KNeighborsClassifier(n_neighbors=best_k))
])
knn.fit(X_train, y_train)
y_pred = knn.predict(X_test)
print(classification_report(y_test, y_pred, target_names=iris.target_names))

訓練データからk個の最近傍点を距離計算で特定し、予測するインスタンスベース学習手法である。k値を1から30の範囲で変化させ、k分割交差検証で最適値を選択する。kが小さいとノイズに敏感になり、kが大きいと決定境界が滑らかになる一方で局所的な構造を捉えにくくなる。k-NNは距離計算を用いるため標準化を行う。本例では Pipeline を用いるため、標準化による検証データから訓練データへの情報混入を避けられる。

演習4．決定木・ランダムフォレスト

import pandas as pd
from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.tree import DecisionTreeClassifier
from sklearn.metrics import classification_report

# Irisデータセットの読み込み
iris = load_iris()
X = pd.DataFrame(iris.data, columns=iris.feature_names)
y = iris.target

# データの分割（分類問題なのでクラス比率を保つ）
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, random_state=42, stratify=y
)

# 決定木による学習と予測
dt = DecisionTreeClassifier(random_state=42)
dt.fit(X_train, y_train)
dt_pred = dt.predict(X_test)
print("Decision Tree:")
print(classification_report(y_test, dt_pred, target_names=iris.target_names))

# ランダムフォレストによる学習と予測
rf = RandomForestClassifier(n_estimators=100, random_state=42)
rf.fit(X_train, y_train)
rf_pred = rf.predict(X_test)
print("\nRandom Forest:")
print(classification_report(y_test, rf_pred, target_names=iris.target_names))

# 特徴量重要度の算出と表示
feature_importance = pd.DataFrame({
    'feature': iris.feature_names,
    'importance': rf.feature_importances_
})
print("\nFeature importance:")
print(feature_importance.sort_values('importance', ascending=False))

決定木はデータを再帰的に分割して判断ルールを生成する手法である。ランダムフォレストは複数の決定木を組み合わせたアンサンブル学習手法であり、複数の決定木の予測を統合することで単一の決定木より過学習を抑制しやすい。両手法でモデルを構築し、classification_report で性能を比較する。特徴量重要度で各特徴量の予測への寄与を評価する。これらの手法は特徴量の大小関係による分岐を用いるため、標準化は行わない。

演習5．k-means

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import silhouette_score

# Irisデータセットの読み込み
iris = load_iris()
X = pd.DataFrame(iris.data, columns=iris.feature_names)

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

# エルボー法とシルエット分析のための値を計算
inertias = []
silhouette_scores = []
k_values = list(range(2, 10))
for k in k_values:
    kmeans = KMeans(n_clusters=k, random_state=42, n_init=10)
    kmeans.fit(X_scaled)
    inertias.append(kmeans.inertia_)
    silhouette_scores.append(silhouette_score(X_scaled, kmeans.labels_))

# エルボー法とシルエット係数のプロット（両者を並べて比較する）
plt.figure(figsize=(12, 5))

# エルボー法（クラスタ内平方和の減少が緩やかになる転換点を読み取る）
plt.subplot(121)
plt.plot(k_values, inertias, 'bo-')
plt.xlabel('Number of clusters (k)')
plt.ylabel('Inertia')
plt.title('Elbow Method')
plt.grid(True)

# シルエット係数（最大となるkが最適クラスタ数の候補）
plt.subplot(122)
plt.plot(k_values, silhouette_scores, 'ro-')
plt.xlabel('Number of clusters (k)')
plt.ylabel('Silhouette score')
plt.title('Silhouette Analysis')
plt.grid(True)

plt.tight_layout()
plt.show()

# クラスタ数ごとのシルエット係数の一覧表示
print("クラスタ数ごとのシルエット係数:")
for k, score in zip(k_values, silhouette_scores):
    print(f"  k={k}: {score:.3f}")

# シルエットスコアが最大となるクラスタ数を選択
best_k = k_values[int(np.argmax(silhouette_scores))]

# 最適なクラスタ数でモデルを再学習
kmeans = KMeans(n_clusters=best_k, random_state=42, n_init=10)
clusters = kmeans.fit_predict(X_scaled)

# 結果の出力
print(f"Best number of clusters: {best_k}")
print(f"Silhouette score: {silhouette_score(X_scaled, clusters):.3f}")

# クラスタリング結果の基本統計量
cluster_stats = X.copy()
cluster_stats['Cluster'] = clusters
print("\nCluster statistics:")
print(cluster_stats.groupby('Cluster').mean())

データ点の割り当てと中心の再計算を収束まで反復する教師なし学習手法である。本例ではエルボー法は可視化による参考とし、最適クラスタ数の決定にはシルエット係数を用いる。KMeans.inertia_ は各点と所属クラスタ中心との距離の二乗和であり、クラスタ内平方和を表す。シルエット係数は同一クラスタ内の凝集度と他クラスタとの分離度を示す。random_state を固定して再現性を確保する。k-meansは距離計算を用いるため標準化を行う。

演習6．階層的クラスタリング

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
from sklearn.cluster import AgglomerativeClustering
from sklearn.preprocessing import StandardScaler
from scipy.cluster.hierarchy import dendrogram, linkage, cophenet
from scipy.spatial.distance import pdist

# Irisデータセットの読み込み
iris = load_iris()
X = pd.DataFrame(iris.data, columns=iris.feature_names)

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

# 階層的クラスタリングの実行（ウォード法）
hierarchical = AgglomerativeClustering(n_clusters=3, linkage='ward')
clusters = hierarchical.fit_predict(X_scaled)

# 結果の表示
print("階層的クラスタリング結果:")
print("各クラスタのサンプル数:", np.bincount(clusters))

# クラスタごとの基本統計量
cluster_stats = X.copy()
cluster_stats['Cluster'] = clusters
print("\nクラスタごとの平均値:")
print(cluster_stats.groupby('Cluster').mean())

# デンドログラムの作成
linkage_matrix = linkage(X_scaled, method='ward', metric='euclidean')
plt.figure(figsize=(10, 7))
dendrogram(linkage_matrix)
plt.title('Hierarchical Clustering Dendrogram')
plt.xlabel('Sample index')
plt.ylabel('Distance')
plt.show()

# コフェネティック相関係数の計算
c, coph_dists = cophenet(linkage_matrix, pdist(X_scaled))
print(f"\nCophenetic correlation coefficient: {c:.3f}")

ボトムアップ（凝集型）で個別データ点から開始して近接クラスタを段階的に統合する手法である。ウォード法は、クラスタを統合したときのクラスタ内平方和の増加が小さくなるように統合を進める方法である。デンドログラムで階層構造を樹形図として可視化する。コフェネティック相関係数で元の距離とデンドログラム上の距離の相関を測定する。値が1に近いほどデンドログラムが元のデータ構造を保持している。階層的クラスタリングは距離に基づくため標準化を行う。

演習7．DBSCAN

import numpy as np
import pandas as pd
from sklearn.datasets import load_iris
from sklearn.cluster import DBSCAN
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import silhouette_score

# Irisデータセットの読み込み
iris = load_iris()
X = pd.DataFrame(iris.data, columns=iris.feature_names)

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

# DBSCANによるクラスタリング
dbscan = DBSCAN(eps=0.5, min_samples=5)
clusters = dbscan.fit_predict(X_scaled)

# クラスタ評価（クラスタラベル -1 はノイズ点を示す）
n_clusters = len(set(clusters)) - (1 if -1 in clusters else 0)
n_noise = list(clusters).count(-1)
print(f"Number of clusters: {n_clusters}")
print(f"Number of noise points: {n_noise}")

# シルエットスコアの計算（ノイズ点を除外。クラスタが2つ以上ある場合のみ）
mask = clusters != -1
if len(set(clusters[mask])) >= 2:
    silhouette_avg = silhouette_score(X_scaled[mask], clusters[mask])
    print(f"Silhouette score: {silhouette_avg:.3f}")
else:
    print("有効なクラスタが2つ未満のため、シルエット係数は計算できません。")

密度ベースのクラスタリング手法である。データ点を、近傍に一定数以上のデータ点を持つコアポイント、境界ポイント、ノイズポイントの3種類に分類する。epsは近傍の半径、min_samplesはコアポイント判定に必要な近傍点数である。クラスタ数の事前指定が不要で、外れ値を自動識別できる。結果はepsとmin_samplesに依存し、値によってはクラスタが1つに集約されたり全点がノイズになるため、値を変えて結果を比較する。シルエット係数は2つ以上のクラスタがあるときに定義されるため、本例ではノイズ点を除いた有効クラスタが2つ以上あるときに計算する。DBSCANは距離計算を用いるため標準化を行う。

演習8．PCA

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

# Irisデータセットの読み込み
iris = load_iris()
X = pd.DataFrame(iris.data, columns=iris.feature_names)
y = iris.target

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

# PCAの実行
pca = PCA()
X_pca = pca.fit_transform(X_scaled)

# 寄与率、累積寄与率、固有値の表示
explained_variance_ratio = pca.explained_variance_ratio_
cumulative_variance_ratio = np.cumsum(explained_variance_ratio)
eigenvalues = pca.explained_variance_

print("Explained variance ratio:", explained_variance_ratio)
print("Cumulative explained variance ratio:", cumulative_variance_ratio)
print("Eigenvalues:", eigenvalues)

# 主成分の方向（固有ベクトル）の表示
print("\n主成分の方向（固有ベクトル）:")
for i, components in enumerate(pca.components_, 1):
    print(f"\n第{i}主成分:")
    for j, component in enumerate(components):
        print(f"  {iris.feature_names[j]}: {component:.4f}")

# 主成分負荷量の表示（元の変数と主成分得点の相関）
pc_names = [f"PC{i}" for i in range(1, X_pca.shape[1] + 1)]
loadings = np.corrcoef(X_scaled.T, X_pca.T)[:X_scaled.shape[1], X_scaled.shape[1]:]
loading_df = pd.DataFrame(loadings, index=iris.feature_names, columns=pc_names)
print("\n主成分負荷量（元の変数と主成分得点の相関）:")
print(loading_df)

# 可視化
plt.figure(figsize=(12, 4))

# スクリープロット（主成分番号に対する寄与率のプロット）
plt.subplot(131)
plt.plot(range(1, len(explained_variance_ratio) + 1), explained_variance_ratio, 'bo-')
plt.xlabel('Principal Component')
plt.ylabel('Explained Variance Ratio')
plt.title('Scree Plot')
plt.grid(True)

# 累積寄与率のプロット
plt.subplot(132)
plt.plot(range(1, len(cumulative_variance_ratio) + 1), cumulative_variance_ratio, 'ro-')
plt.xlabel('Number of Components')
plt.ylabel('Cumulative Explained Variance Ratio')
plt.title('Cumulative Variance')
plt.grid(True)

# 第1・第2主成分でのデータ分布
plt.subplot(133)
scatter = plt.scatter(X_pca[:, 0], X_pca[:, 1], c=y, cmap='viridis')
plt.xlabel('First Principal Component')
plt.ylabel('Second Principal Component')
plt.title('Data Distribution in PC Space')
plt.colorbar(scatter)

plt.tight_layout()
plt.show()

データの分散が最大となる軸を算出して次元削減を行う統計的手法である。寄与率で各主成分が全体の分散を説明する割合を算出する。累積寄与率を基準に保持する主成分数を選択する。固有値は各主成分の分散の大きさを示す。主成分の方向は PCA.components_ で確認でき、主成分負荷量は元の変数と主成分得点の相関として確認できる。スクリープロットで主成分の重要度を可視化する。PCAは分散・共分散の大きさに影響されるため標準化を行う。

演習9．t-SNE

import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
from sklearn.manifold import TSNE
from sklearn.preprocessing import StandardScaler

# Irisデータセットの読み込み
iris = load_iris()
X = pd.DataFrame(iris.data, columns=iris.feature_names)
y = iris.target

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

# t-SNEの実行
tsne = TSNE(
    n_components=2,
    perplexity=30,
    init='pca',
    learning_rate='auto',
    random_state=42
)
X_tsne = tsne.fit_transform(X_scaled)

# 最終的なKLダイバージェンス（高次元と低次元の確率分布の差異）
print(f"KL divergence: {tsne.kl_divergence_:.4f}")

# 結果の可視化
plt.figure(figsize=(10, 6))
scatter = plt.scatter(X_tsne[:, 0], X_tsne[:, 1], c=y, cmap='viridis')
plt.colorbar(scatter)
plt.title('t-SNE visualization of Iris dataset')
plt.xlabel('First t-SNE component')
plt.ylabel('Second t-SNE component')
plt.show()

確率分布の一致度を最適化する非線形次元削減手法である。近傍点間の関係（局所構造）を保持しながら、クラスタ構造の可視化に適した低次元表現を生成する。KL-divergenceで高次元空間と低次元空間の確率分布の一致度を評価し、学習後の kl_divergence_ 属性で最終値を確認できる。perplexityは各データ点の近傍点数の目安を制御するパラメータであり、本例では150サンプルに対して perplexity=30 とする。random_state を固定しないと実行ごとに結果が異なる。t-SNEは距離に基づいて近傍関係を扱うため標準化を行う。

演習10．因子分析

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
from sklearn.decomposition import FactorAnalysis
from sklearn.preprocessing import StandardScaler

# Irisデータセットの読み込み
iris = load_iris()
X = pd.DataFrame(iris.data, columns=iris.feature_names)
y = iris.target

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

# 因子分析の実行（rotation='varimax' で解釈しやすい負荷量を得る）
fa = FactorAnalysis(n_components=2, rotation='varimax', random_state=42)
X_fa = fa.fit_transform(X_scaled)

# 因子負荷量の表示（components_ は形状 (因子数, 変数数)）
print("因子負荷量:")
for i, feature_name in enumerate(iris.feature_names):
    print(f"\n{feature_name}:")
    for j, loading in enumerate(fa.components_[:, i], 1):
        print(f"  因子{j}: {loading:.4f}")

# 共通性の計算（因子方向に負荷量の二乗和をとる）
communalities = np.sum(fa.components_**2, axis=0)
print("\n共通性:")
for feature_name, comm in zip(iris.feature_names, communalities):
    print(f"{feature_name}: {comm:.4f}")

# 因子得点間の相関の計算
factor_corr = np.corrcoef(X_fa.T)
print("\n因子得点間の相関:")
print(factor_corr)

# 可視化
plt.figure(figsize=(12, 4))

# 因子負荷量のプロット
plt.subplot(131)
plt.scatter(fa.components_[0], fa.components_[1])
for i, feature_name in enumerate(iris.feature_names):
    plt.annotate(feature_name, (fa.components_[0, i], fa.components_[1, i]))
plt.xlabel('Factor 1')
plt.ylabel('Factor 2')
plt.title('Factor Loading Plot')
plt.grid(True)
plt.axhline(y=0, color='k', linestyle='--', alpha=0.3)
plt.axvline(x=0, color='k', linestyle='--', alpha=0.3)

# 因子得点の分布
plt.subplot(132)
scatter = plt.scatter(X_fa[:, 0], X_fa[:, 1], c=y, cmap='viridis', alpha=0.6)
plt.xlabel('Factor Score 1')
plt.ylabel('Factor Score 2')
plt.title('Distribution of Factor Scores')
plt.colorbar(scatter)

# 共通性の棒グラフ
plt.subplot(133)
plt.bar(range(len(communalities)), communalities)
plt.xticks(range(len(iris.feature_names)), iris.feature_names, rotation=45, ha='right')
plt.ylabel('Communality')
plt.title('Communalities')
plt.grid(True, axis='y')

plt.tight_layout()
plt.show()

観測変数の背後にある潜在的な共通因子を抽出する統計的手法である。因子負荷量（各観測変数と各潜在因子との関係を表す係数）で各観測変数と潜在因子との関係を評価する。scikit-learnでは components_ 属性（形状は (因子数, 変数数)）から因子負荷量が得られ、行が因子、列が変数に対応する。本例では rotation='varimax' を指定して解釈しやすい負荷量を得る。共通性は各変数について因子方向に負荷量の二乗和をとって算出する。因子分析は変数間の共分散・相関構造を扱うため、標準化を行う。

演習11．線形判別分析（LDA）

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import classification_report

# Irisデータセットの読み込み
iris = load_iris()
X = pd.DataFrame(iris.data, columns=iris.feature_names)
y = iris.target

# データの分割（分類問題なのでクラス比率を保つ）
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, random_state=42, stratify=y
)

# 標準化
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

# LDAの実行
lda = LinearDiscriminantAnalysis()
lda.fit(X_train_scaled, y_train)
X_lda = lda.transform(X_train_scaled)

# 分類用係数の表示（coef_ は形状 (クラス数, 特徴量数)）
print("分類用係数:")
for class_index, coef in enumerate(lda.coef_):
    print(f"\nクラス {class_index} ({iris.target_names[class_index]}):")
    for feature_name, c in zip(iris.feature_names, coef):
        print(f"  {feature_name}: {c:.4f}")

# 分類性能の評価
y_pred = lda.predict(X_test_scaled)
print("\n分類性能:")
print(f"正答率: {lda.score(X_test_scaled, y_test):.4f}")
print("\n分類レポート:")
print(classification_report(y_test, y_pred, target_names=iris.target_names))

# 各判別関数の寄与率
print("\n各判別関数の寄与率:")
explained_ratios = lda.explained_variance_ratio_
for i, ratio in enumerate(explained_ratios, 1):
    print(f"判別関数{i}: {ratio:.4f}")

# 可視化
plt.figure(figsize=(12, 4))

# 第1・第2判別関数でのデータ分布
plt.subplot(131)
scatter = plt.scatter(X_lda[:, 0], X_lda[:, 1], c=y_train, cmap='viridis')
plt.xlabel('Discriminant Function 1')
plt.ylabel('Discriminant Function 2')
plt.title('Data in Discriminant Space')
plt.colorbar(scatter, ticks=[0, 1, 2], label='Class')

# 判別関数の寄与率
plt.subplot(132)
plt.bar(range(1, len(explained_ratios) + 1), explained_ratios)
plt.xlabel('Discriminant Function')
plt.ylabel('Explained Variance Ratio')
plt.title('Explained Ratio of LDA Axes')
plt.grid(True, axis='y')

# 第1判別関数のヒストグラム
plt.subplot(133)
for class_label in np.unique(y_train):
    plt.hist(X_lda[y_train == class_label, 0], alpha=0.5,
             label=iris.target_names[class_label], bins=20)
plt.xlabel('Discriminant Function 1')
plt.ylabel('Frequency')
plt.title('Class Separation by DF1')
plt.legend()
plt.grid(True, axis='y')

plt.tight_layout()
plt.show()

分類と次元削減の両機能を持つ教師あり学習手法である。クラス間分散を最大化し、クラス内分散を最小化する軸を求める。coef_（形状 (クラス数, 特徴量数)）は各クラスを分離する分類用の重みであり、kクラスで最大k-1本得られる判別軸そのものとは異なる。判別軸に射影された値は transform() が返す判別空間に対応する。分類性能を評価し、判別空間でのデータ分布を可視化する。ヒストグラムでクラスが分離されているほど判別関数が有効である。LDAは分散・共分散構造を扱い、係数の比較も特徴量スケールの影響を受けるため標準化を行う。

演習12．特異値分解（SVD）とQR分解

import numpy as np
import matplotlib.pyplot as plt

# サンプルデータの生成
np.random.seed(42)
n_rows, n_cols = 100, 50
rank_true = 3

# 直交基底を作り、指定した特異値を持つランク3行列を生成する
# ノイズを加える前の行列 X は、特異値 10, 5, 2 を持つ
U_random = np.random.normal(0, 1, (n_rows, rank_true))
V_random = np.random.normal(0, 1, (n_cols, rank_true))
U_true, _ = np.linalg.qr(U_random)
V_true, _ = np.linalg.qr(V_random)

singular_values_true = np.array([10.0, 5.0, 2.0])
X = U_true @ np.diag(singular_values_true) @ V_true.T

# ノイズの追加
noise_scale = 0.05
X_noisy = X + np.random.normal(0, noise_scale, (n_rows, n_cols))

# SVDの実行
U, s, Vt = np.linalg.svd(X_noisy, full_matrices=False)

# 特異値の表示
print("指定した非ゼロ特異値:")
print(singular_values_true)

print("\nSVDで得られた特異値（上位10個）:")
for i, value in enumerate(s[:10], 1):
    print(f"σ_{i}: {value:.4f}")

# 累積寄与率の計算（特異値の二乗は分散に比例するため、二乗の比で算出する）
cumulative_variance = np.cumsum(s**2) / np.sum(s**2)
print(f"\n上位3成分の累積寄与率: {cumulative_variance[rank_true - 1]:.4f}")

# 低ランク近似の誤差を計算
errors = []
rank_values = list(range(1, min(n_rows, n_cols) + 1))
for r in rank_values:
    X_approx = U[:, :r] @ np.diag(s[:r]) @ Vt[:r, :]
    errors.append(np.linalg.norm(X_noisy - X_approx, 'fro'))

# 可視化用の色範囲を固定
subset = X_noisy[:20, :20]
vmin, vmax = subset.min(), subset.max()

# 可視化
plt.figure(figsize=(15, 4))

# 特異値のプロット
plt.subplot(141)
plt.plot(range(1, len(s) + 1), s, 'o-')
plt.yscale('log')
plt.xlabel('Component number')
plt.ylabel('Singular value')
plt.title('Singular Values')
plt.grid(True)

# 累積寄与率
plt.subplot(142)
plt.plot(range(1, len(s) + 1), cumulative_variance, 'o-')
plt.xlabel('Number of components')
plt.ylabel('Cumulative ratio')
plt.title('Cumulative Variance Ratio')
plt.grid(True)

# 近似誤差
plt.subplot(143)
plt.plot(rank_values[:20], errors[:20], 'o-')
plt.xlabel('Rank used')
plt.ylabel('Approximation error')
plt.title('Low-rank Approximation Error')
plt.grid(True)

# オリジナルデータの可視化
plt.subplot(144)
plt.imshow(subset, cmap='coolwarm', aspect='auto', vmin=vmin, vmax=vmax)
plt.title('Original data subset')
plt.colorbar()

plt.tight_layout()
plt.show()

# 低ランク近似の比較
plt.figure(figsize=(15, 4))
ranks_demo = [1, 3, 10]
for i, r in enumerate(ranks_demo):
    X_approx = U[:, :r] @ np.diag(s[:r]) @ Vt[:r, :]
    plt.subplot(1, len(ranks_demo), i + 1)
    plt.imshow(X_approx[:20, :20], cmap='coolwarm', aspect='auto', vmin=vmin, vmax=vmax)
    plt.title(f'Rank {r} approximation')
    plt.colorbar()

plt.tight_layout()
plt.show()

# QR分解
np.random.seed(42)
n_rows, n_cols = 5, 3
A = np.random.randn(n_rows, n_cols)

# QR分解の実行（reduced QR：Qの列が互いに直交する）
Q, R = np.linalg.qr(A, mode='reduced')

# 結果の表示
print("元の行列 A:")
print(A)
print("\nQ行列（列直交行列）:")
print(Q)
print("\nR行列（上三角行列）:")
print(R)

# 検証
print("\n検証結果:")
print("Q^T Q（単位行列に近いことを確認）:")
print(np.round(Q.T @ Q, decimals=10))
print("\nQR（元の行列Aに近いことを確認）:")
reconstruction_error = np.linalg.norm(A - Q @ R)
print(f"再構成誤差: {reconstruction_error:.10f}")

# QR分解の応用：最小二乗法による線形回帰
print("\n\nQR分解の応用：最小二乗法による線形回帰")

# サンプルデータの生成
n_samples = 100
X_data = np.random.uniform(0, 10, n_samples)
y_data = 2 * X_data + 1 + np.random.normal(0, 1, n_samples)

# 設計行列の作成（各列が回帰式の説明変数と定数項に対応する行列）
A_reg = np.vstack([X_data, np.ones(n_samples)]).T

# QR分解を使用した最小二乗法
Q_reg, R_reg = np.linalg.qr(A_reg, mode='reduced')
beta = np.linalg.solve(R_reg, Q_reg.T @ y_data)

print(f"\n推定された傾き: {beta[0]:.4f}")
print(f"推定された切片: {beta[1]:.4f}")

# 結果の可視化
plt.figure(figsize=(8, 6))
plt.scatter(X_data, y_data, alpha=0.5, label='Data')
X_sort = np.sort(X_data)
plt.plot(X_sort, beta[0] * X_sort + beta[1], 'r-', linewidth=2, label='Regression line by QR')
plt.plot(X_sort, 2 * X_sort + 1, 'g--', linewidth=2, label='True relation')
plt.xlabel('X')
plt.ylabel('y')
plt.title('Linear Regression using QR Decomposition')
plt.legend()
plt.grid(True)
plt.show()

特異値分解（SVD）は行列を3つの行列（U、Σ、V^T）の積に分解する手法である。特異値で各成分の重要度を評価する。本例では、直交基底を用いて、ノイズを加える前の行列が特異値10、5、2を持つランク3行列になるように生成する。その後に小さいノイズを加えるため、観測行列 X_noisy は厳密にはフルランクになるが、上位3つの特異値が支配的になる。特異値の二乗は各成分が説明する分散に比例するため、累積寄与率は特異値の二乗の比として算出する。低ランク近似では小さい特異値を無視してデータを圧縮する。近似誤差はフロベニウスノルム（行列の全要素の二乗和の平方根）で測定する。

QR分解は行列をQとRの積に分解する手法である。本例のように縦長の長方行列に対して mode='reduced' を用いる場合、Qは正方の直交行列ではなく、列が互いに直交する列直交行列である。このとき Q.T @ Q は単位行列に近くなり、Q @ R は元の行列Aを再構成する。応用として最小二乗法による線形回帰を実装し、QR分解で回帰係数を求める。QR分解を用いることで、正規方程式を直接解く方法より数値的に安定した計算ができる（正規方程式はA^TAの条件数が大きい場合に数値誤差が拡大しやすいため、QR分解による方法が用いられる）。