ラビット・チャレンジ - Stage 2. 機械学習

提出したレポートです。

絶対書きすぎですが、行間を埋めたくなるので仕方ない。

---

Rabbit Challenge - Stage 2. 機械学習

機械学習の課題

機械学習

機械学習 とは、データからパターンや規則性を学習し、その学習結果を用いて新しいデータに対して予測や分類を行う技術である。

訓練誤差

訓練誤差 とは、機械学習モデルが訓練データに対してどれだけ正確に予測を行ったかを示す指標である。訓練誤差が小さいほど、モデルは訓練データに対して良い適合を示す。

モデル $f(\boldsymbol{x}; \theta)$ が与えられて、損失関数 $L(y, f(\boldsymbol{x}; \theta))$ を用いるとき、訓練誤差は以下のように定義される。

$$ E_{train}(\theta) = \frac{1}{N} \sum_{i=1}^{N} L(y_i, f(\boldsymbol{x}_i; \theta)) $$

ここで、$N$ は訓練データのサンプル数、$y_i$ はサンプル $i$ の真のラベル、$\boldsymbol{x}_i$ はサンプル $i$ の特徴ベクトルである。

汎化誤差

汎化誤差 とは、機械学習モデルが未知のデータに対してどれだけ正確に予測を行えるかを示す指標である。汎化誤差が小さいほど、モデルは新しいデータに対して良い適合を示す。

モデル $f(\boldsymbol{x}; \theta)$ が与えられて、損失関数 $L(y, f(\boldsymbol{x}; \theta))$ を用いるとき、汎化誤差は以下のように定義される。

$$ E_{gen}(\theta) = \mathbb{E}_{(\boldsymbol{x}, y) \sim P_{data}} [L(y, f(\boldsymbol{x}; \theta))] $$

ここで、$P_{data}$ はデータの真の分布を表し、$\mathbb{E}$ は期待値を示す。

過剰適合・過少適合

過剰適合 とは、機械学習モデルが訓練データに対して非常に良い適合を示す一方で、未知のデータに対しては悪い適合を示す現象である。

これは、訓練誤差が非常に小さい一方で、汎化誤差が大きい場合に発生する。

原因としては、以下が挙げられる。

• モデルが過度に複雑であり、訓練データのノイズや偶然のパターンまで学習してしまう
• 訓練データが少ない場合、モデルは限られたデータに過度に適応してしまう

過少適合 とは、機械学習モデルが訓練データに対しても未知のデータに対しても良い適合を示さない現象である。

原因としては、以下が挙げられる。

• モデルがあまりにも単純であり、データの複雑なパターンを捉えきれない
• 入力データに重要な特徴量が欠けている場合、モデルは十分な情報を得られず、適切な学習ができない

正則化

正則化 とは、機械学習モデルの複雑さを制御し、過剰適合を防ぐための手法である。正則化は、損失関数にペナルティ項を追加することで実現される。

訓練誤差 $E_{train}(\theta)$ に正則化項 $\Omega(\theta)$ を追加した正則化損失関数は以下のように定義される。

$$ E_{reg}(\theta) = E_{train}(\theta) + \lambda \Omega(\theta) $$

ここで、$\lambda$ は正則化の強さを制御するハイパーパラメータであり、$\Omega(\theta)$ はモデルの複雑さを測る関数である。

• L2 正則化（リッジ回帰）：大きなパラメータ値を抑制して、モデルを滑らかにする。

$$ \Omega(\theta) = \sum_{j} \theta_j^2 $$

• L1 正則化（ラッソ回帰）：不要な特徴量のパラメータをゼロにすることで、特徴選択の効果を持つ。

$$ \Omega(\theta) = \sum_{j} |\theta_j| $$

次元の呪い

次元の呪い とは、データの次元数が増加するにつれて、機械学習モデルの性能が低下する現象を指す。

高次元空間では、データポイントが疎になり、距離や密度の概念が直感的でなくなるため、モデルの学習が困難になる。

深層学習

深層学習 とは、多層のニューラルネットワークを用いてデータから特徴を自動的に学習し、複雑なパターンを捉える機械学習の一分野である。

深層学習の適用範囲は、以下のように多岐にわたる。

• 画像認識
• 音声認識
• 自然言語処理

---

性能指標

混同行列

混同行列は、分類モデルの性能を評価するための表であり、実際のクラスと予測されたクラスの関係を示す。

	予測: ポジティブ	予測: ネガティブ
実際: ポジティブ	True Positive (TP)	False Negative (FN)
実際: ネガティブ	False Positive (FP)	True Negative (TN)

正解率

正解率 (Accuracy) とは、分類モデルが正しく予測したサンプルの割合を示す指標である。

$$ \text{Accuracy} = \frac{TP + TN}{TP + TN + FP + FN} $$

クラスの分布が偏っている場合、正解率は誤解を招く可能性があるため、他の指標と併用することが推奨される。

適合率

適合率 (Precision) とは、モデルがポジティブと予測したサンプルのうち、実際にポジティブであった割合を示す指標である。

$$ \text{Precision} = \frac{TP}{TP + FP} $$

以下のように、誤ってポジティブと予測することを避けたい場合に重要である。

• スパムメール検出
• 医療診断
• 異常検知システム

再現率

再現率 (Recall) とは、実際にポジティブであるサンプルのうち、モデルが正しくポジティブと予測した割合を示す指標である。

$$ \text{Recall} = \frac{TP}{TP + FN} $$

以下のように、ポジティブサンプルを見逃すことを避けたい場合に重要である。

• 疾病スクリーニング
• クレジットカード不正検出
• セキュリティ侵害検出

F 値

F 値 (F1 Score) とは、適合率と再現率の調和平均を示す指標であり、両者のバランスを評価するために用いられる。

$$ \text{F1 Score} = 2 \times \frac{\text{Precision} \times \text{Recall}}{\text{Precision} + \text{Recall}} $$

クラスの分布が偏っている場合や、適合率と再現率の両方を考慮したい場合に有用である。

ROC 曲線と AUC

ROC 曲線 (Receiver Operating Characteristic Curve) とは、分類モデルの性能を評価するためのグラフであり、偽陽性率 (FPR) に対する真陽性率 (TPR) をプロットしたものである。

• 偽陽性率 (FPR)

$$ \text{FPR} = \frac{FP}{FP + TN} $$

• 真陽性率 (TPR)

$$ \text{TPR} = \frac{TP}{TP + FN} $$

AUC (Area Under the Curve) とは、ROC 曲線の下の面積を示す指標であり、分類モデルの全体的な性能を評価するために用いられる。

AUC の値は 0 から 1 の範囲であり、1 に近いほどモデルの性能が高いことを示す。

micro 平均と macro 平均

micro 平均 とは、マルチクラス分類において、全てのクラスの TP、FP、FN を合計してから適合率、再現率、F 値を計算する方法である。

$$ \text{Precision}_{micro} = \frac{\sum_{c} TP_c}{\sum_{c} (TP_c + FP_c)} \\ \text{Recall}_{micro} = \frac{\sum_{c} TP_c}{\sum_{c} (TP_c + FN_c)} $$

クラスの不均衡がある場合に、全体的な性能を評価するのに適している。

