機械学習のPython実現ガイド

機械学習の用語

基本概念

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

教師あり学習手法

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

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

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

次元削減・特徴抽出手法

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

行列分解手法

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

評価指標

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

実装のための技術

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

実装結果の検証・評価

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

Pythonでの実現

演習

以下の演習では、Irisデータセット（4特徴量、3クラス、150サンプルから構成される分類用の標準データセット）を使用する。距離や内積を用いるアルゴリズム（SVM、k-NN、PCA、t-SNE、DBSCAN、因子分析、LDA）では標準化を行う。決定木およびランダムフォレストはスケールに依存しないため標準化を行わない。各演習で共通する処理（データ読み込み、データ分割、標準化、classification_report による評価、random_state=42 の指定）は同一の記述に揃えている。

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

from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split, cross_val_score
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)

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

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

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

# 5分割交差検証による性能評価
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})")

線形モデルとシグモイド関数を組み合わせた分類手法である。最尤推定により重みパラメータを最適化し、テストデータで予測する。classification_report で精度、適合率、再現率、F1スコアを算出する。混同行列で予測結果を真陽性、偽陽性、真陰性、偽陰性に分類して表示する。k分割交差検証で汎化性能を評価する。出力の「+/-」以降は標準偏差の2倍であり、スコアのばらつきを示す。scikit-learn 1.8以降、LogisticRegression の penalty および n_jobs 引数は非推奨であり将来削除される予定である。本例はこれらを指定していないため警告なく動作する。正則化の種類を指定する場合は、penalty ではなく l1_ratio（L2なら0.0、L1なら1.0）を用いる。

演習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.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)

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

# モデルの構築
svm = SVC(kernel='rbf', random_state=42)

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

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

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

# 最適パラメータと評価
print(f"Best parameters: {grid_search.best_params_}")
y_pred = grid_search.predict(X_test_scaled)
print(classification_report(y_test, y_pred))

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

演習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.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)

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

# k値の範囲でモデル評価
k_range = range(1, 31)
k_scores = []
for k in k_range:
    knn = KNeighborsClassifier(n_neighbors=k)
    scores = cross_val_score(knn, X_train_scaled, y_train, cv=5)
    k_scores.append(scores.mean())

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

# 最適モデルの構築と評価
knn = KNeighborsClassifier(n_neighbors=best_k)
knn.fit(X_train_scaled, y_train)
y_pred = knn.predict(X_test_scaled)
print(classification_report(y_test, y_pred))

訓練データからk個の最近傍点を距離計算で特定し、予測するインスタンスベース学習手法である。k値を1から30の範囲で変化させ、k分割交差検証で最適値を選択する。kが小さいとノイズに敏感になり、kが大きいと決定境界が不明瞭になる。距離計算を用いるため標準化を行う。

演習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)

# 決定木による学習と予測
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))

# ランダムフォレストによる学習と予測
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))

# 特徴量重要度の算出と表示
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)
y = iris.target

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

# エルボー法とシルエット分析のための値を計算
inertias = []
silhouette_scores = []
K = range(2, 10)
for k in K:
    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=(8, 6))
plt.plot(K, inertias, 'bo-')
plt.xlabel('Number of clusters (k)')
plt.ylabel('Inertia')
plt.title('Elbow Method')
plt.grid(True)
plt.show()

# シルエットスコアが最大となるクラスタ数を選択
best_k = K[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 = pd.DataFrame(X)
cluster_stats['Cluster'] = clusters
print("\nCluster statistics:")
print(cluster_stats.groupby('Cluster').mean())

データ点の割り当てと中心の再計算を収束まで反復する教師なし学習手法である。本例ではエルボー法は可視化による参考とし、最適クラスタ数の決定にはシルエット係数を用いる。エルボー法はクラスタ内分散の減少が緩やかになる転換点を、シルエット係数は同一クラスタ内の凝集度と他クラスタとの分離度を示す。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)
y = iris.target

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

# 階層的クラスタリングの実行
hierarchical = AgglomerativeClustering(n_clusters=3)
clusters = hierarchical.fit_predict(X_scaled)

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

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

# デンドログラムの作成
linkage_matrix = linkage(X_scaled, method='ward')
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)
y = iris.target

# 標準化
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つに集約されたり全点がノイズになることがある。シルエット係数の計算ではノイズ点を除外する。距離計算を用いるため標準化を行う。

演習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}")

# 可視化
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()

