はじめに

　機械学習で登場する確率分布について色々な角度から理解したいシリーズです。

　この記事では、マハラノビス距離について、図を使って解説します。

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

はじめに
2次元マハラノビス距離の作図
参考文献
おわりに

2次元マハラノビス距離の作図

　2次元変数のマハラノビス距離(2D Mahalanobis' distance)のグラフを作成します。
　マハラノビス距離については「マハラノビス距離の性質の導出 - からっぽのしょこ」、3次元変数の場合については「【Python】3次元マハラノビス距離の作図 - からっぽのしょこ」を参照してください。

　利用するライブラリを読み込みます。

# ライブラリを読込
import numpy as np
#from scipy.stats import wishart
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

マハラノビス距離の定義式

　まずは、マハラノビス距離の定義を数式で確認します。

　 $D$ 次元の変数を $\mathbf{x}$ 、平均ベクトルを $\boldsymbol{\mu}$ 、分散共分散行列を $\boldsymbol{\Sigma}$ とします。

$\displaystyle \mathbf{x} = \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_D \end{bmatrix} ,\ \boldsymbol{\mu} = \begin{bmatrix} \mu_1 \\ \mu_2 \\ \vdots \\ \mu_D \end{bmatrix} ,\ \boldsymbol{\Sigma} = \begin{bmatrix} \sigma_{1,1} & \sigma_{1,2} & \cdots & \sigma_{1,D} \\ \sigma_{2,1} & \sigma_{2,2} & \cdots & \sigma_{2,D} \\ \vdots & \vdots & \ddots & \vdots \\ \sigma_{D,1} & \sigma_{D,2} & \cdots & \sigma_{D,D} \end{bmatrix}$

　 $\mu_d$ は $x_d$ の平均、 $\sigma_d$ は $x_d$ の標準偏差、 $\sigma_d^2 = \sigma_{d,d}$ は $x_d$ の分散、 $\sigma_{i,j} = \sigma_{j,i}$ は $x_i, x_j$ の共分散です。平均パラメータ $\mu_d$ は実数、分散パラメータ $\sigma_d^2$ は正の実数、共分散パラメータ $\sigma_{i,j}$ は実数、また $\boldsymbol{\Sigma}$ は正定値行列を満たす必要があります。

　点 $\boldsymbol{\mu}$ を中心とする点 $\mathbf{x}$ のマハラノビス距離は、精度行列(分散共分散行列の逆行列)を用いた二次形式のルートで定義されます。

$\displaystyle \Delta = \sqrt{ (\mathbf{x} - \boldsymbol{\mu})^{\top} \boldsymbol{\Sigma}^{-1} (\mathbf{x} - \boldsymbol{\mu}) }$

　変数とパラメータの関係やパラメータと二次形式の関係については「変数ベクトルと平均ベクトル・分散共分散行列の関係の導出 - からっぽのしょこ」を参照してください。

　この式をグラフで表現していきます。

マハラノビス距離の図

　次は、2次元空間におけるマハラノビス距離のグラフを作成します。
　この例では、2次元の図で可視化するため、変数の次元が2の場合のみ処理できます。また、Pythonのインデックスに合わせて、添字を0から割り当てます。

　パラメータを設定します。

# 次元数を指定:(固定)
D = 2

# 平均ベクトルを指定
mu_d = np.array([0.0, 0.0])

# 分散共分散行列を指定
sigma_dd = np.array(
    [[9.0, -2.3], 
     [-2.3, 4.0]]
)
#sigma_dd = np.diag(np.ones(D)) # 単位行列
#sigma_dd = wishart.rvs(df=D, scale=np.diag(np.ones(D)), size=1) # ウィシャート乱数
print(sigma_dd.round(2))

[[ 9.  -2.3]
 [-2.3  4. ]]

　次元数を $D = 2$ として、平均ベクトル $\boldsymbol{\mu} = (\mu_0, \mu_1)$ 、分散共分散行列 $\boldsymbol{\Sigma} = ( (\sigma_0^2, \sigma_{0,1}), (\sigma_{1,0}, \sigma_1^2) )$ の値を指定します。

　距離の計算用の座標を作成します。

# 各軸の範囲の調整値を指定
sgm_num = 3.0

# 各軸の範囲を設定
x0_min = mu_d[0] - sgm_num * np.sqrt(sigma_dd[0, 0])
x0_max = mu_d[0] + sgm_num * np.sqrt(sigma_dd[0, 0])
x1_min = mu_d[1] - sgm_num * np.sqrt(sigma_dd[1, 1])
x1_max = mu_d[1] + sgm_num * np.sqrt(sigma_dd[1, 1])
print(x0_min.round(2), x0_max.round(2))
print(x1_min.round(2), x1_max.round(2))

# 各軸の値を作成
x0 = np.linspace(start=x0_min, stop=x0_max, num=50) # 0軸方向の点の数を指定
x1 = np.linspace(start=x1_min, stop=x1_max, num=50) # 1軸方向の点の数を指定
print(x0[:5].round(2))
print(x1[:5].round(2))

# 格子点を作成
X0, X1 = np.meshgrid(x0, x1)
print(X0[:5, :5].round(2))
print(X1[:5, :5].round(2))

# 格子点の形状を設定
grid_shape = X0.shape
print(grid_shape)

# 座標を作成
X = np.stack([X0.flatten(), X1.flatten()], axis=1)
print(X[:5].round(2))
print(X.shape)

-9.0 9.0
-6.0 6.0
[-9.   -8.63 -8.27 -7.9  -7.53]
[-6.   -5.76 -5.51 -5.27 -5.02]
[[-9.   -8.63 -8.27 -7.9  -7.53]
 [-9.   -8.63 -8.27 -7.9  -7.53]
 [-9.   -8.63 -8.27 -7.9  -7.53]
 [-9.   -8.63 -8.27 -7.9  -7.53]
 [-9.   -8.63 -8.27 -7.9  -7.53]]
[[-6.   -6.   -6.   -6.   -6.  ]
 [-5.76 -5.76 -5.76 -5.76 -5.76]
 [-5.51 -5.51 -5.51 -5.51 -5.51]
 [-5.27 -5.27 -5.27 -5.27 -5.27]
 [-5.02 -5.02 -5.02 -5.02 -5.02]]
(50, 50)
[[-9.   -6.  ]
 [-8.63 -6.  ]
 [-8.27 -6.  ]
 [-7.9  -6.  ]
 [-7.53 -6.  ]]
(2500, 2)

　2次元空間における点 $\mathbf{x} = (x_0, x_1)$ の座標を作成します。
　この例では、各軸( $d$ 番目の軸・ $x_d$ 軸)の範囲を平均 $\mu_d$ を中心に標準偏差 $\sigma_d$ の sgm_num 倍とします。
　第0軸(x軸)と第1軸(y軸)の平面における格子状の点 $(x_0, x_1)$ を np.meshgrid() で作成します。点(0・1軸の値)の数が曲線の滑らかさに影響します。
　作図処理が重い場合は、値(点)の数を調整してください。

　2次元空間におけるマハラノビス距離を計算します。

# 精度行列を計算
inv_sigma_dd = np.linalg.inv(sigma_dd)

# マハラノビス距離を計算
dist_vec = np.array(
    [np.sqrt((x_d - mu_d).T @ inv_sigma_dd @ (x_d - mu_d)) for x_d in X]
)
print(dist_vec[:5].round(2))
print(dist_vec.shape)

[5.4  5.29 5.18 5.08 4.97]
(2500,)

　 $\mathbf{x}$ のマハラノビス距離 $\Delta$ を計算します。

　for文を使わない場合は、次のようにも処理できます。

# マハラノビス距離を計算
dist_vec = np.sqrt(
    np.diag((X - mu_d) @ inv_sigma_dd @ (X - mu_d).T)
)
print(dist_vec[:5].round(2))
print(dist_vec.shape)

[5.4  5.29 5.18 5.08 4.97]
(2500,)

　X の行(各点 $\mathbf{x}$ )に関して全ての組み合わせで計算して、対角要素(同じ点)の要素を取り出します。

　2次元空間におけるマハラノビス距離の等高線図を作成します。

# ラベル用の文字列を作成
def_label = '$\\Delta = \\sqrt{(x - \\mu)^{T} \\Sigma^{-1} (x - \\mu)}$'
mu_str    = '(' + ', '.join(str(val.round(2)) for val in mu_d) + ')'
sigma_str = '(' + ', '.join('(' + ', '.join(str(val.round(2)) for val in vec) + ')' for vec in sigma_dd) + ')'
param_label  = '$D = 2$\n'
param_label += '$\\mu = ' + mu_str + '$\n'
param_label += '$\\Sigma = ' + sigma_str + '$'

# 2Dマハラノビス距離を作図
fig, ax = plt.subplots(figsize=(8, 6), dpi=100, facecolor='white', constrained_layout=True)
cs = ax.contour(X0, X1, dist_vec.reshape(grid_shape), 
                cmap='viridis') # 等高線
ax.set_xlim(xmin=x0_min, xmax=x0_max)
ax.set_ylim(ymin=x1_min, ymax=x1_max)
ax.set_xlabel('$x_0$')
ax.set_ylabel('$x_1$')
ax.set_title(param_label, loc='left')
fig.suptitle("Mahalanobis' Distance", fontsize=20)
fig.colorbar(cs, ax=ax, shrink=0.8, label=def_label)
ax.grid()
ax.set_aspect('equal', adjustable='box')
plt.show()

　平均ベクトルの点 $\boldsymbol{\mu}$ を中心として、同じ距離の点(位置)を結んだ曲線でマハラノビス距離の形状を、線の色で距離の大きさを示します。
　中心から離れるほど距離が大きくなります。

　2次元空間におけるマハラノビス距離の曲面図を作成します。

# グラデーションの範囲を設定
dist_min = 0.0
dist_max = np.ceil(dist_vec.max())

# 軸サイズを設定
x0_size = x0_max - x0_min
x1_size = x1_max - x1_min

# 2Dマハラノビス距離を作図
fig, ax = plt.subplots(figsize=(8, 8), dpi=100, facecolor='white', constrained_layout=True, 
                       subplot_kw={'projection': '3d'})
ax.contour(X0, X1, dist_vec.reshape(grid_shape), offset=0.0, 
           cmap='viridis', vmin=dist_min, vmax=dist_max) # 等高線
cs = ax.plot_surface(X0, X1, dist_vec.reshape(grid_shape), 
                     cmap='viridis', vmin=dist_min, vmax=dist_max, alpha=0.9) # 曲面
ax.set_xlim(xmin=x0_min, xmax=x0_max)
ax.set_ylim(ymin=x1_min, ymax=x1_max)
ax.set_xlabel('$x_0$')
ax.set_ylabel('$x_1$')
ax.set_zlabel('$\\Delta$')
ax.set_title(param_label, loc='left')
fig.suptitle("Mahalanobis' Distance", fontsize=20)
fig.colorbar(cs, ax=ax, shrink=0.8, pad=0.1, label=def_label)
ax.set_box_aspect([1.0, x1_size/x0_size, 0.6]) # 高さ(横サイズに対する比)を指定
#ax.view_init(elev=90, azim=270) # 0・1軸平面
#ax.view_init(elev=0, azim=270) # 0・2軸平面
#ax.view_init(elev=0, azim=0) # 1・2軸平面
plt.show()

　マハラノビス距離の大きさを高さ(z軸の値)で示します。
　目安として、底面(x・y軸平面)や曲面上に等高線図を表示しています。

　以上で、マハラノビス距離のグラフを作成できました。

パラメータと形状の関係

　最後は、パラメータの値を少しずつ変化させて、距離の形状の変化をグラフで確認します。
　アニメーションの作図コードは「Probability-Distribution/code/mahalanobis_distance/2d_parameter.py at main · anemptyarchive/Probability-Distribution · GitHub」を参照してください。

平均パラメータの影響

　平均パラメータとマハラノビス距離の関係をアニメーションで確認します。

　フレームごとの平均パラメータ $\boldsymbol{\mu}$ をバツ印、推移を破線で示します。

　点 $\boldsymbol{\mu}$ は、マハラノビス距離が $\Delta = 0$ の点であり、等高線の中心です。 $\boldsymbol{\mu}$ が変化すると、形状は変わらずに全体が移動します。

分散パラメータの影響

　分散パラメータとマハラノビス距離の関係をアニメーションで確認します。

　フレームごとの分散パラメータ $\sigma_d^2$ を、軸線に平行で中心からの長さが標準偏差パラメータ $\sigma_d$ の線分で示します。全ての共分散が $\sigma_{i,j} = 0$ のとき、 $\Delta = 1$ の線と一致します。
　比較用に、ユークリッド距離(全ての方向に等間隔の距離・分散共分散行列が単位行列のときのマハラノビス距離)を網目の曲面と破線で表示しています。

　分散 $\sigma_d^2$ が大きく(精度 $\lambda_{d,d}$ が小さく)なるほど、同じ位置のマハラノビス距離が大きくなり、等高線が拡がり(同じ距離の位置が遠のき)ます。
　マハラノビス距離とユークリッド距離の関係については「マハラノビス距離の性質の導出」、マハラノビス距離と標準偏差の関係については「変数ベクトルと平均ベクトル・分散共分散行列の関係の導出」を参照してください。

共分散パラメータの影響

　共分散パラメータとマハラノビス距離の関係をアニメーションで確認します。

　フレームごとの共分散パラメータ $\sigma_{i,j}$ の影響を、楕円の軸で示します。分散共分散行列 $\boldsymbol{\Sigma}$ の固有値 $\lambda_d$ と固有ベクトル $\mathbf{u}_d$ を用いて描画できます。 $x_d, \mu_d$ と $\lambda_d, \mathbf{u}_d$ の次元番号(軸番号) $d$ は無関係です。固有値と精度パラメータは(同じ記号を使いましたが全くの)無関係です。(なんか軸の色が入れ替わりますが今回は無視しました…)
　比較用に、ユークリッド距離を網目の曲面と破線で表示しています。

　分散共分散行列 $\boldsymbol{\Sigma}$ によって楕円の形状が変わり、軸の長さや傾きが変わります。 $x_i, x_j$ の共分散が正の値 $\sigma_{i,j} \gt 0$ のとき「 $i$ 軸が正の値 $x_i \gt 0$ と $j$ 軸が正の値 $x_j \gt 0$ 」と「 $i$ 軸が負の値 $x_i \lt 0$ と $j$ 軸が負の値 $x_j \lt 0$ 」の( $i, j$ 軸の符号が同じ)方向に、負の値 $\sigma_{i,j} \lt 0$ のとき「 $i$ 軸が正の値 $x_i \gt 0$ と $j$ 軸が負の値 $x_j \lt 0$ 」と「 $i$ 軸が負の値 $x_i \lt 0$ と $j$ 軸が正の値 $x_j \gt 0$ 」の( $i, j$ 軸の符号が異なる)方向に伸びます。
　分散共分散行列と固有値・固有ベクトルの関係「2.3.0：分散共分散行列と固有値・固有ベクトルの関係の導出【PRMLのノート】 - からっぽのしょこ」を参照してください。

　この記事では、2変数の場合のマハラノビス距離を図で確認しました。次の記事では、3変数の場合を図で確認します。

参考文献

パターン認識と機械学習上

作者:C.M. ビショップ
丸善出版

Amazon

おわりに

　線形回帰モデルの最小二乗法の解説にて3次元のマハラノビス距離のグラフを作りまして、その作図処理を次回の記事にまとめました。ただ2次元版(ポピュラーなグラフ)を飛ばして3次元版(あまり見ないグラフ)を載せてもなと、あとからこの記事を書きました。2次元版の内容はこれまでも(たぶん)扱った内容なので書いてて楽しくないし、このブログのコンセプト的にも図解だけだと物足りないしで、数式パートの記事もついでに書いたり書き直したりしたのが前回の記事です。
　2次元版が書けてからは、1次元変数のマハラノビス距離はどう図示できるのか気になってます。計算自体は標準化っぽい処理になるんですかね。今回は挑戦しませんでした。

【次の内容】

からっぽのしょこ

読んだら書く！書いたら読む！同じ事は二度調べ(たく)ない

【Python】2次元マハラノビス距離の作図