macro 平均 とは、マルチクラス分類において、各クラスの適合率、再現率、F 値を個別に計算し、その平均を取る方法である。

$$ \text{Precision}_{macro} = \frac{1}{C} \sum_{c} \frac{TP_c}{TP_c + FP_c} \\ \text{Recall}_{macro} = \frac{1}{C} \sum_{c} \frac{TP_c}{TP_c + FN_c} $$

クラスごとの性能を均等に評価するのに適している。

MSE・RMSE・MAE・決定係数

平均二乗誤差 (MSE: Mean Squared Error) とは、回帰モデルの予測値と実際の値との誤差の二乗の平均を示す指標である。

$$ \text{MSE} = \frac{1}{N} \sum_{i=1}^{N} (y_i - \hat{y}_i)^2 $$

平方平均二乗誤差 (RMSE: Root Mean Squared Error) とは、MSE の平方根を取ったものであり、元のデータの単位に戻した指標である。

$$ \text{RMSE} = \sqrt{\text{MSE}} = \sqrt{\frac{1}{N} \sum_{i=1}^{N} (y_i - \hat{y}_i)^2} $$

平均絶対誤差 (MAE: Mean Absolute Error) とは、回帰モデルの予測値と実際の値との誤差の絶対値の平均を示す指標である。

$$ \text{MAE} = \frac{1}{N} \sum_{i=1}^{N} |y_i - \hat{y}_i| $$

決定係数 (R²: Coefficient of Determination) とは、回帰モデルがデータの分散をどれだけ説明できるかを示す指標である。

$$ R^2 = 1 - \frac{\sum_{i=1}^{N} (y_i - \hat{y}_i)^2}{\sum_{i=1}^{N} (y_i - \bar{y})^2} $$

---

線形回帰モデル

線形回帰

線形回帰 とは、教師あり学習の一種であり、説明変数と目的変数の関係を線形関数でモデル化する手法である。

いま、

• 説明変数 $\boldsymbol{x} = (x_1, x_2, \ldots, x_m)^T$
• 目的変数 $y$
• 重み $\boldsymbol{w} = (w_1, w_2, \ldots, w_m)^T$
• バイアス $b$

とすると、線形回帰モデルは以下のように表される。

$$ \hat{y} = f(\boldsymbol{x}; \boldsymbol{w}, b) = w_1 x_1 + w_2 x_2 + \ldots + w_m x_m + b = \boldsymbol{w}^T \boldsymbol{x} + b $$

重み $\boldsymbol{w}$ とバイアス $b$ は、訓練データに基づいて最適化される。一般的には、最小二乗法を用いてパラメータを推定する。

例：線形単回帰 $y = 3x + 2$

例として、線形単回帰 $y = 3x + 2$ を予測する。

サンプル数を $N$ とする。

また、$i$ 番目のサンプルの説明変数と目的変数をそれぞれ $x_i$、$y_i$ とする。

このとき、真のデータの生成過程は以下のようになる。

$$ y_i = 3 x_i + 2 + \epsilon_i $$

ここで、$\epsilon_i$ は平均 0、分散 $\sigma^2$ の正規分布に従うノイズである。

import numpy as np

N = 100
"""サンプル数"""

w_true = 3
"""真の重みの値"""

b_true = 2
"""真のバイアスの値"""

sigma = 0.5
"""ノイズの標準偏差"""

def generate_true_data(x):
    """真のデータを生成する"""
    # 平均 0、標準偏差 sigma の正規分布に従うノイズを生成する
    epsilon = np.random.normal(0, sigma)
    return w_true * x + b_true + epsilon

この訓練データに対して、予測モデルを以下のように定義する。

$$ \hat{y}_i = w x_i + b $$

重み $w$ とバイアス $b$ を求める。

まず、損失関数には、以下のように二乗誤差損失を用いる。

$$ L(y_i, \hat{y}_i) = (y_i - \hat{y}_i)^2 = \{ y_i - (w x_i + b) \}^2 $$

経験損失は、損失関数の平均として、以下のように定義される。

$$ \hat{R}_N (w, b) = \frac{1}{N} \sum_{i=1}^{N} L(y_i, \hat{y}_i) = \frac{1}{N} \sum_{i=1}^{N} \{ y_i - (w x_i + b) \}^2 $$

この経験損失が最小となるような、重み $\hat{w}$ とバイアス $\hat{b}$ を求める。これが「学習」に相当する。

$$ (\hat{w}, \hat{b}) = \arg\min_{w, b} \hat{R}_N (w, b) $$

経験損失 $\hat{R}_N (w, b)$ を $w, b$ でそれぞれ偏微分して、0 となる値が最小値を与える $\hat{w}, \hat{b}$ である。

$$ \begin{align*} \frac{\partial \hat{R}_N (w, b)}{\partial w} &= \frac{\partial}{\partial w} \left( \frac{1}{N} \sum_{i=1}^{N} \{ y_i - (w x_i + b) \}^2 \right) \\ &= \frac{\partial}{\partial w} \left( \frac{1}{N} \sum_{i=1}^{N} ( y_i^2 + w^2 x_i^2 + b^2 - 2 w x_i y_i + 2 b w x_i - 2 b y_i ) \right) \\ &= \frac{1}{N} \sum_{i=1}^{N} ( 2 w x_i^2 - 2 x_i y_i + 2 b x_i ) \\ &= -\frac{2}{N} \sum_{i=1}^{N} x_i \{ y_i - (w x_i + b) \} \\ \frac{\partial \hat{R}_N (w, b)}{\partial b} &= \frac{\partial}{\partial b} \left( \frac{1}{N} \sum_{i=1}^{N} \{ y_i - (w x_i + b) \}^2 \right) \\ &= \frac{\partial}{\partial b} \left( \frac{1}{N} \sum_{i=1}^{N} ( y_i^2 + w^2 x_i^2 + b^2 - 2 w x_i y_i + 2 b w x_i - 2 b y_i ) \right) \\ &= \frac{1}{N} \sum_{i=1}^{N} ( 2 b + 2 w x_i - 2 y_i ) \\ &= -\frac{2}{N} \sum_{i=1}^{N} \{ y_i - (w x_i + b) \} \end{align*} $$

これらをそれぞれ $0$ とおくと、以下の連立方程式が得られる。

$$ \begin{cases} \sum_{i=1}^{N} x_i \{ y_i - (\hat{w} x_i + \hat{b}) \} = 0 \\ \sum_{i=1}^{N} \{ y_i - (\hat{w} x_i + \hat{b}) \} = 0 \end{cases} $$

これを解くと、以下のように重み $\hat{w}$ とバイアス $\hat{b}$ が求まる。

$$ \hat{w} = \frac{N \sum_{i=1}^{N} x_i y_i - \sum_{i=1}^{N} x_i \sum_{i=1}^{N} y_i}{N \sum_{i=1}^{N} x_i^2 - (\sum_{i=1}^{N} x_i)^2} $$

$$ \hat{b} = \frac{\sum_{i=1}^{N} y_i - \hat{w} \sum_{i=1}^{N} x_i}{N} $$

import numpy as np

def find_w_hat(x_array, y_array):
    """予測の重み w_hat を求める"""
    sum_x = np.sum(x_array)
    sum_y = np.sum(y_array)
    sum_x2 = np.sum(x_array**2)
    sum_xy = np.sum(x_array * y_array)
    return (N * sum_xy - sum_x * sum_y) / (N * sum_x2 - sum_x**2)