データの分散が最大となる軸を算出して次元削減を行う統計的手法である。寄与率で各主成分が全体の分散を説明する割合を算出する。累積寄与率を基準に保持する主成分数を選択する。スクリープロットで主成分の重要度を可視化する。距離計算を用いるため標準化を行う。

演習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, 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は各データ点の近傍点数の目安を制御するパラメータである（t-SNE原論文では5から50の範囲が推奨されている）。random_state を固定しないと実行ごとに結果が異なる。距離計算を用いるため標準化を行う。

演習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('第1因子')
plt.ylabel('第2因子')
plt.title('因子負荷量プロット')
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('第1因子得点')
plt.ylabel('第2因子得点')
plt.title('因子得点の分布')
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('共通性')
plt.title('各変数の共通性')
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)

# 標準化
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 i, coef in enumerate(lda.coef_, 1):
    print(f"\nクラス{i}:")
    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('第1判別関数')
plt.ylabel('第2判別関数')
plt.title('判別空間でのデータ分布')
plt.colorbar(scatter, ticks=[0, 1, 2], label='Class')

# 判別関数の寄与率
plt.subplot(132)
plt.bar(range(1, len(explained_ratios) + 1), explained_ratios)
plt.xlabel('判別関数')
plt.ylabel('寄与率')
plt.title('判別関数の寄与率')
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('第1判別関数')
plt.ylabel('頻度')
plt.title('第1判別関数によるクラス分離')
plt.legend()
plt.grid(True, axis='y')

plt.tight_layout()
plt.show()

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

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

import numpy as np
import matplotlib.pyplot as plt

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

# ランク3の行列を生成（ランクとは行列の線形独立な行または列の最大数）
U_true = np.random.normal(0, 1, (n_rows, 3))
V_true = np.random.normal(0, 1, (n_cols, 3))
singular_values_true = np.array([10, 5, 2])
X = U_true @ np.diag(singular_values_true) @ V_true.T

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

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

# 特異値の表示
print("特異値（上位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[2]:.4f}")

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

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

# 特異値のプロット
plt.subplot(141)
plt.plot(range(1, len(s) + 1), s, 'o-')
plt.yscale('log')
plt.xlabel('成分番号')
plt.ylabel('特異値')
plt.title('特異値のスクリープロット')
plt.grid(True)

# 累積寄与率
plt.subplot(142)
plt.plot(range(1, len(s) + 1), cumulative_variance, 'o-')
plt.xlabel('成分数')
plt.ylabel('累積寄与率')
plt.title('累積寄与率')
plt.grid(True)

# 近似誤差
plt.subplot(143)
plt.plot(ranks[:20], errors[:20], 'o-')
plt.xlabel('使用する特異値の数')
plt.ylabel('近似誤差')
plt.title('低ランク近似の誤差')
plt.grid(True)

# オリジナルデータの可視化
plt.subplot(144)
plt.imshow(X_noisy[:20, :20], cmap='coolwarm', aspect='auto')
plt.title('オリジナルデータ（部分）')
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')
    plt.title(f'ランク{r}での近似')
    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分解の実行
Q, R = np.linalg.qr(A)

# 結果の表示
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)
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='データ点')
X_sort = np.sort(X_data)
plt.plot(X_sort, beta[0] * X_sort + beta[1], 'r-', linewidth=2, label='QR分解による回帰直線')
plt.plot(X_sort, 2 * X_sort + 1, 'g--', linewidth=2, label='真の関係')
plt.xlabel('X')
plt.ylabel('y')
plt.title('QR分解を用いた線形回帰')
plt.legend()
plt.grid(True)
plt.show()

特異値分解（SVD）は行列を3つの行列（U、Σ、V^T）の積に分解する手法である。特異値で各成分の重要度を評価する。特異値の二乗は各成分が説明する分散に比例するため、累積寄与率は特異値の二乗の比として算出する（PCAの固有値に対応する）。低ランク近似で小さい特異値を無視してデータを圧縮する。近似誤差はフロベニウスノルム（行列の全要素の二乗和の平方根）で測定する。QR分解は行列を直交行列Qと上三角行列Rの積に分解する手法である。直交行列QはQ^TQ = Iを満たす。Q^TQが単位行列に近いこと、QRの積が元の行列を再構成できることを検証する。応用として最小二乗法による線形回帰を実装し、QR分解で回帰係数を求める。QR分解を用いることで、正規方程式を直接解く方法より数値的に安定した計算ができる（正規方程式はA^TAの条件数が大きい場合に数値誤差が拡大しやすいため、QR分解による方法が用いられる）。