Python3ではじめるシステムトレード：判別分析

判別分析は歴史が古くフィッシャーの線形判別分析が起源です。その後多くの一般化が行われました。そのため、内容が難しく感じ、混乱も生じます。統計学実践ワークブックを参考にまとめてみました。PythonのコードはChatGPTの力をかなり借りています。本記事は学習過程のメモとして書かれています。間違い等がありましたら、ご指摘いただけると助かります。

勉強会をしています。

判別分析

1：５０

ーフィッシャーの判別分析ー

フィッシャーの判別分析（LDA）は、クラスラベル付データを用いて、それらの情報を標本平均,標本分散共分散行列に集約し、クラス（または群）の情報を明示的に使用して判別に適した射影軸を選択することで,未知のデータのクラスを判別します。

べクトル$\boldsymbol{w}\in \mathbb{R}^p$による射影を考え,射影後のサンプルを

$y_i^{(j)}=\boldsymbol{w}^T\boldsymbol{x}_i^{(j)},j=1,2,i=1,\cdots,n$

とします。一般に「重みベクトル」と呼ばれるベクトル $\boldsymbol{w}^T$はデータポイント（ベクトル）を新しい座標・系に「移す」または「変換する」ための重要な操作です。この新しい座標系は、データポイントが異なる群に分類されるように選ばれます。その際に、それらの値が各群の「中心(重心)」を代表する点として標本平均（または平均ベクトル）を用います。

この1次元の標本$y_i$から射影後の空間を考えると,求めた群$G_1$,$G_2$の標本平均は,群1,2に属する標本

\{\boldsymbol{x}_i^{(j)}\}_{i=0}^{n_j}

の平均べクトルを

\boldsymbol{\bar{x}}_{(j)}=\frac{1}{n_j}\sum_{i=1}^{n_j}\boldsymbol{x}_i^{(j)}

として,

\bar{y}^{(j)}=\frac{1}{n_j}\sum_{i=1}^{n_j}y_i^{(j)}=\boldsymbol{w}^T\boldsymbol{\bar{x}}^{(j)}, j=1,2

となります。2つの群の平均は離れていたほうがよく,一方で,各群のなかでの分散は小さいほうがよいので,以下の基準を最大にする射影方向を求めます。

\lambda(\boldsymbol{w})=\frac{群間分散}{群内分散}

ここで,群間分散は

(\bar{y}^{(1)}-\bar{y}^{(2)})^2=(\boldsymbol{w}^T(\boldsymbol{\bar{x}}^{(1)}-\boldsymbol{\bar{x}}^{(2)}))^2

で定義します。群内分散は

\frac{1}{n_j-1}\sum _{i=1}^{n_j}(y_j^{(j)}-\bar{y}^{(j)})^2=\boldsymbol{w}^TS_j\boldsymbol{w}

で定義します。ただし$S_j$は群$G_j$に含まれる標本${ x_i^{(j)}}_{i=1}^{n_j}$の標本分散共分散行列です。

群内分散は未知ですので、分散の不偏推定量を用いて共通の母分散を$\sigma_1^2=\sigma_2^2=\sigma^2$とします。共通の標本分散共分散行列

S=\frac{1}{n_1+n_2-2}((n_1-1)S_1+(n_2-1)S_2)

を用いて$\boldsymbol{w}^TS\boldsymbol{w}$で定義します。$S$の算出には２つの不偏分散を使っているために自由度を２つ使います。よって、評価関数

\lambda(\boldsymbol{w})=\frac{\boldsymbol{w}^T(\boldsymbol{\bar{x}}^{(1)} - \boldsymbol{\bar{x}}^{(2)})^2}{\boldsymbol{w}^TS\boldsymbol{w}}

を最大にする射影方向を求めればよく,最適な射影方向は,

\hat{\boldsymbol{w}}=S^{-1}(\boldsymbol{\bar{x}}^{(1)}-\boldsymbol{\bar{x}}^{(2)})

で与えられます。

import numpy as np
import matplotlib.pyplot as plt

# 2群のデータ点を生成（各群で10個の2D点）
n1 = n2 = 10
x1 = np.random.randn(n1, 2) + np.array([0, 0])
x2 = np.random.randn(n2, 2) + np.array([1, 1])

# 各群の平均を計算
mean1 = np.mean(x1, axis=0)
mean2 = np.mean(x2, axis=0)

# 各群の共分散を計算
cov1 = np.cov(x1, rowvar=False)
cov2 = np.cov(x2, rowvar=False)

# 群内共分散行列（Sw）を計算
S = (cov1*(n1-1) + cov2*(n2-1))/ (n1 + n2 - 2)

# 重みベクトル（w）を計算
w = np.linalg.inv(S).dot(mean1 - mean2)

# 元のデータ点をプロット
plt.figure(figsize=(3, 3))
plt.scatter(x1[:, 0], x1[:, 1], c='b', label='x 1')
plt.scatter(x2[:, 0], x2[:, 1], c='g', label='x 2')

# 射影軸の端点を計算
axis_start = np.min(np.vstack([x1.dot(w), x2.dot(w)])) * w / np.dot(w, w) 
axis_end = np.max(np.vstack([x1.dot(w), x2.dot(w)])) * w / np.dot(w, w) 

# 射影軸をプロット
plt.plot([axis_start[0], axis_end[0]], [axis_start[1], axis_end[1]]\
         , 'r-', label='Projection axis')