def find_b_hat(x_array, y_array, w_hat):
    """予測のバイアス b_hat を求める"""
    sum_x = np.sum(x_array)
    sum_y = np.sum(y_array)
    return (sum_y - w_hat * sum_x) / N

では、実際にこの予測モデルを使ってみる。

import numpy as np
import matplotlib.pyplot as plt

# 訓練データを生成する
x_train = np.random.uniform(0, 10, N)
y_train = np.array([generate_true_data(x) for x in x_train])

# 予測モデルのパラメータを求める
w_hat = find_w_hat(x_train, y_train)
b_hat = find_b_hat(x_train, y_train, w_hat)
print(f"Estimated w: {w_hat}, Estimated b: {b_hat}")

# 予測モデルを用いて予測を行う
y_pred = w_hat * x_train + b_hat

# 結果をプロットする
plt.scatter(x_train, y_train, label="Training Data", color="blue")
plt.plot(x_train, y_pred, label="Predicted Line", color="red")
plt.xlabel("x")
plt.ylabel("y")
plt.title("Linear Regression Prediction")
plt.legend()
plt.show()

例：線形重回帰 $y = 2x_1 + 3x_2 + 5$

例として、線形重回帰 $y = 2x_1 + 3x_2 + 5$ を予測する。

前の例では、学習部分を手計算で行ったが、scikit-learn を用いると簡単に実装できる。

import numpy as np
from sklearn.linear_model import LinearRegression

N = 100
"""サンプル数"""

w_true = np.array([2, 3])
"""真の重みの値"""

b_true = 5
"""真のバイアスの値"""

sigma = 0.5
"""ノイズの標準偏差"""

def generate_true_data(x1, x2):
    """真のデータを生成する"""
    # 平均 0、標準偏差 sigma の正規分布に従うノイズを生成する
    epsilon = np.random.normal(0, sigma)
    return w_true[0] * x1 + w_true[1] * x2 + b_true + epsilon

# 訓練データを生成する
x1_train = np.random.uniform(0, 10, N)
x2_train = np.random.uniform(0, 10, N)
y_train = np.array([generate_true_data(x1, x2) for x1, x2 in zip(x1_train, x2_train)])
X_train = np.column_stack((x1_train, x2_train))

# 予測モデルを作成し、学習を行う
model = LinearRegression()
model.fit(X_train, y_train)

# 学習したパラメータを表示する
print(f"Estimated w: {model.coef_}, Estimated b: {model.intercept_}")

# 予測モデルを用いて予測を行う
y_pred = model.predict(X_train)

# 結果を表示する
for i in range(5):
    print(f"True y: {y_train[i]}, Predicted y: {y_pred[i]}")

---

非線形回帰モデル

非線形回帰 とは、説明変数と目的変数の関係を非線形関数でモデル化する手法である。非線形回帰モデルは、線形回帰モデルでは捉えきれない複雑なパターンを学習することができる。

本稿では特に「入力に対して非線形だが、パラメータに対して線形なモデル」を非線形回帰として扱う。

基底関数展開

基底関数展開 とは、非線形関数を複数の基底関数の線形結合として表現する手法である。

これにより、非線形回帰モデルを線形回帰モデルの枠組みで扱うことができる。

線形回帰の行列表現

非線形回帰モデルを基底関数展開すれば、線形回帰のように扱えると言われても、とても同じようには見えない。

そのため、既にやった線形回帰を行列を用いて表現して一般化することで、基底関数展開を用いた非線形回帰モデルも同様に扱えるという理解につなげる。

まず、線形回帰モデルを行列表現で表す。

訓練データが $N$ 個あり、各サンプルの説明変数が $m$ 次元であるとする。

各サンプルの説明変数ベクトルを

$$ \boldsymbol{x}_i = (x_{i1}, x_{i2}, \ldots, x_{im})^T \quad (i = 1, \ldots, N) $$

とおく。

バイアス項を含めるため、各サンプルに定数 $1$ を付加した特徴ベクトルを

$$ \tilde{\boldsymbol{x}}_i = (1, x_{i1}, x_{i2}, \ldots, x_{im})^T $$

と定義する。

このとき、訓練データ全体の 説明変数行列（設計行列、計画行列） $\mathbf{X}$ は、

$$ \mathbf{X} = \begin{bmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1m} \\ 1 & x_{21} & x_{22} & \cdots & x_{2m} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{N1} & x_{N2} & \cdots & x_{Nm} \end{bmatrix} $$

で与えられる。

この説明変数行列 $\mathbf{X}$ を用いることで、線形回帰モデルは簡潔に表現できる。

各サンプルに対する目的変数を $y_i \ (i = 1, 2, \ldots, N)$ とし、それらをまとめた目的変数ベクトルを

$$ \boldsymbol{y} = (y_1, y_2, \ldots, y_N)^T $$

と定義する。

また、バイアス項を含めたパラメータベクトルを

$$ \boldsymbol{w} = (b, w_1, w_2, \ldots, w_m)^T $$

とおく。

このとき、線形回帰モデルは各サンプルごとの式

$$ y_i = b + \sum_{j=1}^{m} w_j x_{ij} + \varepsilon_i $$

をまとめて、行列を用いて次のように表せる。

$$ \boldsymbol{y} = \mathbf{X}\boldsymbol{w} + \boldsymbol{\varepsilon} $$

ここで $\boldsymbol{\varepsilon} = (\varepsilon_1, \ldots, \varepsilon_N)^T$ はノイズベクトルである。

次に、最尤推定との関係を確認する。

各ノイズが平均 $0$、分散 $\sigma^2$ の正規分布に従うと仮定すると、

$$ \varepsilon_i \sim \mathcal{N}(0, \sigma^2) $$

が成り立つ。このとき、目的変数ベクトル $\boldsymbol{y}$ の条件付き分布は

$$ p(\boldsymbol{y} \mid \mathbf{X}, \boldsymbol{w}, \sigma^2) = \mathcal{N}(\boldsymbol{y} \mid \mathbf{X}\boldsymbol{w}, \sigma^2 \mathbf{I}) $$

と書ける。

対数尤度を最大化することは、定数項を除けば、

$$ \| \boldsymbol{y} - \mathbf{X}\boldsymbol{w} \|^2 $$

を最小化することと等価であり、これは 最小二乗法 と一致する。

したがって、

$$ \boldsymbol{w}^\ast = \arg\min_{\boldsymbol{w}} \| \boldsymbol{y} - \mathbf{X}\boldsymbol{w} \|^2 $$

を解くことで、線形回帰モデルのパラメータを推定する。

この目的関数を $\boldsymbol{w}$ で偏微分して、$0$ とおくと、

$$ -2\mathbf{X}^T(\boldsymbol{y} - \mathbf{X}\boldsymbol{w}) = \boldsymbol{0} $$

が得られる。これを整理すると、

$$ \mathbf{X}^T \mathbf{X} \boldsymbol{w} = \mathbf{X}^T \boldsymbol{y} $$

となる。この式は 正規方程式 と呼ばれる。

$\mathbf{X}^T \mathbf{X}$ が正則であると仮定すると、パラメータの最尤推定値は

