「入門はじめての多変量解析」をPythonで写経 Vol.24 ～ ５章「はじめての判別分析」③マハラノビス距離による判別
５章「はじめての判別分析」
書籍の著者 石村貞夫 先生、石村光資郎 先生
書籍「入門はじめての多変量解析」５章「はじめての判別分析」の Python写経活動記録 です。
多変量解析の入門を Python と一緒に学ぶ写経シリーズです。
判別分析はいわゆる「分類」タスクです。
書籍は次の２種類の判別分析を掲載しています。
1️⃣ 線形判別分析（２変数、２群）
2️⃣ マハラノビス距離による判別 ◀ 今回はこちら！
この記事は、マハラノビス距離とは何者か をざっくり概観して、 マハラノビス距離で境界線を引く判別方法 を学びます。
Python ライブラリの活用も進めます！
ChatGPT 活用型学習で進めてまいります！
では書籍を開いて多変量解析の旅に出かけましょう🚀
はじめに
このブログシリーズは、書籍「入門はじめての多変量解析」（東京図書、「テキスト」と呼びます）の Python 写経を通じて得た「多変量解析の楽しさ」をご紹介します。
書籍の紹介と引用表記はリンク先の記事に掲載しています。
５章 はじめての判別分析
この記事は５章の以下のSectionを取り扱います。
5.7 １変数のマハラノビスの距離
5.8 ２変数のマハラノビスの距離
5.9 マハラノビスの距離による境界線
記事に用いるデータは、テキストに掲載されたデータそのものを引用しています。
データ件数の少ないものはコード上でデータを登録し、データ件数の多いものはCSVファイル化してデータを読み込みしています。
この記事で用いるライブラリをインポートします。
### インポート
# 数値計算
import numpy as np
import pandas as pd
import sympy
# マハラノビスの距離
from scipy.spatial import distance
from sklearn.covariance import EmpiricalCovariance
# 統計
import scipy.stats as stats
# LaTeX表示
from IPython.display import Math
# 描画
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo' # または import japanize_matplotlib
マハラノビス距離のイメージ
マハラノビス距離は「距離の一種」です。
データの「分散・共分散」を考慮した距離の概念です。
次の等高線は、２変数、２グループのマハラノビス距離を示しています。
赤と青の２つのグループの中心（バツ点）から「楕円」の等高線が引かれています。
１本の等高線は中心から等しい距離にあります。
この距離がマハラノビス距離なのです。
例えば赤いグループの 0.1 の楕円は、中心から 0.1 の距離にあります。
マーカーAが 2.7 あたり、マーカーB が 3.5 あたりの青い点に注目しましょう。
青いグループの中心からのマハラノビス距離は 3.0 と 4.0 の間です。
一方、赤いグループの中心からのマハラノビス距離は 1.0 と 2.0 の間です。
青いグループに属するデータ点ですが、マハラノビス距離ではかると、赤いグループのほうが近いのです！
日常的に使用する距離は、中心からコンパスで描くような「同心円」が等高線になります。
このような距離をユークリッド距離と呼ぶこともあります。
一方で、マハラノビス距離は分散・共分散を用いることで、データのバラツキの方向を考慮して、「楕円」の等高線になっています。
そしてそして…
２つのグループの中心からの距離が等しいところが「境界線」の候補です！
では、１変数のマハラノビス距離から進めていきましょう！
１変数のマハラノビス距離 p.204～
■ データの準備
テキスト p.204 表 5.7.1 のメタボ検査結果データをお借りします。
### 検査結果データ p.204 表5.7.1 ※体型 0:標準体型、1:肥満体型
data2 = pd.DataFrame(
{'被験者No.': range(1, 21),
'検査結果': [22, 20, 23, 23, 17, 24, 23, 18, 22, 19,
40, 35, 27, 22, 25, 26, 19, 31, 30, 36],
'体型': np.hstack([np.zeros(10), np.ones(10)]).astype(int)})
data2【実行結果】
体型＝０が標準体型のグループ１、体型＝１が肥満体型のグループ２です。
両グループの標本サイズは 10 です。
このデータを教師データにして、次の R さんの検査結果からどちらのグループに属するかを検討します。
# Rさんの測定値
R = 25
print('Rさんの検査結果:', 25)【実行結果】
後の処理で利用するグループごとの平均・分散・標準偏差を計算します。
分散・標準偏差の分母は $${N-1}$$ です。
## グループごとの統計量の算出
# グループごとの平均値
mean_g1, mean_g2 = data2.groupby(['体型'])['検査結果'].mean().values
print(f'平均 G1={mean_g1:7.4f}, G2={mean_g2:7.4f}')
# グループごとの分散
var_g1, var_g2 = data2.groupby(['体型'])['検査結果'].var().values
print(f'分散 G1={var_g1:7.4f}, G2={var_g2:7.4f}')
# グループごとの標準偏差
std_g1, std_g2 = data2.groupby(['体型'])['検査結果'].std(ddof=1).values
print(f'標準偏差 G1={std_g1:7.4f}, G2={std_g2:7.4f}')【実行結果】
グループ２の方がバラツキが大きいです。
■ データをながめる
可視化して２つのグループの特徴を確認しましょう。
こちらは１次元のプロットです。
テキスト p.204 図 5.7.1 ライクなチャートです。
### p.204 ２つのグループのグラフ表現 図5.7.1 その１ １次元プロット版
# データ点の描画
sns.swarmplot(x=data2['検査結果'],
hue=np.where(data2['体型'] == 1, 'グループ$G_2$', 'グループ$G_1$'),
palette=['tab:blue', 'tomato'], s=10, alpha=0.9)
# Rさんのデータ点の描画
plt.plot(R, 0.1, '*', color='tab:green', ms=15, label='R')
# グループ１平均の垂直線の描画
plt.axvline(mean_g1, color='tab:blue', lw=1, ls='--', label='$G_1$平均')
# グループ２平均の垂直線の描画
plt.axvline(mean_g2, color='tomato', lw=1, ls='--', label='$G_2$平均')
# 修飾
plt.ylim(-0.5, 0.5)
plt.legend(bbox_to_anchor=(1.3, 1));【実行結果】
Rさん（星印）はグループ１の平均とグループ２の平均のほぼ真ん中に位置しているようです。
テキストによると：
グループ $${G_1}$$ の平均との差 $${= 25 - 21.1 = 3.9}$$
グループ $${G_2}$$ の平均との差 $${= 29.1 - 25 = 4.1}$$
Rさんはグループ１の平均の方が近いです。
ただ、グループ１はバラツキが小さく、グループ２はバラツキが大きいです。
平均との差の比較がいいのかどうか、迷います。
続いてヒストグラム形式です。
### p.204 ２つのグループのグラフ表現 図5.7.1 その２ ヒストグラム版
# bins（階級）の設定
bins = np.arange(data2['検査結果'].min()-0.5, data2['検査結果'].max()+1.5)
# ヒストグラムの描画
sns.histplot(x=data2['検査結果'], multiple='stack', ec='white', bins=bins,
hue=np.where(data2['体型']==1, 'グループ$G_2$', 'グループ$G_1$'),
palette=['tab:blue', 'tomato'], alpha=0.7)
# Rさんのデータ点の描画
sns.scatterplot(x=[R], y=[1.2], marker='*', s=500, color='tab:green')
# グループ１平均の垂直線の描画
plt.axvline(mean_g1, color='tab:blue', lw=1, ls='--')
# グループ２平均の垂直線の描画
plt.axvline(mean_g2, color='tomato', lw=1, ls='--')
# Rさんの注釈の表示
plt.annotate(xy=[25.4, 1.25], xytext=[26, 1.6], text='Rさん', fontsize=14,
arrowprops=dict(facecolor='black', arrowstyle='->'))
# 修飾
plt.yticks([0, 1, 2, 3]);【実行結果】
Rさんはグループ１に近いような気もしますし、グループ２に含まれているような気もします。
テキスト p.205 図 5.7.2 の分布図で２つのグループのバラツキの違いを確認しましょう。
２つのグループが正規分布に従うと仮定して、データの平均・標準偏差をパラメータとする正規分布の確率密度関数を描きます。
### バラツキの異なる分布 p.205 図5.7.2
## 描画データ: 正規分布の確率密度関数の算出
# グループG1の分布
x_vals1 = np.linspace(10, 30, 101)
y_vals1 = stats.norm.pdf(x_vals1, loc=mean_g1, scale=std_g1)
# グループG2の分布
x_vals2 = np.linspace(10, 50, 101)
y_vals2 = stats.norm.pdf(x_vals2, loc=mean_g2, scale=std_g2)
## 描画
# 描画領域の設定
plt.figure(figsize=(8, 4))
# グループG1の正規分布の確率密度関数の描画
plt.plot(x_vals1, y_vals1, color='tab:blue', label='グループ$G_1$の分布')
# グループG1の平均の垂直点線の描画
plt.vlines(mean_g1, 0, stats.norm.pdf(mean_g1, loc=mean_g1, scale=std_g1),
color='black', ls='--', lw=0.5)
# グループG2の正規分布の確率密度関数の描画
plt.plot(x_vals2, y_vals2, color='tomato', label='グループ$G_2$の分布')
# グループG2の平均の垂直点線の描画
plt.vlines(mean_g2, 0, stats.norm.pdf(mean_g2, loc=mean_g2, scale=std_g2),
color='black', ls='--', lw=0.5)
# グループG1の1σの垂直線の描画
plt.axvline(mean_g1-std_g1, color='tab:blue', lw=1, ls='--',
label=f'$G_1$ 平均$\pm1\sigma({std_g1:.2f})$')
plt.axvline(mean_g1+std_g1, color='tab:blue', lw=1, ls='--')
# グループG2の1σの垂直線の描画
plt.axvline(mean_g2-std_g2, color='tomato', lw=1, ls='--',
label=f'$G_2$ 平均$\pm1\sigma({std_g2:.2f})$')
plt.axvline(mean_g2+std_g2, color='tomato', lw=1, ls='--')
# Rさんの25の垂直線の描画
plt.axvline(R, color='tab:green', ls='--', label='R')
# 修飾
plt.ylim(0, 0.18)
plt.xticks([mean_g1, 25, mean_g2])
plt.legend()
# 不要な軸・軸目盛りの削除
plt.yticks([])
plt.gca().spines['right'].set_visible(False)
plt.gca().spines['top'].set_visible(False)
plt.gca().spines['left'].set_visible(False);【実行結果】
グループ１の分布はバラツキが小さく、グループ２はバラツキが大きいことがよく分かります。
青と赤の点線は各グループの平均から $${1\sigma}$$（１標準偏差）離れている位置です。
Rさん（緑線）はグループ１の $${\pm 1 \sigma}$$ には含まれていませんが、グループ２の $${\pm 1 \sigma}$$ には含まれています。
ますます、混迷を極めてきました…
ひとまず、データの KDE プロットを鑑賞して落ち着きましょう。
### 参考: KDEプロット
sns.kdeplot(x=data2['検査結果'],
hue=np.where(data2['体型']==1, 'グループ$G_2$', 'グループ$G_1$'))
plt.axvline(R, color='tab:green', ls='--', label='R')
plt.text(R+0.5, 0.07, '← R 25', fontsize=12);【実行結果】
ここからマハラノビス距離にどんどん近づいていきます。
■ データの標準化で位置をはかる
テキストはデータの標準化を利用して「グループ１の平均とグループ２の平均から等距離の位置」を探ります。
📊 データの標準化
$$
データ \longmapsto \cfrac{データ - 平均}{標準偏差}
$$
グループ１の平均・標準偏差で標準化結果 と、グループ２の平均・標準偏差で標準化結果が等しくなる$${x}$$ を探ります。
$$
\cfrac{x-G_1の平均}{G_1の標準偏差} = \cfrac{G_2の平均-x}{G_2の標準偏差} \\
$$
Pythonで方程式を解きましょう。
### データの標準化を用いて各グループの平均から等距離の位置を求める p.206
# sympyの変数定義
x = sympy.symbols('x')
# 方程式を解く
solved = sympy.solve([(x - mean_g2) / std_g2 - (mean_g1 - x) / std_g1])
# 結果の表示
print(f'x = {float(solved[x]):.4f}')【実行結果】
テキストの $${23.3}$$ と合っています。
$${23.3}$$ が標準化したデータの視点でグループ１の平均とグループ２の平均から等距離にあることが分かりました。
R さんの $${25}$$ は $${23.3}$$ より大きいので、右側のグループ２に属すると推測されます。
可視化して確かめます。
テキスト p.207 図 5.7.3 に相当します。
### バラツキの異なる分布 p.205 図5.7.2
## 描画
# 描画領域の設定
plt.figure(figsize=(7, 3))
# グループG1の正規分布の確率密度関数の描画
plt.plot(x_vals1, y_vals1, color='tab:blue', label='グループ$G_1$の分布')
# グループG2の正規分布の確率密度関数の描画
plt.plot(x_vals2, y_vals2, color='tomato', label='グループ$G_2$の分布')
# グループ1とグループ2の境界線の描画
plt.axvline(x, color='tab:red', ls='--', lw=1)
# 両グループの平均から等距離点の描画
plt.plot([x], [0], 'o', color='tab:red', label='中心')
# Rさんのデータ点の描画
plt.plot([R], [0], 'D', color='tab:green', label='R')
# Rさんの注釈の表示
plt.annotate(xy=[25.45, 0.002], xytext=[27, 0.03], text='Rさん:25', fontsize=14,
arrowprops=dict(facecolor='tab:red', arrowstyle='->'))
# 修飾
plt.ylim(-0.01, 0.18)
plt.xticks([round(x, 1)])
plt.legend()
plt.gca().spines['bottom'].set_position('zero')
# 不要な軸・軸目盛りの削除
plt.yticks([])
plt.gca().spines['right'].set_visible(False)
plt.gca().spines['top'].set_visible(False)
plt.gca().spines['left'].set_visible(False);【実行結果】
R さんが $${23.3}$$ より大きく、グループ２に近いことが分かりました。
■ １変数のマハラノビス距離
ここまでは各グループの平均から等距離の位置を起点にして R さんの所属グループを見つけに行きました。
ここからは、各グループの平均からの距離を比べて R さんの所属グループを検討します。
各グループの平均からの距離 $${D}$$ も「データの標準化」で求めます。
R さんの検査結果 $${25}$$ について、グループ１の平均からの距離を $${D_1}$$、グループ２の平均から距離を $${D_2}$$ として「ある距離」を求めます。
$$
D_1 = \cfrac{|R - \mu_{G_1}|}{\sigma_{G_1}} = \cfrac{|25- 21.1|}{2.42} = 1.61 \\
\\
D_2 = \cfrac{|R - \mu_{G_2}|}{\sigma_{G_2}} = \cfrac{|25- 29.1|}{6.57} = 0.62
$$
記号 $${|\cdot|}$$ は絶対値です。
$${D_2}$$ の方が小さいので、R さんの検査結果はグループ２の平均の方が近い、したがって、Rさんはグループ２に所属するということです。
Python で計算しましょう。
### 標準化を取り入れた距離 p.207
# 標準化データを用いた距離の算出
D1 = abs(R - mean_g1) / std_g1
D2 = abs(R - mean_g2) / std_g2
# 結果の表示
print(f'D₁:{D1:.3f}, D₂:{D2:.3f}')【実行結果】
この距離 $${D_1, D_2}$$ をマハラノビス距離と呼んだり、この距離の二乗をマハラノビス距離と呼んだりしています。
◆
📊 １変数のマハラノビス距離（テキストの定義）
テキストは「二乗」をマハラノビス距離としています。
$$
D^2 = \cfrac{(x-平均)^2}{分散}=\cfrac{(x-\bar{x})^2}{s^2}
$$
$${\bar{x}}$$ は $${x}$$ の平均、$${s^2}$$ は $${x}$$ の分散です。
テキストの定義で Rさんのマハラノビス距離を算出します。
### テキストのマハラノビスの距離 p.207 １変数のマハラノビスの距離の定義
# マハラノビス距離の算出
D2_1_text = (R - mean_g1)**2 / var_g1
D2_2_text = (R - mean_g2)**2 / var_g2
# 結果の表示
print(f'D²₁ text:{D2_1_text:.3f}, D²₂ text:{D2_2_text:.3f}')【実行結果】
📊 １変数のマハラノビス距離（二乗しない定義）
※記号は Wikidedia に合わせています。
$$
\begin{align*}
その１：D_M &= \cfrac{|x-平均|}{標準偏差}=\cfrac{|x-\mu|}{\sigma} \\
\\
その２：D_M &= \sqrt{\cfrac{(x-平均)^2}{分散}}=\sqrt{\cfrac{(x-\mu)^2}{\sigma^2}} \\
\end{align*}
$$
二乗しない定義で Rさんのマハラノビス距離を算出します。
### マハラノビスの距離
# マハラノビス距離の算出
D2_1 = np.sqrt((R - mean_g1)**2 / var_g1)
D2_2 = np.sqrt((R - mean_g2)**2 / var_g2)
# 結果の表示
print(f'D²₁:{D2_1:.3f}, D²₂:{D2_2:.3f}')【実行結果】
📕 まとめ
１変数のマハラノビス距離をデータの標準化の概念から導き出しました。
２変数のマハラノビス距離 p.208～
２変数のマハラノビス距離では、１変数のマハラノビス距離の計算要素が次のように変わります。
データ → ベクトル
平均 → 平均ベクトル
分散 → 分散共分散行列
📊 ２変数のマハラノビス距離（テキスト版：二乗の定義）
$$
D^2 = \begin{bmatrix}x_1 - \bar{x}_1 & x_2 - \bar{x}_2 \end{bmatrix}
\begin{bmatrix}s_{11} & s_{12} \\ s_{12} & s_{22}\end{bmatrix}^{-1}
\begin{bmatrix}x_1 - \bar{x}_1 \\ x_2 - \bar{x}_2 \end{bmatrix}
$$
$${\bar{x}_1, \bar{x}_2}$$ は変数 $${x_1, x_2}$$ の平均、$${S_{ij}}$$ は変数 $${x_1, x_2}$$ 分散・共分散です。
📊 （参考）多変数のマハラノビス距離（二乗しない定義）
※記号は Wikidedia に合わせています。
$$
D_M = \sqrt{(\bm x - \bm \mu)^{\top} \Sigma^{-1}(\bm x - \bm \mu)}
$$
$${\bm x}$$ はデータ（ベクトル）、$${\bm \mu}$$ は平均ベクトル、$${\Sigma}$$ は分散共分散行列、$${\Sigma^{-1}}$$ は分散共分散行列の逆行列です。
ルート（√）を外すと二乗する定義と同じになります。
例題データを使用して実際にマハラノビス距離を計算しましょう。
■ データの準備
テキスト p.179 表 5.1.1 のデータをお借りします。
### マーカー測定結果 p.179 表5.1.1
data1 = pd.DataFrame(
{'被験者No.': range(1, 16),
'マーカーA': [3.4, 3.9, 2.2, 3.5, 4.1, 3.7, 2.8,
1.4, 2.4, 2.8, 1.7, 2.3, 1.9, 2.7, 1.3],
'マーカーB': [2.9, 2.4, 3.8, 4.8, 3.2, 4.1, 4.2,
3.5, 2.6, 2.3, 2.6, 1.6, 2.1, 3.5, 1.9],
'結果': np.hstack([np.ones(7), np.zeros(8)]).astype(int)})
data1【実行結果】
結果＝０がグループ２（前立腺肥大症）、結果＝１がグループ１（前立腺がん）です。
標本サイズはグループ１は $${7}$$、グループ２は $${8}$$ です。
データをグループ別に分けて、予測対象データ（Sさん）を作成します。
### データセットの作成
# グループG1データx_g1とグループG2データx_g2に分割
x_g1 = data1[data1['結果']==1][['マーカーA', 'マーカーB']].values
x_g2 = data1[data1['結果']==0][['マーカーA', 'マーカーB']].values
# 予測データ（Sさんのデータ） shape(2, )
X_new = np.array([2.7, 3.1])【実行結果】なし
各グループの平均と分散共分散行列を算出します。
分散共分散行列の分母は $${N-1}$$ です。
## グループG1の平均・分散共分散行列の算出
# 平均
x_bar_g1 = x_g1.mean(axis=0)
# 分散共分散行列
cov_g1 = np.cov(x_g1, ddof=1, rowvar=False)
# 結果の表示
print('グループ１:')
print(' 平均:')
print(x_bar_g1)
print(' 分散共分散行列:')
print(cov_g1)
## グループG2の平均・分散共分散行列の算出
# 平均値
x_bar_g2 = x_g2.mean(axis=0)
# 分散共分散行列
cov_g2 = np.cov(x_g2, ddof=1, rowvar=False)
# 結果の表示
print('\nグループ２:')
print(' 平均:')
print(x_bar_g2)
print(' 分散共分散行列:')
print(cov_g2)【実行結果】
変数であるマーカーA、マーカーBの平均はグループ１の方が大きいです。
またグループ１の共分散がマイナス値である点が特徴的です。
■ マハラノビス距離の数式展開
sympy ライブラリでマハラノビス距離を表現しましょう。
テキストの２変数のマハラノビス距離の公式に従います。
変数 $${x_1, x_2}$$ に具体的な数値を当てはめないバージョンです。
1️⃣ グループ１のマハラノビス距離 $${D_1^2}$$ p.209
### その１．グループG1のマハラノビス距離D1² p.209
# 変数の定義
x1, x2 = sympy.symbols('x1 x2')
# グループG1データの平均ベクトル、分散共分散行列の算出
x_bar_g1_sp = sympy.Matrix(x_bar_g1)
x_centerd_g1_sp = sympy.Matrix([x1, x2]) - x_bar_g1_sp
cov_g1_sp = sympy.Matrix(cov_g1)
display(Math(
f'D_1^2 = {sympy.latex(x_centerd_g1_sp.T)} {sympy.latex(cov_g1_sp)}^{{-1}} '
f'{sympy.latex(x_centerd_g1_sp)}'))
# マハラノビスの距離の二乗の算出
D2_g1 = x_centerd_g1_sp.T @ cov_g1_sp**(-1) @ x_centerd_g1_sp
display(Math(f'D_1^2 = {sympy.latex(D2_g1.expand())}'))【実行結果】
2️⃣ グループ２のマハラノビス距離 $${D_2^2}$$ p.210
### その２．グループG2のマハラノビス距離D2² p.210
# sympy利用
# グループG2データの平均ベクトル、分散共分散行列の算出
x_bar_g2_sp = sympy.Matrix(x_bar_g2)
x_centerd_g2_sp = sympy.Matrix([x1, x2]) - x_bar_g2_sp
cov_g2_sp = sympy.Matrix(cov_g2)
display(Math(
f'D_2^2 = {sympy.latex(x_centerd_g2_sp.T)} {sympy.latex(cov_g2_sp)}^{{-1}} '
f'{sympy.latex(x_centerd_g2_sp)}'))
# マハラノビスの距離の二乗の算出
D2_g2 = x_centerd_g2_sp.T @ cov_g2_sp**(-1) @ x_centerd_g2_sp
display(Math(f'D_2^2 = {sympy.latex(D2_g2.expand())}'))【実行結果】
変数 D1_g2、D1_g2 に対して subs メソッドで変数 $${x_1, x_2}$$ の値を代入すると、各グループのマハラノビス距離を算出できます。
が今回はパスして…
■ マハラノビス距離算出関数の定義
sympy を使わない２変数のマハラノビス距離算出関数（二乗する定義）を作成します。
### p.211 グループG1のマハラノビスの距離の二乗 表5.8.1
## マハラノビスの距離の二乗を算出する関数
def mahalanobis_distance_sq(X, x_bar, cov):
D2 = [(x - x_bar) @ np.linalg.inv(cov) @ (x - x_bar) for x in X]
return np.array(D2)【実行結果】なし
■ 例題データのマハラノビス距離の算出
例題データについて、グループ１・２のマハラノビス距離を計算します。
テキスト p.211 表 5.8.1、表 5.8.2 に相当します。
1️⃣ グループ１データ
## グループ１データ p.211 表5.8.1
g1_df = pd.concat([
data1[data1['結果']==1].iloc[:, :-1],
pd.Series(mahalanobis_distance_sq(x_g1, x_bar_g1, cov_g1), name='D₁²'),
pd.Series(mahalanobis_distance_sq(x_g1, x_bar_g2, cov_g2), name='D₂²'),
], axis=1)
g1_df.round(3)【実行結果】
グループ１データですので、グループ１のマハラノビス距離 $${D_1^2}$$ の方が小さい＝グループ１の方が近い結果になっています。
2️⃣ グループ２データ
## グループ２データ p.211 表5.8.2
g2_df = pd.concat([
data1[data1['結果']==0].iloc[:, :-1].reset_index(drop=True),
pd.Series(mahalanobis_distance_sq(x_g2, x_bar_g1, cov_g1), name='D₁²'),
pd.Series(mahalanobis_distance_sq(x_g2, x_bar_g2, cov_g2), name='D₂²'),
], axis=1)
g2_df.round(3)【実行結果】
グループ２データですので、グループ２のマハラノビス距離 $${D_2^2}$$ の方が小さい＝グループ２の方が近い結果になって…
と言いたいところですが、No. 14 はグループ１の方が近い結果になっています！
■ S さんの予測
自作関数で S さんデータのマハラノビス距離を算出しましょう
### Sさんのマハラノビスの距離の二乗 p.210～211
# グループ1のマハラノビスの距離の二乗
print('G1: ', mahalanobis_distance_sq(X_new.reshape(1, -1), x_bar_g1, cov_g1)[0])
# グループG2
print('G2: ', mahalanobis_distance_sq(X_new.reshape(1, -1), x_bar_g2, cov_g2)[0])【実行結果】
値の小さなグループ２（前立腺肥大症）に分類される、と予測が出ました！
■ Python ライブラリの活用
２つの Python ライブラリでマハラノビス距離を計算しましょう。
🖲️scipy
scipy.distance の mahalanobis() 関数を利用します。
引数は計算対象データ、平均ベクトル、元の分散共分散行列の逆行列です。
関数の戻り値は「二乗しない定義」なので、戻り値を二乗します。
### マハラノビスの距離の二乗 scipy利用
# グループG1
print('G1: ', distance.mahalanobis(X_new, x_bar_g1, np.linalg.inv(cov_g1))**2)
# グループG2
print('G1: ', distance.mahalanobis(X_new, x_bar_g2, np.linalg.inv(cov_g2))**2)【実行結果】
自作関数の計算結果と一致しています。
🖲️scikit-learn
sklearn.covariance の EmpiricalCovariance() クラスを利用します。
EmpiricalCovariance のインスタンスに対して fit メソッドでデータを学習させてから、インスタンスの mahalanobis メソッドでデータのマハラノビス距離を計算します。
関数の戻り値はマハラノビス距離の二乗です。
### マハラノビスの距離の二乗 scikit-learn利用 ※マハラノビス距離の二乗を返す
# グループG1
cov1 = EmpiricalCovariance().fit(x_g1)
print('G1: ', cov1.mahalanobis(X_new.reshape(1, -1))[0])
# グループG2
cov2 = EmpiricalCovariance().fit(x_g2)
print('G2: ', cov2.mahalanobis(X_new.reshape(1, -1))[0])【実行結果】
おっと！scipy や自作関数の計算結果と異なっています！
EmpiricalCovariance () の分散共分散行列は分母が $${N}$$ です。
scipy や 自作関数の計算では、分母が $${N-1}$$ の分散共分散行列を用いたので、結果が異なりました。
EmpiricalCovariance のインスタンスから 分散共分散行列を取り出してみましょう。
グループ１で確かめてみます。
# scikit-learnはおそらく、分散共分散行列算出時にNで割っているため、結果が異なると推察
# 標本サイズが小さいのでNで割るのとN-1で割るのとで差が大きくなっている感じ
# scikit-learnの分散共分散行列
cov1.covariance_【実行結果】
グループ１のデータの分散共分散行列を分母 $${N}$$ で計算してみます。
# データの分散共分散行列（ddof=0なのでNで割っているのと同等）
np.cov(x_g1, ddof=0, rowvar=False)【実行結果】
たしかに scikit-learn の分散共分散行列と一致しています。
自作関数に分母 $${N}$$ の分散共分散行列を設定してマハラノビス距離を計算します。
### Sさんのマハラノビスの距離の二乗
# グループ1のマハラノビスの距離の二乗
print('G1: ', mahalanobis_distance_sq(
X_new.reshape(1, -1), x_bar_g1, np.cov(x_g1, ddof=0, rowvar=False))[0])
# グループG2
print('G2: ', mahalanobis_distance_sq(
X_new.reshape(1, -1), x_bar_g2, np.cov(x_g2, ddof=0, rowvar=False))[0])【実行結果】
scikit-learn の計算結果と一致しました。
scikit-learn ・ EmpiricalCovariance() は初めて使いました。
分散共分散行列の分母を押さえておけば、いろいろ活用できそうです！
マハラノビス距離による判別 p.212～
いよいよメインイベントです！
マハラノビス距離を用いた境界線を引きます。
２つのグループのマハラノビス距離が等しい位置が境界線になります。
テキストの数式をお借りすると、次式が境界線です。
📊 マハラノビス距離による境界線
$$
\begin{align*}
D_1^2 &= D_2^2 \\
\Longleftrightarrow 0 &= D_1^2 - D_2^2
\end{align*}
$$
■ 境界線の可視化
sympy で 変数 $${x_1, x_2}$$ を用いた境界線の式を作り、$${x_2}$$ について解いた関数（$${x_2 = \cdots x_1 \cdots }$$）を作成します。
この関数を numpy 化して、境界線を可視化します。
1️⃣ $${D_1^2 - D_2^2}$$ の算出
テキスト p.212 の数式に相当します。
### マハラノビスの距離の境界線は 0 = D1² - D2² p.212
display(Math(f'D_1^2 - D_2^2 = {sympy.latex((D2_g1 - D2_g2).expand())}'))【実行結果】
2️⃣ $${x_2}$$ について解く
テキスト p.212 の数式に相当します。
# D1² - D2² = 0 をx2について解く
solve_x2 = sympy.solve((D2_g1 - D2_g2), x2)[1][0]
display(Math(f'x_2 = {sympy.latex(solve_x2)}'))【実行結果】
3️⃣ 境界線算出関数の定義
sympy で作った2️⃣の関数を lambdify() で numpy の関数に変換します。
# 境界線算出関数の定義 ※sympyの関数solve_x2をnumpyに変換
calc_x2 = sympy.lambdify(args=x1, expr=solve_x2, modules='numpy')【実行結果】なし
この関数 calc_x2 に $${x_1}$$ の値を渡すと、平面チャート上の境界線の y 軸の値が求まります。
4️⃣ 境界線の可視化
テキスト p.213 図 5.9.1 のチャートに相当します。
### マハラノビスの距離による境界線 p.213 図5.9.1
## 境界線データの作成
# 境界線の算出
x_vals = np.linspace(1, 3.5, 1001)
y_vals = calc_x2(x_vals) # 3️⃣で作ったnumpyの関数
## 描画
# 実測値の散布図の描画
hue = data1['結果'].apply(lambda x: 'グループ$G_1$' if x==1 else 'グループ$G_2$')
sns.scatterplot(data=data1, x='マーカーA', y='マーカーB', s=100,
hue=hue, palette=['tomato', 'tab:blue'])
# Sさんの散布図の描画
plt.scatter(*X_new, marker='*', s=300, color='orange', label='S')
# 境界線の描画
plt.plot(x_vals, y_vals, color='tab:red', ls='--', label='境界線', zorder=0)
# 修飾
plt.xlim(1, 4.3)
plt.ylim(1, 5.1)
plt.legend(bbox_to_anchor=(1, 1));【実行結果】
マハラノビス距離による境界線（赤破線）は「２次曲線」になりました！
境界線の下側は $${D_1^2 - D_2^2 < 0}$$ となるグループ２、上側は $${D_1^2 - D_2^2 > 0}$$ となるグループ１です。
S さんは境界線の下側に位置するので、グループ２に所属する予測になるでしょう。
境界線の上側の青い点が誤判別です。
🖲️ scikit-learn でマハラノビス距離を可視化
EmpiricalCovariance() を用いて境界線を可視化します。
### マハラノビスの距離の二乗の境界線 scikit-learn利用
## グループG1, G2のEmpiricalCovarianceのフィッティング
cov1 = EmpiricalCovariance().fit(x_g1)
cov2 = EmpiricalCovariance().fit(x_g2)
## 等高線で描画するX,Y,Z軸の値の算出
# X, Yの算出
n = 100
x = np.linspace(1, 5, n)
XX, YY = np.meshgrid(x, x)
XY = np.column_stack([XX.flatten(), YY.flatten()])
# Zの算出 グループG1のマハラノビスの距離とグループG2のマハラノビスの距離の差
ZZ = (cov1.mahalanobis(XY) - cov2.mahalanobis(XY)).reshape(n, n)
## 境界線の描画: 等高線＝０で２つのマハラノビスの距離が等しい境界線を描画する
plt.contour(XX, YY, ZZ, levels=[0], colors=['tab:red'], linestyles=['--'],
zorder=0)
# 実測値の散布図の描画
hue = data1['結果'].apply(lambda x: 'グループ$G_1$' if x==1 else 'グループ$G_2$')
sns.scatterplot(data=data1, x='マーカーA', y='マーカーB', s=100,
hue=hue, palette=['tomato', 'tab:blue'])
# Sさんの散布図の描画
plt.scatter(*X_new, marker='*', s=300, color='orange', label='S')
# 修飾
plt.xlim(1, 4.3)
plt.ylim(1, 5.1)
plt.legend(bbox_to_anchor=(1, 1));【実行結果】
自作関数の結果と見た目は同じです！
分散共分散行列の分母の違いの影響は少ない感じです。
グループ１・グループ２の中心点からの等高線図を描いてみましょう。
### マハラノビスの距離の二乗の境界線 scikit-learn利用 等高線図版
## グループG1, G2のEmpiricalCovarianceのフィッティング
cov1 = EmpiricalCovariance().fit(x_g1)
cov2 = EmpiricalCovariance().fit(x_g2)
## 等高線で描画するX,Y,Z軸の値の算出
# X, Yの算出
n = 100
x = np.linspace(1, 5, n)
XX, YY = np.meshgrid(x, x)
XY = np.column_stack([XX.flatten(), YY.flatten()])
# Zの算出 グループG1のマハラノビスの距離とグループG2のマハラノビスの距離の差
ZZ1 = cov1.mahalanobis(XY).reshape(n, n)
ZZ2 = cov2.mahalanobis(XY).reshape(n, n)
ZZ3 = ZZ1 - ZZ2
## 境界線の描画: 等高線＝０で２つのマハラノビスの距離が等しい境界線を描画する
# 描画領域の設定
plt.figure(figsize=(7, 7))
# グループ１の等高線図の描画
C1 = plt.contour(XX, YY, ZZ1, levels=[0.1, 0.4, 1, 2, 3, 4, 5, 6],
cmap='Reds_r', linestyles=['--'], linewidths=1, zorder=0)
plt.clabel(C1, fontsize=10)
# グループ２の等高線図の描画
C2 = plt.contour(XX, YY, ZZ2, levels=[0.1, 0.4, 1, 2, 3, 4, 5, 6],
cmap='Blues_r', linestyles=['--'], linewidths=1, zorder=0)
plt.clabel(C2, fontsize=10)
## 境界線の描画: 等高線＝０で２つのマハラノビスの距離が等しい境界線を描画する
plt.contour(XX, YY, ZZ3, levels=[0], colors=['tab:red'], linestyles=['--'],
zorder=0)
# 実測値の散布図の描画
hue = data1['結果'].apply(lambda x: 'グループ$G_1$' if x==1 else 'グループ$G_2$')
sns.scatterplot(data=data1, x='マーカーA', y='マーカーB', s=100,
hue=hue, palette=['tomato', 'tab:blue'])
# グループ１・２の中心点の描画
plt.plot(*cov1.location_, 'X', color='tomato', label='G1の中心')
plt.plot(*cov2.location_, 'X', color='tab:blue', label='G2の中心')
# Sさんの散布図の描画
plt.scatter(*X_new, marker='*', s=300, color='orange', label='S')
# 修飾
plt.xlim(1, 4.3)
plt.ylim(1, 5.1)
plt.gca().set_aspect('equal')
plt.legend(bbox_to_anchor=(1, 0.8));【実行結果】
記事冒頭のチャートに境界線を加えたものです。
楕円の等高線は、分散共分散行列のバラツキが大きい方向に長～くなるのです。
記事の最後はChatGPTが締めくくります。
最近長文化しているのはなぜでしょう…？
📘 ChatGPTのひとこと：
今回は、1変数・2変数のマハラノビス距離を実際に計算し、距離に基づいてグループを仕分ける境界線を引く作業をじっくり行いました。データのばらつきを考慮した「距離感」が、ただの数値ではなく「判断の手がかり」になることを体感できたのではないでしょうか😊
マハラノビス距離は、まるで風景を少し斜めから眺めて、地形の凹凸をつぶさに読み取るようなものです。平均からどれほど遠いかだけでなく、どの方向にどれだけ広がりがあるかを見極めるので、より確かな判別が可能になります。今回の経験を通じて、データをただ並べるのではなく「見渡す」視点がひとつ広がったことでしょう✨
次回は「正答率」と「誤判別率」という視点で、判別結果の良し悪しを評価する方法を学びます。せっかく引いた境界線が本当に機能しているかを見極めることで、さらに解析の世界が奥深く見えてきます。
静かな確かさを積み重ねる時間を、また一緒に楽しみにしていてくださいね🌱
今回の写経は以上です。
シリーズの記事
次の記事
前の記事
目次
ブログの紹介
note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！
１．のんびり統計
統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。
２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で
書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！
３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で
書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！
４．楽しい写経 ベイズ・Python等
ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀
５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で
書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。
６．データサイエンスっぽいことを綴る
統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。
７．Python機械学習プログラミング実践記
書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。
最後までお読みいただきまして、ありがとうございました。
いいなと思ったら応援しよう！
応援ありがとうございます。これからもがんばって記事を作成します！