plt.xlabel('features 1')
plt.ylabel('features 2')
plt.legend(bbox_to_anchor=(1, 1), loc='upper left')
plt.show()

標本$y$は$y=\boldsymbol{w}^T\boldsymbol{x}$として計算され、スカラーです。この値はベクトル$\boldsymbol{w}$の延長線上に射影された点の長さを示しています。この点の標準基底上の位置は$y$と$\boldsymbol{w}$の内積から求めます。実際にその様子を見てみましょう。

# データ点を射影
y1 = x1.dot(w)
y2 = x2.dot(w)

# 射影軸の端点を計算
axis_start = np.min(np.vstack([y1, y2])) * w
axis_end = np.max(np.vstack([y1, y2])) * w

# 射影軸をプロット
plt.figure(figsize=(4, 4))
plt.plot([axis_start[0], axis_end[0]], [axis_start[1]\
                    , axis_end[1]], 'r-', label='Projection axis')

# 元のデータ点をプロット
plt.scatter(x1[:, 0], x1[:, 1], c='b', label='x 1')
plt.scatter(x2[:, 0], x2[:, 1], c='g', label='x 2')

# 射影されたデータ点をプロット
for xx, yy in zip(x1, y1):
    plt.plot([xx[0], w[0]*yy], [xx[1], w[1]*yy], 'b--')

for xx, yy in zip(x2, y2):
    plt.plot([xx[0], w[0]*yy], [xx[1], w[1]*yy], 'g--')

plt.axis('equal')
plt.grid(True)
plt.xlabel('features 1')
plt.ylabel('features 2')
plt.legend(bbox_to_anchor=(1, 1), loc='upper left')
plt.show()

新しいサンプル$\boldsymbol{x_{new}}$は$y=\boldsymbol{\boldsymbol{\hat{w}}}^T\boldsymbol{x}$と同じ形で射影されます.このサンプルがどちらの群
に属するかを定めるには,射影軸上の各群に属するサンプルの標本平均を求め,その中点

$(\boldsymbol{\hat{w}}^T\boldsymbol{\bar{x}}_1+\boldsymbol{\hat{w}}^T\boldsymbol{\bar{x}}_2)\frac{1}{2}=(\boldsymbol{\bar{x}}^{(1)}-\boldsymbol{\bar{x}}^{(2)})^TS^{-1}(\boldsymbol{\bar{x}}^{(1)}+\boldsymbol{\bar{x}}^{(2)})\frac{1}{2}$

と$y_{new}$を比較すればよいことになります.
すなわち,

$f(\boldsymbol{x_{new}})=\boldsymbol{\hat{w}}^T\boldsymbol{x_{new}}-\frac{1}{2}(\boldsymbol{\bar{x}}^{(1)}-\boldsymbol{\bar{x}}^{(2)})^TS^{-1}(\boldsymbol{\bar{x}}^{(1)}+\boldsymbol{\bar{x}}^{(2)})=\boldsymbol{\hat{w}}^T[\boldsymbol{x_{new}}-\frac{1}{2}(\boldsymbol{\bar{x}}^{(1)}+\boldsymbol{\bar{x}}^{(2)})]$

の値が正ならば$\boldsymbol{x}_{new}$を群$G_1$, 負ならば群$G_2$に分類すればよいのです.この関数$f(\boldsymbol{x})$
をフィッシャーの線形判別関数と呼びます.

判別関数 $ f(\boldsymbol{x}) $ が基本的に $ \boldsymbol{w}^T \boldsymbol{x} + \boldsymbol{b} $ の形で与えられていると仮定すると、 $ w $ は射影ベクトル（重みベクトル）であり、$\boldsymbol{b}$はバイアス項です。$-\frac{1}{2}(\boldsymbol{\bar{x}}^{(1)}-\boldsymbol{\bar{x}}^{(2)})^TS^{-1}(\boldsymbol{\bar{x}}^{(1)}+\boldsymbol{\bar{x}}^{(2)})$はバイアス項です。

nn1 = nn2 = 5
xn1 = np.random.randn(nn1, 2) + np.array([0, 0])
xn2 = np.random.randn(nn2, 2) + np.array([1, 1])

# データ点を射影
yn1 = xn1.dot(w)
yn2 = xn2.dot(w)

# 射影軸の端点を計算
axis_start = np.min(np.r_[y1, y2, yn1, yn2]) * w
axis_end = np.max(np.r_[y1, y2, yn1, yn2]) * w

# 射影軸をプロット
plt.figure(figsize=(4, 4))
plt.plot([axis_start[0], axis_end[0]], [axis_start[1]\
                    , axis_end[1]], 'r-', label='Projection axis')

# 射影されたデータ点をプロット
plt.scatter(yn1*w[0],yn1*w[1],c='b', label='xn 1')
plt.scatter(yn2*w[0],yn2*w[1],c='r', label='xn 2')
plt.scatter((mean1.dot(w)*w[0]+mean2.dot(w)*w[0])/2,(mean1.dot(w)*w[1]+mean2.dot(w)*w[1])/2\
            ,marker="|",s=400,label='descriminant')
plt.xlabel('features 1')
plt.ylabel('features 2')
plt.legend(bbox_to_anchor=(1, 1), loc='upper left')
plt.show()

D^2=\lambda(\boldsymbol{\hat{w}})=(\boldsymbol{x}^{(1)}-\boldsymbol{\bar{x}}^{(2)})^TS^{-1}(\boldsymbol{x}^{(1)}-\boldsymbol{\bar{x}}^{(2)})