$$ \boldsymbol{w}^\ast = (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \boldsymbol{y} $$

と求められる。

この推定されたパラメータ $\boldsymbol{w}^\ast$ を用いることで、訓練データに対する 予測値ベクトル は

$$ \hat{\boldsymbol{y}} = \mathbf{X}\boldsymbol{w}^\ast $$

と表される。

すなわち、線形回帰における予測は、

1. 説明変数行列 $\mathbf{X}$ を構成し
2. 正規方程式により $\boldsymbol{w}^\ast$ を求め
3. 行列積 $\mathbf{X}\boldsymbol{w}^\ast$ を計算する

という一連の操作によって得られる。

ここまで長くなってしまったが、この議論により、非線形回帰モデルを基底関数展開によって線形回帰のように扱うための下地が整った。

基底関数展開による非線形回帰の行列表現

先に方針を述べると、前項の記号を

• $\mathbf{X}$ → $\boldsymbol{\Phi}$
• $\tilde{\boldsymbol{x}}$ → $\boldsymbol{\phi}(\boldsymbol{x})$

のように置き換えていくことで、同じ議論ができることを確認する。

入力 $\boldsymbol{x}$ に対して、あらかじめ定めた $M$ 個の基底関数

$$ \phi_1(\boldsymbol{x}), \phi_2(\boldsymbol{x}), \ldots, \phi_M(\boldsymbol{x}) $$

を用い、基底関数ベクトルを、

$$ \boldsymbol{\phi}(\boldsymbol{x}) = (\phi_1(\boldsymbol{x}), \phi_2(\boldsymbol{x}), \ldots, \phi_M(\boldsymbol{x}))^T $$

と定義する。

このとき、基底関数展開を用いた回帰モデルは、

$$ y_i = \boldsymbol{w}^T \boldsymbol{\phi}(\boldsymbol{x}_i) + \varepsilon_i $$

と表される。ここで $\boldsymbol{w}$ は基底関数に対応するパラメータベクトルである。

次に、線形回帰と同様に行列表現を導入する。

各サンプルに対する基底関数ベクトルを行として並べた 基底関数行列（設計行列） $\boldsymbol{\Phi}$ を、

$$ \boldsymbol{\Phi} = \begin{bmatrix} \boldsymbol{\phi}(\boldsymbol{x}_1)^T \\ \boldsymbol{\phi}(\boldsymbol{x}_2)^T \\ \vdots \\ \boldsymbol{\phi}(\boldsymbol{x}_N)^T \end{bmatrix} = \begin{bmatrix} \phi_1(\boldsymbol{x}_1) & \cdots & \phi_M(\boldsymbol{x}_1) \\ \phi_1(\boldsymbol{x}_2) & \cdots & \phi_M(\boldsymbol{x}_2) \\ \vdots & \ddots & \vdots \\ \phi_1(\boldsymbol{x}_N) & \cdots & \phi_M(\boldsymbol{x}_N) \end{bmatrix} $$

と定義する。

このとき、基底関数展開を用いた非線形回帰モデルは、線形回帰と同様に

$$ \boldsymbol{y} = \boldsymbol{\Phi}\boldsymbol{w} + \boldsymbol{\varepsilon} $$

とまとめて表される。

ノイズが独立同分布な正規分布

$$ \varepsilon_i \sim \mathcal{N}(0, \sigma^2) $$

に従うと仮定すると、最尤推定は線形回帰の場合と同様に、次の最小二乗問題に帰着する。

$$ \boldsymbol{w}^\ast = \arg\min_{\boldsymbol{w}} \| \boldsymbol{y} - \boldsymbol{\Phi}\boldsymbol{w} \|^2 $$

この問題の解は、正規方程式

$$ \boldsymbol{\Phi}^T \boldsymbol{\Phi}\boldsymbol{w} = \boldsymbol{\Phi}^T \boldsymbol{y} $$

を満たし、$\boldsymbol{\Phi}^T \boldsymbol{\Phi}$ が正則であれば、

$$ \boldsymbol{w}^\ast = (\boldsymbol{\Phi}^T \boldsymbol{\Phi})^{-1}\boldsymbol{\Phi}^T \boldsymbol{y} $$

として求められる。

この推定されたパラメータ $\boldsymbol{w}^\ast$ を用いることで、訓練データに対する 予測値ベクトル は

$$ \hat{\boldsymbol{y}} = \boldsymbol{\Phi}\boldsymbol{w}^\ast $$

と表される。

例：多項式回帰 $y = 0.5 x^2 − 1.0 x + 2$

例として、多項式回帰 $y = 0.5 x^2 − 1.0 x + 2$ を予測する。

細かい計算式などは既に触れているため、Numpy での実装にのみ集中する。

まず、真のデータ生成過程。

import numpy as np

N = 100
"""サンプル数"""

w2_true = 0.5
"""真の 2 次係数"""

w1_true = -1.0
"""真の 1 次係数"""

b_true = 2.0
"""真の定数項"""

sigma = 0.5
"""ノイズの標準偏差"""

def generate_true_data(x):
    """真のデータを生成する"""
    # 平均 0、標準偏差 sigma の正規分布に従うノイズを生成する
    epsilon = np.random.normal(0, sigma)
    return w2_true * (x**2) + w1_true * x + b_true + epsilon

次に、推定。

import numpy as np

def make_design_matrix_poly(x_array, degree):
    """
    多項式回帰の設計行列を作る（バイアス項込み）

    例: degree=2 -> [1, x, x^2]
    """
    # shape: (N, degree+1)
    return np.vstack([x_array**d for d in range(degree + 1)]).T

def fit_by_normal_equation(Phi, y):
    """正規方程式で w_hat を求める"""

    A = Phi.T @ Phi
    b = Phi.T @ y

    # 逆行列よりも solve の方が数値的に安定する
    if np.linalg.matrix_rank(A) < A.shape[0]:
        raise ValueError("Phi^T Phi が特異（rank 不足）で解けません。データや次数を見直してください。")
    return np.linalg.solve(A, b)

最後に、予測。

import numpy as np
import matplotlib.pyplot as plt

degree = 2
"""多項式の次数 (2 次多項式回帰)"""

def predict_poly(x_array, w, degree):
    """多項式回帰で予測する"""
    Phi_new = make_design_matrix_poly(x_array, degree)
    return Phi_new @ w

# 訓練データを生成する
x_train = np.random.uniform(0.0, 10.0, N)
y_train = np.array([generate_true_data(x) for x in x_train])

# 予測モデルのパラメータを求める
Phi = make_design_matrix_poly(x_train, degree)
w_hat = fit_by_normal_equation(Phi, y_train)

print("True coefficients:")
print(f"  b={b_true:.4f}, w1={w1_true:.4f}, w2={w2_true:.4f}")
print("Estimated coefficients:")
print(f"  b={w_hat[0]:.4f}, w1={w_hat[1]:.4f}, w2={w_hat[2]:.4f}")

# 予測モデルを用いて予測を行う
y_pred_train = predict_poly(x_train, w_hat, degree)

x_line = np.linspace(x_train.min(), x_train.max(), 300)
y_line = predict_poly(x_line, w_hat, degree)

# 結果をプロットする
plt.figure()
plt.scatter(x_train, y_train, label="Training Data")
plt.plot(x_line, y_line, label=f"Polynomial Fit (degree={degree})")
plt.xlabel("x")
plt.ylabel("y")
plt.title("Polynomial Regression (Design Matrix + Normal Equation)")
plt.legend()
plt.show()

ホールドアウト法

ホールドアウト法（hold-out method） とは、データ集合をあらかじめ 訓練データ と テストデータ に分割し、訓練データで学習したモデルをテストデータで評価する手法である。

データ集合を

$$ \mathcal{D} = \mathcal{D}_{\text{train}} \cup \mathcal{D}_{\text{test}}, \quad \mathcal{D}_{\text{train}} \cap \mathcal{D}_{\text{test}} = \varnothing $$

のように分割し、学習は $\mathcal{D}_{\text{train}}$ のみを用いて行う。

その後、テストデータに対する予測値 $\hat{y}$ と真の値 $y$ を比較することで、汎化性能を評価する。

ホールドアウト法は、

• 実装が容易で計算コストが低い
• 大規模データに適している

といった利点がある一方で、

• データの分割方法によって評価結果がばらつく
• データ数が少ない場合、評価が不安定になる

という欠点を持つ。

$k$-分割交差検証法

$k$-分割交差検証法（$k$-fold cross validation） は、データ集合を $k$ 個の部分集合（フォールド）に分割し、そのうち 1 つをテストデータ、残りの $k-1$ 個を訓練データとして学習と評価を繰り返す手法である。

具体的には、各フォールド $j = 1, \ldots, k$ に対して、

1. フォールド $j$ をテストデータとする
2. 残りのデータでモデルを学習する
3. テスト誤差を計算する

という操作を行い、最終的に得られた $k$ 個の評価値の平均をモデル性能とする。

$$ \text{評価値} = \frac{1}{k} \sum_{j=1}^{k} \text{Error}_j $$

この方法により、特定のデータ分割に依存しない、より安定した性能評価が可能となる。

$k$-分割交差検証法は、

• データを有効に活用できる
• 評価の分散が小さくなる

という利点を持つ一方で、

• 学習を $k$ 回行う必要があり計算コストが高い

という欠点がある。

---

検証手法の比較と使い分け

手法	特徴	主な用途
ホールドアウト法	高速・簡単	データが多い場合、初期検証
$k$-分割交差検証法	安定した評価	データが少ない場合、モデル選択

---

ロジスティック回帰

定義

ロジスティック回帰（logistic regression） は、主に 二値分類問題 に用いられる教師あり学習手法である。

説明変数の線形結合を シグモイド関数 によって非線形変換することで、目的変数がクラス 1 に属する確率をモデル化する。

ロジスティック回帰は「回帰」という名称を持つが、連続値を予測する回帰問題ではなく、確率的な分類モデル である点に注意が必要である。

モデルの定式化

各サンプルの説明変数ベクトルを $\boldsymbol{x}_i$ 、重みベクトルを $\boldsymbol{w}$ 、バイアスを $b$ とする。

まず、線形結合を

$$ z_i = \boldsymbol{w}^T \boldsymbol{x}_i + b $$

と定義し、これをシグモイド関数

$$ \sigma(z) = \frac{1}{1 + e^{-z}} $$

に通すことで、クラス 1 に属する確率を

$$ P(y_i = 1 \mid \boldsymbol{x}_i) = \sigma(z_i) $$

とモデル化する。

損失関数と最尤推定

目的変数 $y_i \in \{ 0, 1 \}$ がベルヌーイ分布に従うと仮定すると、尤度は

$$ p(y_i \mid \boldsymbol{x}_i) = \sigma(z_i)^{y_i} (1 - \sigma(z_i))^{1 - y_i} $$

となる。

この対数尤度の 負 を取ったものが、ロジスティック回帰で用いられる 交差エントロピー損失（対数損失） である。

$$ L(\boldsymbol{w}, b) = - \sum_{i=1}^{N} \left[ y_i \log \sigma(z_i) + (1 - y_i)\log(1 - \sigma(z_i)) \right] $$

ロジスティック回帰では、この損失関数を最小化することでパラメータを推定する。

なお、この最適化問題には解析解が存在しないため、勾配降下法 などの反復最適化手法が用いられる。

実装では準ニュートン法（例：LBFGS）が用いられることが多い。

予測

推定されたパラメータ $(\boldsymbol{w}^\ast, b^\ast)$ を用いることで、

• 確率予測

$$ \hat{p}_i = \sigma(\boldsymbol{w}^{\ast T}\boldsymbol{x}_i + b^\ast) $$

• クラス予測

$$ \hat{y}_i = \begin{cases} 1 & (\hat{p}_i \ge 0.5) \\ 0 & (\hat{p}_i < 0.5) \end{cases} $$

が得られる。

例：タイタニックの乗客データ

以下では、scikit-learn による実装例を示す。

細かい前処理の数学的説明は省略し、実装の流れに集中する。

まず、インポート。

import numpy as np
import pandas as pd

from sklearn.datasets import fetch_openml
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score, confusion_matrix

前処理。

# Titanic データセットの読み込み
titanic = fetch_openml("titanic", version=1, as_frame=True)
df = titanic.frame

# 使用する特徴量と目的変数
features = ["pclass", "sex", "age", "fare"]
target = "survived"

# 欠損値処理（簡易的）
df = df[features + [target]].dropna()

# カテゴリ変数の処理（sex を 0/1 に変換）
df["sex"] = df["sex"].map({"male": 0, "female": 1})

X = df[features]
y = df[target].astype(int)

次に、学習データとテストデータに分割する。

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)