をマハラノビス平方距離と呼びます

$\hat{\boldsymbol{w}}=[S^{-1}(\boldsymbol{\bar{x}}^{(1)}-\boldsymbol{\bar{x}}^{(2)})$を

\lambda(\boldsymbol{w})=\frac{[\boldsymbol{w}^T(\boldsymbol{\bar{x}}^{(1)} - \boldsymbol{\bar{x}}^{(2)})]^2}{\boldsymbol{w}^TS\boldsymbol{w}} 

に代入すると

\lambda(\boldsymbol{w})=\frac{\{[S^{-1}(\boldsymbol{\bar{x}}^{(1)}-\boldsymbol{\bar{x}}^{(2)})]^T(\boldsymbol{\bar{x}}^{(1)} - \boldsymbol{\bar{x}}^{(2)})\}^2}{[S^{-1}(\boldsymbol{\bar{x}}^{(1)}-\boldsymbol{\bar{x}}^{(2)})]^TSS^{-1}(\boldsymbol{\bar{x}}^{(1)}-\boldsymbol{\bar{x}}^{(2)})}=

[S^{-1}(\boldsymbol{\bar{x}}^{(1)}-\boldsymbol{\bar{x}}^{(2)})]^T(\boldsymbol{\bar{x}}^{(1)} - \boldsymbol{\bar{x}}^{(2)})=

(\boldsymbol{\bar{x}}^{(1)}-\boldsymbol{\bar{x}}^{(2)})^TS^{-1}(\boldsymbol{\bar{x}}^{(1)} - \boldsymbol{\bar{x}}^{(2)})

となります。

フィッシャーの判別分析 (FDA) は、異なるクラス間の距離を最大化し、同じクラス内の距離を最小化するような射影ベクトルを見つけることを目的としています。この射影ベクトルは、クラスの平均値の差と、クラス内の分散を考慮して求められます。

一方、マハラノビス距離は、データの分布（共分散行列）を考慮した距離を計算するための手法です。この距離は、特定の共分散構造を持つデータセットにおける、2つの点または点と中心との間の「距離」を示します。

共通の標本分散共分散行列$S$を用いて新しいサンプル$\boldsymbol{x}$から群$G_j$の標本平均べクトル$\boldsymbol{\bar{x}}^{(j)}$までのマハラノビス平方距離$D_j^2:=(\boldsymbol{x}-\boldsymbol{\bar{x}}^{(j)})^TS^{-1}(\boldsymbol{x}-\boldsymbol{\bar{x}}^{(j)})$を計算すると,$D_2^2-D_1^2=2f(\boldsymbol{x})$なので,$D_2^2\ge D_1^2$なら$G_1$に,$D_2^2<D_1^2$な
ら$G_2$にサンプルを分類することと等価です。

D_2^2-D_1^2=(\boldsymbol{x}-\boldsymbol{\bar{x}}^{(2)})^TS^{-1}(\boldsymbol{x} - \boldsymbol{\bar{x}}^{(2)})-(\boldsymbol{x}-\boldsymbol{\bar{x}}^{(1)})^TS^{-1}(\boldsymbol{x} - \boldsymbol{\bar{x}}^{(1)})=

(\boldsymbol{x}^T-\boldsymbol{\bar{x}}^{(2)T})S^{-1}(\boldsymbol{x} - \boldsymbol{\bar{x}}^{(2)})-(\boldsymbol{x^T}-\boldsymbol{\bar{x}}^{(1)T})S^{-1}(\boldsymbol{x} - \boldsymbol{\bar{x}}^{(1)})=

\boldsymbol{x}^TS^{-1}\boldsymbol{x}-\boldsymbol{\bar{x}}^{(2)T}S^{-1}\boldsymbol{x} -\boldsymbol{x}^TS^{-1}\boldsymbol{\bar{x}}^{(2)}+\boldsymbol{\bar{x}}^{(2)T}S^{-1} \boldsymbol{\bar{x}}^{(2)}-\boldsymbol{x^T}S^{-1}\boldsymbol{x}+

\boldsymbol{\bar{x}}^{(1)T}S^{-1}\boldsymbol{x} + \boldsymbol{x^T}S^{-1}\boldsymbol{\bar{x}}^{(1)}-\boldsymbol{\bar{x}}^{(1)T}S^{-1}\boldsymbol{\bar{x}}^{(1)}=

-\boldsymbol{\bar{x}}^{(2)T}S^{-1}\boldsymbol{x} -\boldsymbol{x}^TS^{-1}\boldsymbol{\bar{x}}^{(2)}+\boldsymbol{\bar{x}}^{(2)T}S^{-1} \boldsymbol{\bar{x}}^{(2)}+\boldsymbol{\bar{x}}^{(1)T}S^{-1}\boldsymbol{x} + \boldsymbol{x^T}S^{-1}\boldsymbol{\bar{x}}^{(1)}-\boldsymbol{\bar{x}}^{(1)T}S^{-1}\boldsymbol{\bar{x}}^{(1)}=

2(\boldsymbol{\bar{x}}^{(1)T}-\boldsymbol{\bar{x}}^{(2)T})S^{-1}\boldsymbol{x} -(\boldsymbol{\bar{x}}^{(1)} -\boldsymbol{\bar{x}}^{(2)})^TS^{-1} (\boldsymbol{\bar{x}}^{(1)}+\boldsymbol{\bar{x}}^{(2)})=

2[S^{-1}(\boldsymbol{\bar{x}}^{(1)}-\boldsymbol{\bar{x}}^{(2)})]^T\boldsymbol{x} -(\boldsymbol{\bar{x}}^{(1)} -\boldsymbol{\bar{x}}^{(2)})^TS^{-1} (\boldsymbol{\bar{x}}^{(1)}+\boldsymbol{\bar{x}}^{(2)})

上述の計算は行列 $ S $ の逆行列 $ S^{-1} $ が転置と等しいという性質を使っています。 $ (S^{-1})^T = S^{-1} $ であるための条件は、行列 $ S $ が実対称行列(エルミート行列)であることです。

# Step 1: Generate random data for two groups
#np.random.seed(0)
mean1 = [-10, -10]
cov1 = [[1, 0.9], [0.9, 1]]
x1 = np.random.multivariate_normal(mean1, cov1, 10)

mean2 = [10, 10]
cov2 = [[1, 0.9], [0.9, 1]]
x2 = np.random.multivariate_normal(mean2, cov2, 10)
plt.scatter(x1[:, 0], x1[:, 1], c='b', label='x 1')
plt.scatter(x2[:, 0], x2[:, 1], c='g', label='x 2')
# 各群の平均を計算
mean_x1 = np.mean(x1, axis=0)
mean_x2 = np.mean(x2, axis=0)

# 各群の共分散を計算
cov_x1= np.cov(x1, rowvar=False)
cov_x2 = np.cov(x2, rowvar=False)

# 群内共分散行列（Sw）を計算
S= cov_x1*(len(x1)-1) + cov_x2*(len(x2)-1)
S = S / (len(x1) + len(x2) - 2)

w = np.linalg.inv(S).dot(mean_x1 - mean_x2)

# 射影軸の端点を計算
axis_start = np.min(np.vstack([x1.dot(w), x2.dot(w)])) * w / np.dot(w, w) 
axis_end = np.max(np.vstack([x1.dot(w), x2.dot(w)])) * w / np.dot(w, w) 
# 射影軸をプロット
plt.plot([axis_start[0], axis_end[0]], [axis_start[1], axis_end[1]]\
         , 'r-', label='Projection axis')
plt.xlabel('features 1')
plt.ylabel('features 2')

# Step 3: Generate a new sample data point
new_x = np.array([1, 1])
plt.scatter(new_x[0], new_x[1], c='y', label='x new')
plt.legend(bbox_to_anchor=(1, 1), loc='upper left')

# Step 4: Calculate Mahalanobis distance for both groups
D_2_2=(new_x-mean_x2).T.dot(np.linalg.inv(S)).dot(new_x-mean_x2)
D_1_2=(new_x-mean_x1).T.dot(np.linalg.inv(S)).dot(new_x-mean_x1)
# Step 5: Classify the new sample based on Mahalanobis distances
classification = "Group 1" if D_1_2 < D_2_2 else "Group 2"

f=w.T.dot(new_x)-1/2*((mean_x1-mean_x2).T).dot(np.linalg.inv(S).dot(mean_x1+mean_x2))
print('D_1^2-D_2^2:                  ',(D_2_2-D_1_2)/2)
print('classification:               ',classification)
print('fisher discriminant analysis: ',f)

• 一般の線形判別分析(LDA)の前提条件

  • 多変量正規性：独立変数は、グループ変数の各水準に対して正規分布に従います。
  • 分散/共分散の均質性（等分散性）：予測変数の水準を通じてグループ変数間の分散が同じです。
  • 多重共線性：予測変数間の相関が増加すると、予測力が低下する可能性があります。
  • 独立性：データはランダムにサンプリングされると仮定され、データは独立していると仮定されます。

• フィッシャーの線形判別(FDA)の前提条件

  • 多重共線性：予測変数間の相関が増加すると、予測力が低下する可能性があります。
  • 独立性：データはランダムにサンプリングされると仮定され、データは独立していると仮定されます。
    多変量正規性と等分散性はFDAでは仮定されていません。

• 生起確率の相違

群$G_1$,$G_2$に対応する事象の生起確率が異なるときは,各群に属するデータが生じ
る事前確率$\pi_1=1-\pi_2$を導入して$f(\boldsymbol{x})-\log(\pi_2/\pi_1)$の正負を利用して判別を行うことができます。

ここでも$-\log(\pi_2/\pi_1)$はバイアス項です。
条件付き確率 $P(G_1|\boldsymbol{x})$ と $P(G_2|\boldsymbol{x})$は、新しいデータ点 $x$ が与えられたときにそれがクラス $G_1$ または $G_2$ に属する確率を表しています。ベイズの定理を用いると、これらの条件付き確率は次のように表されます：

P(G_1|\boldsymbol{x}) = \frac{P(\boldsymbol{x}|G_1) \pi_1}{P(\boldsymbol{x})}

P(G_2|\boldsymbol{x}) = \frac{P(\boldsymbol{x}|G_2) \pi_2}{P(\boldsymbol{x})}

ここで $P(\boldsymbol{x}|G_1)$ と $P(\boldsymbol{x}|G_2)$ は射影されたデータの密度関数、$\pi_1$ と $\pi_2$ はクラス $G_1$ と $G_2$ の事前確率、$P(\boldsymbol{x})$ は全体のデータの密度関数です。

新しいデータ点 $\boldsymbol{x}$ のクラスは、$P(G_1|\boldsymbol{x})$ と $P(G_2|\boldsymbol{x})$ の比較によって決定されます。特に、$P(G_1|\boldsymbol{x}) > P(G_2|\boldsymbol{x})$ ならば $\boldsymbol{x}$ は $G_1$ に、そうでなければ $G_2$ に分類されます。
この条件は以下のように書き換えられます：

\frac{P(\boldsymbol{x}|G_1) \pi_1}{P(\boldsymbol{x}|G_2) \pi_2} > 1

この不等式に自然対数を適用すると、

\log(P(\boldsymbol{x}|G_1)) + \log(\pi_1) > \log(P(\boldsymbol{x}|G_2)) + \log(\pi_2)

が得られ、これが判別に用いる関数となります。この式が、事前確率を考慮した判別関数の形になります。ここで登場する $\log(\pi_1)$ と $\log(\pi_2)$ は、それぞれのクラスの事前確率の対数を取ったもので、判別関数に直接的な影響を与えます。

訓練データでのクラスの比率と新しいサンプルでのクラスの比率が異なる場合、特に注意が必要です。以下はそのような状況を考慮するためのいくつかのアプローチです。

クラスの不均衡が訓練データに存在する場合

1. リサンプリング: 少数派クラスをオーバーサンプリングするか、多数派クラスをアンダーサンプリングして、クラス間の不均衡を解消する。
2. コスト感受性のある学習: クラス間の不均衡を明示的に考慮するコスト関数を使用する。

新しいサンプルでのクラス比率が異なる場合

1. 事前確率の調整: モデルが新しいサンプルでより一般的に使われるクラス比率に合わせて調整されるよう、事前確率 $ \pi_1 $ と $ \pi_2 $ を設定します。
2. しきい値の調整: 分類のしきい値を調整して、少数派クラスが適切に識別されるようにします。

これらのアプローチは、訓練フェーズとテスト（または運用）フェーズでのクラス比率の不一致を考慮に入れるためのものです。特に、新しいサンプルでのクラス比率が訓練データと異なると予想される場合、事前確率の調整やしきい値の調整が有用です。これにより、モデルが新しいサンプルに対してより頑健な性能を発揮できます。

このようなバイアス項は次のような注意点もあります。射影された $ y $ の分布とその分散がどれだけ大きいかによって、事前確率 $ \pi_1 $ と $ \pi_2 $（およびそれから計算される $ \log\left(\frac{\pi_2}{\pi_1}\right) $）の影響が大きく変わります。

• もし $ y $ の分散が非常に大きければ、それが持つ情報が強く、事前確率の影響は相対的に小さくなる可能性が高いです。
• 逆に $ y $ の分散が非常に小さいと、事前確率が持つ影響が大きくなる可能性があります。

これは、$ y $ の分布がどれだけ「分離」されているか、すなわち各クラスの $ y $ 値がどれだけ重なっているかにも依存します。各クラスの $ y $ 値が明確に分離されている場合、事前確率の影響はそれほど大きくありません。しかし、もし $ y $ 値が大きく重なっている場合、事前確率が重要な役割を果たすことになります。

最終的には、これらの要素は全て、具体的なデータと問題設定に依存します。したがって、モデルの性能を評価し、必要に応じて事前確率を調整することが重要です。

from scipy.stats import multivariate_normal

# 1. データの生成
#np.random.seed(0)
mean1 = [0, 0]
cov1 = [[1, 0.5], [0.5, 1]]
x1 = np.random.multivariate_normal(mean1, cov1, 100)

mean2 = [2, 2]
cov2 = [[1, 0.5], [0.5, 1]]
x2 = np.random.multivariate_normal(mean2, cov2, 100)

# 2. 事前確率の設定
pi_1 = 0.5  # 事前確率 for クラス G_1
pi_2 = 1 - pi_1  # 事前確率 for クラス G_2

# 3. 新しいサンプルデータ点の生成
new_sample = np.array([1, 1])

# 4. P(x|G_1) と P(x|G_2) の計算
P_x_given_G1 = multivariate_normal.pdf(new_sample, mean=mean1, cov=cov1)
P_x_given_G2 = multivariate_normal.pdf(new_sample, mean=mean2, cov=cov2)

# 5. 判別関数の計算
LHS = np.log(P_x_given_G1) + np.log(pi_1)
RHS = np.log(P_x_given_G2) + np.log(pi_2)

# 6. 判別結果の判定
if LHS > RHS:
    classification = "G_1"
else:
    classification = "G_2"

print("Classification result for the new data point x is:", classification)

Classification result for the new data point x is: G_2

• 分散共分散の不均一性

共通の分散共分散行列ではなく,各群の分散構造が異なることを仮定し,$D_j^2=(\boldsymbol{x}-\boldsymbol{\bar{x}}^{(j)})^TS^{-1}(\boldsymbol{x}-\boldsymbol{\bar{x}}^{(j)})$を用いて2次式

$q(\boldsymbol{x})=D_2^2-D_1^2=
(\boldsymbol{x}-\boldsymbol{\bar{x}}^{(2)})^TS^{-1}(\boldsymbol{x}-\boldsymbol{\bar{x}}^{(2)})-
(\boldsymbol{x}-\boldsymbol{\bar{x}}^{(1)})^TS^{-1}(\boldsymbol{x}-\boldsymbol{\bar{x}}^{(1)})$

の正負で判別を行うことを,2次判別分析とよび,関数$q(\boldsymbol{x})$を2次判別関数とよぶ。1次関数による線形判別分析と比較して柔軟な判別境界を表現することができる一方で,　各クラスで異なる分散共分散行列を推定するため,サンプルサイズが小さい場合には数値的な不安定性が間題となったり,データへの過適合が生じる可能性が高くなることに注意が必要です.

# 1. データの生成
mean1 = [0, 0]
cov1 = [[1, 0.5], [0.5, 1]]
x1 = np.random.multivariate_normal(mean1, cov1, 100)

mean2 = [2, 2]
cov2 = [[1.5, 0.2], [0.2, 1.5]]
x2 = np.random.multivariate_normal(mean2, cov2, 100)

# 2. 各クラスの平均と共分散行列を計算
mean1_sample = np.mean(x1, axis=0)
mean2_sample = np.mean(x2, axis=0)

cov1_sample= np.cov(x1, rowvar=False)
cov2_sample = np.cov(x2, rowvar=False)

# 3. 新しいサンプルデータ点の生成
new_sample = np.array([1, 1])

# 4. 各クラスに対するマハラノビス距離を計算
D1_square = (new_sample - mean1_sample).T.dot(np.linalg.inv(cov1_sample)).dot(new_sample - mean1_sample)
D2_square = (new_sample - mean2_sample).T.dot(np.linalg.inv(cov2_sample)).dot(new_sample - mean2_sample)

q_x = D2_square - D1_square

# 5. 判別結果の判定
classification_qda = "G_1" if q_x < 0 else "G_2"

print(classification_qda, q_x)

# Generating a grid of values for plotting decision boundary
x_vals = np.linspace(-3, 5, 200)
y_vals = np.linspace(-3, 5, 200)
xx, yy = np.meshgrid(x_vals, y_vals)
grid_points = np.c_[xx.ravel(), yy.ravel()]

# Calculating decision function for each point in the grid
q_values = []
for point in grid_points:
    q = (point - mean2).T @ np.linalg.inv(cov2) @ (point - mean2) - (point - mean1).T @ np.linalg.inv(cov1) @ (point - mean1)
    q_values.append(q)
q_values = np.array(q_values).reshape(xx.shape)

# Plotting
plt.figure(figsize=(5, 4))
plt.scatter(x1[:, 0], x1[:, 1], c='blue', s=50, edgecolor='k', label='G_1')
plt.scatter(x2[:, 0], x2[:, 1], c='red', s=50, edgecolor='k', label='G_2')
plt.scatter(1, 1, c='green', s=100, marker='x', label='New sample x')
plt.contour(xx, yy, q_values, levels=[0], colors='k')  # Decision boundary where q(x) = 0
plt.title("Quadratic Discriminant Analysis")
plt.xlabel("x1")
plt.ylabel("x2")
plt.legend()
plt.grid(True)
plt.show()

• マハラノビス平方距離(wiki 英語版)

$\mathbb{R}^{N}$ 上の確率分布 $Q$ が与えられたとき、その平均が

\boldsymbol{\mu}=(\mu _{1},\dots ,\mu _{N})^{\mathsf {T}}

であり、正定値の共分散行列が $S$ である場合、点 ${\boldsymbol {x}}=(x_{1},\dots ,x_{N})^{\mathsf {T}}$ と $Q$ との間のマハラノビス距離は
$$d_{M}({\boldsymbol {x}},Q)={\sqrt {({\boldsymbol {x}}-\boldsymbol{ \bar{x}}})^{\mathsf {T}}S^{-1}(\boldsymbol{ {x}}-\boldsymbol {\bar{x}})}$$
となります。ここで全サンプルの平均べクトルは

$$\bar{\boldsymbol{x}}=\left(\sum_{j=1}^gn_j\right)^{-1}\sum_{j=1}^g\sum_{i=1}^{n_j}\boldsymbol{x}_i^{(j)}$$

です。実際には、分布 $Q$ は通常、未知の基本分布から独立同一に得られたサンプル分布です。$S$ はサンプルの共分散行列です。

$\mathbb{R}^{N}$ 内の2つの点 ${\boldsymbol {x}}$ および ${\boldsymbol {y}}$ に対して、それらの間の $Q$ に関するマハラノビス距離は
$$d_{M}({\boldsymbol {x}},{\boldsymbol {y}};Q)={\sqrt {({\boldsymbol {x}}-{\boldsymbol {y}})^{\mathsf {T}}S^{-1}({\boldsymbol {x}}-{\boldsymbol {y}})}}$$
となり、これは
$$d_{M}({\boldsymbol {x}},Q)=d_{M}({\boldsymbol {x}},{\boldsymbol {\mu }};Q)$$
を意味します。
$S$ は正定値であるため、$S^{-1}$ も正定値であり、そのため平方根は常に定義されます。

# 1. データ生成
np.random.seed(42)
mean1 = [0, 0]
cov1 = [[2, 1], [1, 2]]
x1 = np.random.multivariate_normal(mean1, cov1, 100)

mean2 = [10, 10]
cov2 = [[4, -1], [-1, 4]]
x2 = np.random.multivariate_normal(mean2, cov2, 100)

mean3 = [0, 10]
cov3 = [[1, 0], [0, 1]]
x3 = np.random.multivariate_normal(mean3, cov3, 100)
# 2. マハラノビス平方距離の計算 (スクラッチ)
def mahalanobis_scratch(x, mean, cov):
    diff = x - mean
    inv_cov = np.linalg.inv(cov)
    return diff.T @ inv_cov @ diff

# 3群のデータに対してマハラノビス平方距離を計算
distances1_scratch = [mahalanobis_scratch(xi, mean1, cov1) for xi in x1]
distances2_scratch = [mahalanobis_scratch(xi, mean2, cov2) for xi in x2]
distances3_scratch = [mahalanobis_scratch(xi, mean3, cov3) for xi in x3]

# 可視化


# 散布図は2次元データなので、マハラノビス平方距離を色として表現します
plt.figure(figsize=(4, 3))
plt.scatter(x1[:, 0], x1[:, 1], c=distances1_scratch, cmap='viridis')
plt.scatter(x2[:, 0], x2[:, 1], c=distances2_scratch, cmap='viridis')
plt.scatter(x3[:, 0], x3[:, 1], c=distances3_scratch, cmap='viridis')
plt.colorbar(label='Mahalanobis Distance')
plt.title("Scatter plot with Mahalanobis Distances as color")
plt.xlabel("Feature 1")
plt.ylabel("Feature 2")
plt.grid(True)
plt.show()

• 正準判別分析

サンプルが属する群が2より多い場合の多群の判別においても,各群の平均(重心)とサンプル$x$とのマハラノビス平方距離を求めることができます。そして,それが最小となる群にサンプル$x$を分類することで,多群の判別が可能となります。

$G$群の$p$次元データを判別性が高くなるように低次元に射影する方法として正準判別分析を用います.群$G_j,j=1,\cdots,g$
に属するサンプル

\{x_i^{(j)}\}_{i=1}^{n_j}, j=1,\cdots,g

の平均べクトルを$\bar{\boldsymbol{x}}^{(j)}=\frac{1}{n_j}\sum_{i=1}^{n_j}\boldsymbol{x}^{(j)}$,
分散共分散行列を

$S_j=\frac{1}{n_j-1}\sum_{i=1}^{n_j}(\boldsymbol{\bar{x}_i}^{(j)}-\boldsymbol{\bar{x}}^{(j)})(\boldsymbol{x}_i^{(j)}-\boldsymbol{\bar{x}}^{(j)})^T$

$\bar{y}=\boldsymbol{w}^T\bar{\boldsymbol{x}}$とします.群間変動行列と群内変動行列をそれぞれ

$S_B=\sum_{j=1}^{g}n_j(\boldsymbol{\bar{x}}^{(j)}-\boldsymbol{\bar{x}})(\bar{\boldsymbol{x}}^{(j)}-\boldsymbol{\bar{x}})^T$

および

$S_W=\sum_{j=1}^{g}(n_j-1)S_j$

として,フィッシャーの判別分析と同様に群間,群内の分散比を最大にする射影べクトル$\boldsymbol{w} \in \mathbb{R}^p$を求めれば射影軸が得られます。

フィッシャーの線形判別分析は、主に2つのクラス間の分離を最大化する射影軸を特定することに焦点を当てています。一方、正準判別分析は、3つ以上の複数のクラスが存在する状況にも対応しており、複数の射影軸を同時に特定することができます。このプロセスでは、群間の変動（クラス間の差）と群内の変動（クラス内の差）の比率を最大化する方向にデータを射影することで、マハラノビス距離の二乗を最大化します。

固有値分解を用いてこの比率を最大化する射影ベクトル（固有ベクトル）を特定します。たとえば、
$(\boldsymbol{w}_1^T\boldsymbol{x}, \boldsymbol{w}_2^T\boldsymbol{x})$ のようにデータを新しい空間に射影し、この射影された空間でサンプルをプロットすることによって、各クラスに属するサンプルの分布を視覚的に分析することができます。また、射影されたデータを新しい特徴量として利用し、さらにクラスタ分析や他の機械学習タスクに活用することも可能です。

# 正準判別分析

# 全体の平均を計算
overall_mean = np.mean(np.vstack((x1, x2, x3)), axis=0)

# 群間変動行列 S_B を計算
S_B1 = len(x1) * np.outer((mean1 - overall_mean), (mean1 - overall_mean))
S_B2 = len(x2) * np.outer((mean2 - overall_mean), (mean2 - overall_mean))
S_B3 = len(x3) * np.outer((mean3 - overall_mean), (mean3 - overall_mean))
S_B = S_B1 + S_B2 + S_B3

# 群内変動行列 S_W を計算
S_W = np.cov(x1, rowvar=False) * (len(x1) - 1) + \
    np.cov(x2, rowvar=False) * (len(x2) - 1) + \
    np.cov(x3, rowvar=False) * (len(x3) - 1)

# 固有値と固有ベクトルを計算
eig_vals, eig_vecs = np.linalg.eig(np.linalg.inv(S_W).dot(S_B))

# 固有値を降順にソート
sorted_indices = np.argsort(eig_vals)[::-1]
eig_vals_sorted = eig_vals[sorted_indices]
eig_vecs_sorted = eig_vecs[:, sorted_indices]

# 最大の固有値に対応する2つの固有ベクトルを取得
w1 = eig_vecs_sorted[:, 0]
w2 = eig_vecs_sorted[:, 1]

# データを射影
x1_projected = np.vstack((x1 @ w1, x1 @ w2)).T
x2_projected = np.vstack((x2 @ w1, x2 @ w2)).T
x3_projected = np.vstack((x3 @ w1, x3 @ w2)).T

# プロット
plt.figure(figsize=(4, 3))
plt.scatter(x1_projected[:, 0], x1_projected[:, 1], label="Group 1", alpha=0.6)
plt.scatter(x2_projected[:, 0], x2_projected[:, 1], label="Group 2", alpha=0.6)
plt.scatter(x3_projected[:, 0], x3_projected[:, 1], label="Group 3", alpha=0.6)
plt.title("Canonical Discriminant Analysis")
plt.xlabel("Canonical Axis 1")
plt.ylabel("Canonical Axis 2")
plt.legend()
plt.grid(True)
plt.show()

• 白色化変換

白色化変換は、データの共分散行列が単位行列になるようにデータを変換する操作です。具体的には、元のデータの共分散行列の平方根の逆行列を用いてデータを変換します。この変換により、変換後のデータのすべての成分は無相関となり、それぞれの分散が1になります。マハラノビス距離はデータの共分散構造を考慮して距離を測定しますが、白色化変換を適用することで、この共分散構造を取り除き、変換後のデータにおける距離を単純なユークリッド距離として評価することができます。

def whiten(X, mean, cov):
    # 1. Center the data
    X_centered = X - mean

    # 2. Get the covariance matrix is already given as 'cov'

    # 3. Eigenvalue decomposition
    eig_vals, eig_vecs = np.linalg.eigh(cov)

    # 4. Create a diagonal matrix with the inverse square root of eigenvalues
    diag_inv_sqrt = np.diag(1.0 / np.sqrt(eig_vals))

    # 5. Compute the whitening matrix
    whitening_matrix = np.dot(eig_vecs, np.dot(diag_inv_sqrt, eig_vecs.T))

    # 6. Whiten the data
    X_whitened = np.dot(X_centered, whitening_matrix)

    return X_whitened

# Whiten each group of data
x1_whitened = whiten(x1, mean1, cov1)
x2_whitened = whiten(x2, mean2, cov2)
x3_whitened = whiten(x3, mean3, cov3)

# After whitening, compute the means and covariances again
mean1_w = np.mean(x1_whitened, axis=0)
mean2_w = np.mean(x2_whitened, axis=0)
mean3_w = np.mean(x3_whitened, axis=0)
cov1_w = np.cov(x1_whitened, rowvar=False)
cov2_w = np.cov(x2_whitened, rowvar=False)
cov3_w = np.cov(x3_whitened, rowvar=False)

# S_W and S_B for whitened data
S_W_w = cov1_w * len(x1_whitened) + cov2_w * len(x2_whitened) + cov3_w * len(x3_whitened)
S_B1_w = len(x1_whitened) * np.outer((mean1_w - overall_mean), (mean1_w - overall_mean))
S_B2_w = len(x2_whitened) * np.outer((mean2_w - overall_mean), (mean2_w - overall_mean))
S_B3_w = len(x3_whitened) * np.outer((mean3_w - overall_mean), (mean3_w - overall_mean))
S_B_w = S_B1_w + S_B2_w + S_B3_w

# Eigenvalue decomposition for whitened data
eig_vals_w, eig_vecs_w = np.linalg.eig(np.linalg.inv(S_W_w).dot(S_B_w))
sorted_indices_w = np.argsort(eig_vals_w)[::-1]
eig_vals_sorted_w = eig_vals_w[sorted_indices_w]
eig_vecs_sorted_w = eig_vecs_w[:, sorted_indices_w]

# Projection vectors for whitened data
w1_w = eig_vecs_sorted_w[:, 0]
w2_w = eig_vecs_sorted_w[:, 1]

# Project the whitened data
x1_projected_w = np.vstack((x1_whitened @ w1_w, x1_whitened @ w2_w)).T
x2_projected_w = np.vstack((x2_whitened @ w1_w, x2_whitened @ w2_w)).T
x3_projected_w = np.vstack((x3_whitened @ w1_w, x3_whitened @ w2_w)).T

# 1次元データに対してマハラノビス平方距離を計算する関数の修正

def compute_mahalanobis_1D(data, mean, variance):
    distances = []
    for xi in data:
        diff = xi - mean
        dist = (diff**2) / variance
        distances.append(dist)
    return distances

variance1_projected = np.var(x1_projected)
variance2_projected = np.var(x2_projected)
variance3_projected = np.var(x3_projected)

distances1_projected = compute_mahalanobis_1D(x1_projected, np.mean(x1_projected), variance1_projected)
distances2_projected = compute_mahalanobis_1D(x2_projected, np.mean(x2_projected), variance2_projected)
distances3_projected = compute_mahalanobis_1D(x3_projected, np.mean(x3_projected), variance3_projected)

distances1_projected[:5], distances2_projected[:5], distances3_projected[:5]  # それぞれの群の先頭5つのマハラノビス平方距離

# 正準変数を計算
C_x1 = x1 @ eig_vecs_sorted
C_x2 = x2 @ eig_vecs_sorted
C_x3 = x3 @ eig_vecs_sorted

# 2Dプロット
plt.figure(figsize=(5, 4))
plt.scatter(C_x1[:, 0], C_x1[:, 1], label="Group 1", alpha=0.6)
plt.scatter(C_x2[:, 0], C_x2[:, 1], label="Group 2", alpha=0.6)
plt.scatter(C_x3[:, 0], C_x3[:, 1], label="Group 3", alpha=0.6)
plt.title("Canonical Discriminant Analysis - 2D View")
plt.xlabel("Canonical Axis 1")
plt.ylabel("Canonical Axis 2")
plt.legend()
plt.grid(True)
plt.show()