# ロジスティック回帰モデル
model = LogisticRegression(max_iter=1000)
model.fit(X_train_scaled, y_train)

最後に、予測を行ってモデルの性能評価をする。

# 予測
y_pred = model.predict(X_test_scaled)

# 精度評価
accuracy = accuracy_score(y_test, y_pred)
cm = confusion_matrix(y_test, y_pred)

print(f"Accuracy: {accuracy:.4f}")
print("Confusion Matrix:")
print(cm)

実行結果は以下。

Accuracy: 0.7420
Confusion Matrix:
[[142  33]
 [ 48  91]]

• 混同行列

	予測: 死亡 (0)	予測: 生存 (1)
実際: 死亡 (0)	142	33
実際: 生存 (1)	48	91

• Accuracy（正解率）

Accuracy は

$$ \text{Accuracy} = \frac{142 + 91}{142 + 33 + 48 + 91} = \frac{233}{314} \approx 0.742 $$

全体の 約 74% を正しく分類できているため、モデルとしては割と良いように見える。

• Recall（再現率）

$$ \text{Recall} = \frac{\text{TP}}{\text{TP} + \text{FN}} = \frac{91}{91 + 48} \approx 0.655 $$

実際に生存した人のうち、約 65.5% しか正しく生存と予測できていない。

裏を返せば、約 3 割以上の生存者を「死亡」と見誤っている と言える。

• Precision（適合率）

$$ \text{Precision} = \frac{\text{TP}}{\text{TP} + \text{FP}} = \frac{91}{91 + 33} \approx 0.734 $$

「生存」と予測した場合、その 約 73% は正しかった。

生存予測の信頼性はそこそこ高い。

• 性能評価

死亡の正解率が $142 / 175 \approx 81%$ と高いが、生存の見逃し (48人) が比較的多い。これは、全体的に 死亡判定に寄る モデルと言える。

これは、このデータセットには 死亡者の方が多い ため、極論、何も考えずに「全員死亡」と予測しても Accuracy はそこそこの値が出てしまうからだと考えられる。

そのため、$\text{Accuracy} \approx 0.742$ は悪くないが、生存者の $\text{Recall} \approx 0.655$ は改善の余地がある、と評価できる。

本課題では「生存者を見逃さない」ことが重要であるため、Recall を重視した評価や、閾値の調整が今後の課題として考えられる。

• 課題「年齢が 30 歳で男の乗客は生き残れるか？」

# 仮定する乗客の条件
# pclass, sex, age, fare は学習時と同じ特徴量
sample_passenger = pd.DataFrame({
    "pclass": [3],   # 3等客室を仮定
    "sex": [0],      # male = 0
    "age": [30],     # 年齢 30 歳
    "fare": [30.0]   # 平均的な運賃を仮定
})

# 学習時と同じ標準化を適用
sample_scaled = scaler.transform(sample_passenger)

# 生存確率を予測
survival_prob = model.predict_proba(sample_scaled)[0, 1]

# クラス予測（0 or 1）
survival_pred = model.predict(sample_scaled)[0]

print(f"Survival probability: {survival_prob:.4f}")

if survival_pred == 1:
    print("→ この乗客は『生存する』と予測されます。")
else:
    print("→ この乗客は『生存しない』と予測されます。")

実行結果は以下。

Survival probability: 0.0773
→ この乗客は『生存しない』と予測されます。

自分の解釈としては、「男性」「低等級客室」の乗客の生存率が低いという傾向にあるため、生存確率が低く算出されたと考えられる。

また、先の性能評価から、死亡判定に寄った予測に見える。

---

主成分分析

主成分分析（PCA: Principal Component Analysis） とは、多次元データの分散をできるだけ保持したまま、より低次元の空間へデータを写像する 次元削減手法 である。

高次元データに含まれる相関の強い変数を統合し、情報量（分散）を最大化する新たな軸（主成分）を求めることで、

• データの可視化
• ノイズの除去
• 計算コストの削減

などを目的として用いられる。

数学的定式化

サンプル数を $N$ 、特徴量の次元を $m$ とする。

各サンプルの説明変数ベクトルを $\boldsymbol{x}_i = (x_{i1}, x_{i2}, \ldots, x_{im})^T \ (i = 1, \ldots, N)$ とする。

これをそのまま並べるのではなく、各要素から対応する特徴量の平均を引いて、平均を $0$ にする（中心化する）。

第 $j$ 特徴量の平均を

$$ \bar{x}_j = \frac{1}{N} \sum_{i=1}^N x_{ij} \quad (j = 1, \ldots, m) $$

と定義する。

各要素から対応する特徴量の平均を引くと、

$$ \tilde{x}_{ij} = x_{ij} - \bar{x}_j $$

となるので、これを各要素とする 中心化データ行列 は、

$$ \mathbf{X} = \begin{pmatrix} x_{11} - \bar{x}_1 & x_{12} - \bar{x}_2 & \cdots & x_{1m} - \bar{x}_m \\ x_{21} - \bar{x}_1 & x_{22} - \bar{x}_2 & \cdots & x_{2m} - \bar{x}_m \\ \vdots & \vdots & \ddots & \vdots \\ x_{N1} - \bar{x}_1 & x_{N2} - \bar{x}_2 & \cdots & x_{Nm} - \bar{x}_m \end{pmatrix} $$

となる。

このとき、共分散行列は

$$ \mathbf{\Sigma} = \frac{1}{N} \mathbf{X}^T \mathbf{X} $$

で与えられる。

なお、分母を $N$ とする定義と $N-1$ とする定義があるが、本稿では簡単のため $1/N$ を用いている。

主成分分析では、この共分散行列の 固有値問題

$$ \mathbf{\Sigma}\boldsymbol{u}_k = \lambda_k \boldsymbol{u}_k $$

を解くことで、

• 固有値 $\lambda_k$：分散の大きさ
• 固有ベクトル $\boldsymbol{u}_k$：主成分方向

を求める。

固有値の大きい順に主成分を並べ、上位 $d$ 個の主成分を用いることで、次元削減を行う。

主成分への射影

主成分行列を

$$ \mathbf{U}_d = (\boldsymbol{u}_1, \boldsymbol{u}_2, \ldots, \boldsymbol{u}_d) $$

とすると、元のデータは

$$ \boldsymbol{z}_i = \mathbf{U}_d^T \tilde{\boldsymbol{x}}_i $$

として、$d$ 次元空間に射影される。

このとき、$\boldsymbol{z}_i$ が 主成分得点である。

例：乳がん検査データ

以下では、scikit-learn による実装例を示す。

まず、インポート。

import numpy as np
import pandas as pd

from sklearn.datasets import load_breast_cancer
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA

データセットの読み込み。

data = load_breast_cancer()
X = data.data  # 説明変数
y = data.target  # 0: malignant, 1: benign

print(X.shape)  # (569, 30)

このデータセットは、

• サンプル数：569
• 特徴量数：30

からなる高次元データである。

次に、標準化する。

主成分分析は 分散に基づく手法であるため、各特徴量のスケールを揃えることが重要である。

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

ここから主成分分析を行うが、課題にある通り 2 次元に次元圧縮を行う。

pca = PCA(n_components=2)
X_pca = pca.fit_transform(X_scaled)

print(X_pca.shape)  # (569, 2)

選んだ 2 個の成分が適切か確認するため、分散説明率を見る。

分散説明率とは、各主成分が元データ全体の分散のうち、どの程度を説明しているかを表す指標である。数式は省略。

主成分分析では、分散が大きい方向ほど重要な情報を含むと考えられるため、分散説明率を用いて次元削減後に保持される情報量を評価する。

explained_variance_ratio = pca.explained_variance_ratio_

print("Explained variance ratio:")
print(explained_variance_ratio)
print("Cumulative explained variance:")
print(np.cumsum(explained_variance_ratio))

結果。

Explained variance ratio:
[0.44272026 0.18971182]
Cumulative explained variance:
[0.44272026 0.63243208]

第 1 成分で 44% もあり、累積で約 63% あるので、十分そう。

最後に、結果をプロットする。

import matplotlib.pyplot as plt

plt.figure()
plt.scatter(X_pca[y == 0, 0], X_pca[y == 0, 1], label="Malignant", alpha=0.7)
plt.scatter(X_pca[y == 1, 0], X_pca[y == 1, 1], label="Benign", alpha=0.7)

plt.xlabel("First Principal Component")
plt.ylabel("Second Principal Component")
plt.title("PCA of Breast Cancer Dataset")
plt.legend()
plt.show()

次元をかなり減らしたが、それでも Malignant (悪性) と Benign (良性) はある程度分離できている。

これは、乳がん検査データにおいて、少数の主成分が重要な情報を集約している ことを示している。

なお、主成分分析は教師なし学習であり、分類精度そのものを最適化する手法ではない点に注意が必要である。

---

k 近傍法

定義

k近傍法 (k-NN) は、入力 $\boldsymbol{x}$ に対して、訓練データの中から距離が近い順に $k$ 個のサンプル（近傍）を選び、そのラベルの多数決によりクラスを予測する分類手法である。

パラメータを学習するというより、訓練データを保持し、予測時に参照する（lazy learning） 点が特徴である。

計算（距離と多数決）

訓練データ $(\boldsymbol{x}_i, y_i) \ (i=1,\dots,N)$ と、予測したい入力 $\boldsymbol{x}$ があるとする。

ここでは代表的なユークリッド距離に基づき、距離が小さい順に $k$ 個の近傍集合 $N_k(\boldsymbol{x})$ を求める。

$$ d(\boldsymbol{x}, \boldsymbol{x}_i) = \| \boldsymbol{x} - \boldsymbol{x}_i \| $$

二値分類 ( $y \in \{0,1\}$ ) の場合、近傍のラベルの多数決により

$$ \hat{y} = \mathrm{mode} \{ y_i \mid i \in N_k(\boldsymbol{x}) \} $$

として予測する。

なお、距離比較では平方根を取らない 二乗距離 でも大小関係が同じであるため、

$$ d^2(\boldsymbol{x}, \boldsymbol{x}_i)=\sum_j (x_j - x_{ij})^2 $$

を用いて計算を簡略化できる。

例

課題にあった Jupyter Notebook を少し修正して、k 近傍法による分類と決定領域を可視化する。

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

def gen_data(seed=42):
    rng = np.random.default_rng(seed)
    x0 = rng.normal(size=50).reshape(-1, 2) - 1.0
    x1 = rng.normal(size=50).reshape(-1, 2) + 1.0
    X_train = np.concatenate([x0, x1])
    y_train = np.concatenate([np.zeros(25), np.ones(25)]).astype(int)
    return X_train, y_train

def distance(x1, X2):
    return np.sum((X2 - x1) ** 2, axis=1)

def knc_predict(n_neighbors, X_train, y_train, X_test):
    y_pred = np.empty(len(X_test), dtype=y_train.dtype)
    for i, x in enumerate(X_test):
        distances = distance(x, X_train)
        nearest_index = np.argsort(distances)[:n_neighbors]
        mode = stats.mode(y_train[nearest_index], keepdims=True).mode[0]
        y_pred[i] = int(mode)
    return y_pred

def plt_result(X_train, y_train, y_pred, grid_min=-5, grid_max=5, grid_n=100):
    xx0, xx1 = np.meshgrid(
        np.linspace(grid_min, grid_max, grid_n),
        np.linspace(grid_min, grid_max, grid_n),
    )
    plt.scatter(X_train[:, 0], X_train[:, 1], c=y_train)
    plt.contourf(
        xx0,
        xx1,
        y_pred.reshape(grid_n, grid_n).astype(float),
        alpha=0.2, levels=np.linspace(0, 1, 3),
    )
    plt.xlabel("x1")
    plt.ylabel("x2")
    plt.title("k-NN Decision Region")
    plt.show()

# 実行
X_train, y_train = gen_data(seed=42)

n_neighbors = 3
xx0, xx1 = np.meshgrid(np.linspace(-5, 5, 100), np.linspace(-5, 5, 100))
X_test = np.array([xx0, xx1]).reshape(2, -1).T

y_pred = knc_predict(n_neighbors, X_train, y_train, X_test)
plt_result(X_train, y_train, y_pred)

この結果、訓練データの分布に応じて、非線形な決定境界が得られる。

一般に、$k$ を小さくすると境界は複雑になり過学習しやすく、$k$ を大きくすると滑らかになり汎化しやすい傾向がある。

また、scikit-learn だと以下の通り。

from sklearn.neighbors import KNeighborsClassifier

# 訓練データの生成（前節と同じ）
X_train, y_train = gen_data(seed=42)

# k の設定
k = 3

# sklearn の k-NN 分類器
knn_sklearn = KNeighborsClassifier(
    n_neighbors=k,
    metric="euclidean",
)

knn_sklearn.fit(X_train, y_train)

# グリッド点の生成
xx0, xx1 = np.meshgrid(
    np.linspace(-5, 5, 100),
    np.linspace(-5, 5, 100),
)
X_test = np.array([xx0, xx1]).reshape(2, -1).T

# sklearn による予測
y_pred_sklearn = knn_sklearn.predict(X_test)

# 可視化
plt.figure()
plt.scatter(X_train[:, 0], X_train[:, 1], c=y_train)
plt.contourf(
    xx0,
    xx1,
    y_pred_sklearn.reshape(100, 100).astype(float),
    alpha=0.2, levels=np.linspace(0, 1, 3),
)
plt.xlabel("x1")
plt.ylabel("x2")
plt.title("k-NN Decision Region (scikit-learn)")
plt.show()

結果は同じようなので、Numpy での実装は上手くできていることが分かった。

実際は、ライブラリを使用するのが望ましいだろう。

---

k-means

定義

k-means（k平均法） は、教師なし学習における代表的な クラスタリング手法である。

データ点を $K$ 個のクラスタに分割し、各クラスタの代表点（重心）を更新しながら、クラスタ内のばらつきを小さくするように反復最適化を行う。

ただし、k-means は初期値に依存し、大域最適解が保証されるわけではない。

数学的定式化

データ点を $\boldsymbol{x}_i \ (i=1,\dots,N)$ とし、クラスタ数を $K$ とする。

k-means は、各点の割り当て $c_i \in {1,\dots,K}$ と各クラスタの重心 $\boldsymbol{\mu}_k$ を変数として、目的関数（クラスタ内平方和）を最小化する。

$$ J = \sum_{i=1}^{N} \left\| \boldsymbol{x}_i - \boldsymbol{\mu}_{c_i} \right\|^2 $$

計算（アルゴリズム）

k-means は次の 2 ステップを反復する。

1. 割り当て (Assignment)

各点 $\boldsymbol{x}_i$ を最も近い重心へ割り当てる：

$$ c_i = \arg \min_{k} \left\| \boldsymbol{x}_i - \boldsymbol{\mu}_k \right\|^2 $$

2. 重心更新 (Update)

各クラスタの重心を、そのクラスタに属する点の平均で更新する：

$$ \boldsymbol{\mu}_k = \frac{1}{|C_k|} \sum_{i: c_i=k} \boldsymbol{x}_i $$

この 2 ステップを収束（割り当てが変わらないなど）まで繰り返す。

例

インポート。

import numpy as np
import matplotlib.pyplot as plt

データ生成。

def gen_data(seed=42, n_per_cluster=100):
    rng = np.random.default_rng(seed)
    x1 = rng.normal(size=(n_per_cluster, 2)) + np.array([-5.0, -5.0])
    x2 = rng.normal(size=(n_per_cluster, 2)) + np.array([ 5.0, -5.0])
    x3 = rng.normal(size=(n_per_cluster, 2)) + np.array([ 0.0,  5.0])
    return np.vstack((x1, x2, x3))

# データ作成
X_train = gen_data(seed=42, n_per_cluster=100)

# データ描画
plt.figure()
plt.scatter(X_train[:, 0], X_train[:, 1], alpha=0.7)
plt.title("Generated Data")
plt.xlabel("x1")
plt.ylabel("x2")
plt.show()

NumPy で k-means を作成する。

def kmeans_numpy(X, n_clusters=3, seed=42, max_iter=100, tol=1e-6):
    """
    NumPyでk-means
    - X: (N, d)
    - returns: labels (N,), centers (K, d), inertia (float)
    """
    rng = np.random.default_rng(seed)
    N, d = X.shape
    K = n_clusters

    # 初期中心：データ点からK個ランダムに選ぶ
    centers = X[rng.choice(N, size=K, replace=False)].copy()
    labels = np.full(N, -1, dtype=int)

    for _ in range(max_iter):
        # 距離^2行列 (N, K)
        dists2 = np.sum((X[:, None, :] - centers[None, :, :]) ** 2, axis=2)

        new_labels = np.argmin(dists2, axis=1)

        # 収束：割当が変わらない
        if np.array_equal(new_labels, labels):
            labels = new_labels
            break
        labels = new_labels

        # 中心更新
        new_centers = np.zeros_like(centers)
        for k in range(K):
            members = X[labels == k]
            if len(members) == 0:
                # 空クラスタ対策：ランダム点で再初期化
                new_centers[k] = X[rng.integers(0, N)]
            else:
                new_centers[k] = members.mean(axis=0)

        # 中心移動が小さければ収束
        shift = np.linalg.norm(new_centers - centers)
        centers = new_centers
        if shift < tol:
            break

    inertia = np.sum((X - centers[labels]) ** 2)
    return labels, centers, inertia

実行。

n_clusters = 3
labels_np, centers_np, inertia_np = kmeans_numpy(X_train, n_clusters=n_clusters, seed=42)

print("NumPy inertia:", inertia_np)
print("NumPy centers:\n", centers_np)

NumPy inertia: 566.527212719175
NumPy centers:
 [[ 5.05883416 -5.01867893]
 [-5.0147462  -5.04615525]
 [-0.11107121  4.97482521]]

プロットする。

def plot_kmeans_result(X, centers, grid_min=-10, grid_max=10, grid_n=200, title="k-means"):
    xx0, xx1 = np.meshgrid(
        np.linspace(grid_min, grid_max, grid_n),
        np.linspace(grid_min, grid_max, grid_n),
    )
    grid = np.c_[xx0.ravel(), xx1.ravel()]

    # グリッド点を最近傍中心に割当
    dists2 = np.sum((grid[:, None, :] - centers[None, :, :]) ** 2, axis=2)
    pred = np.argmin(dists2, axis=1).reshape(xx0.shape)

    plt.figure()
    plt.contourf(xx0, xx1, pred, alpha=0.2, levels=np.arange(len(centers)+1)-0.5)
    plt.scatter(X[:, 0], X[:, 1], alpha=0.7)
    plt.scatter(centers[:, 0], centers[:, 1], marker="x", s=200)
    plt.title(title)
    plt.xlabel("x1")
    plt.ylabel("x2")
    plt.show()

plot_kmeans_result(X_train, centers_np, title="k-means (NumPy)")

こちらも scikit-learn と比較する。

from sklearn.cluster import KMeans

kmeans = KMeans(n_clusters=3, random_state=42, n_init=10)
labels_sk = kmeans.fit_predict(X_train)
centers_sk = kmeans.cluster_centers_
inertia_sk = kmeans.inertia_

print("sklearn inertia:", inertia_sk)
print("sklearn centers:\n", centers_sk)

plot_kmeans_result(X_train, centers_sk, title="k-means (scikit-learn)")

良さそう。

---

パターン認識

パターン認識（pattern recognition） とは、データの特徴量に基づいて、対象を分類・識別・検索するための枠組みである。

距離計算

パターン認識においては、データ点同士の「似ている・似ていない」を数値化するために 距離（distance） を定義する。

以下では、代表的な距離指標を示す。

ユークリッド距離（Euclidean distance）

2 つのベクトル $\boldsymbol{x}, \boldsymbol{y} \in \mathbb{R}^m$ に対して、

$$ d_{\text{Euc}}(\boldsymbol{x}, \boldsymbol{y}) = \sqrt{\sum_{i=1}^{m} (x_i - y_i)^2} $$

で定義される。

• 直感的で最も一般的
• k-NN、k-means などで広く用いられる
• スケールの影響を受けやすい

マンハッタン距離（Manhattan distance）

$$ d_{\text{Man}}(\boldsymbol{x}, \boldsymbol{y}) = \sum_{i=1}^{m} |x_i - y_i| $$

• 格子状の移動距離に相当
• 外れ値の影響がユークリッド距離より小さい場合がある

Lp 距離（Minkowski distance）

ユークリッド距離やマンハッタン距離を一般化した距離で、

$$ d_p(\boldsymbol{x}, \boldsymbol{y}) = \left( \sum_{i=1}^{m} |x_i - y_i|^p \right)^{1/p} \quad (p \ge 1) $$

• $p=1$：マンハッタン距離
• $p=2$：ユークリッド距離

として、それぞれが特別な場合に対応する。

コサイン距離（Cosine distance）

ベクトルの 向き に注目した距離であり、

$$ \text{cosine similarity} = \frac{\boldsymbol{x}^T \boldsymbol{y}} {\| \boldsymbol{x} \| \| \boldsymbol{y} \|} $$

を用いて、

$$ d_{\text{cos}}(\boldsymbol{x}, \boldsymbol{y}) = 1 - \frac{\boldsymbol{x}^T \boldsymbol{y}} {\| \boldsymbol{x} \| \| \boldsymbol{y} \|} $$

と定義される。

• 文書ベクトルや特徴量の大きさが重要でない場合に有効
• 情報検索・自然言語処理で頻用される

マハラノビス距離（Mahalanobis distance）

特徴量間の相関を考慮した距離であり、共分散行列 $\mathbf{\Sigma}$ を用いて、

$$ d_{\text{Mah}}(\boldsymbol{x}, \boldsymbol{y}) = \sqrt{(\boldsymbol{x} - \boldsymbol{y})^T \mathbf{\Sigma}^{-1} (\boldsymbol{x} - \boldsymbol{y})} $$

と定義される。

• 分散の大きい方向を抑制し、相関を考慮できる
• PCA やガウス分布を仮定する手法と相性が良い
• 共分散行列の推定が必要

kd-tree

kd-tree（k-dimensional tree） は、多次元空間の点を効率的に管理するための 空間分割データ構造である。

• 各ノードで 1 次元を選び、データを二分割する
• 再帰的に空間を分割することで木構造を構成する

kd-tree を用いることで、全探索では $O(N)$ かかる最近傍探索を、平均的にはより高速に行える。

ただし、

• 次元が高くなると性能が低下する（次元の呪い）

という欠点がある。

近似最近傍探索（Approximate Nearest Neighbor）

高次元データや大規模データセットでは、厳密な最近傍探索が計算量的に困難になる。

この問題に対して、近似最近傍探索（ANN） が用いられる。

• 厳密な最近傍でなく、「十分近い点」を高速に探索
• 精度と計算量のトレードオフを許容

代表的な考え方には、

• kd-tree の枝刈り
• 局所性鋭敏ハッシュ（LSH）
• グラフベース探索

などがある。