はじめに

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

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

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

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

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

　3次元変数のマハラノビス距離(3D Mahalanobis' distance)のグラフを作成します。
　マハラノビス距離については「マハラノビス距離の性質の導出 - からっぽのしょこ」、2次元変数の場合については「【Python】2次元マハラノビス距離の作図 - からっぽのしょこ」、3変数関数のグラフについては「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}) }$

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

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

マハラノビス距離の図

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

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

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

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

# 分散共分散行列を指定
sigma_dd = np.array(
    [[16.0, -1.5, 0.0], 
     [-1.5, 9.0, 0.9], 
     [0.0, 0.9, 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))

[[16.  -1.5  0. ]
 [-1.5  9.   0.9]
 [ 0.   0.9  4. ]]

　次元数を $D = 3$ として、平均ベクトル $\boldsymbol{\mu} = (\mu_0, \mu_1, \mu_2)$ 、分散共分散行列 $\boldsymbol{\Sigma} = ( (\sigma_0^2, \sigma_{0,1}, \sigma_{0,2}), (\sigma_{1,0}, \sigma_1^2, \sigma_{1,2}), (\sigma_{2,0}, \sigma_{2,1}, \sigma_2^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])
x2_min = mu_d[2] - sgm_num * np.sqrt(sigma_dd[2, 2])
x2_max = mu_d[2] + sgm_num * np.sqrt(sigma_dd[2, 2])
print(x0_min.round(2), x0_max.round(2))
print(x1_min.round(2), x1_max.round(2))
print(x2_min.round(2), x2_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軸方向の点の数を指定
x2 = np.linspace(start=x2_min, stop=x2_max, num=100) # 2軸方向の線の数を指定
print(x0[:5].round(2))
print(x1[:5].round(2))
print(x2[:5].round(2))

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

# 0・1軸の形状を設定
grid_shape = X0.shape
grid_size  = X0.size
print(grid_shape)
print(grid_size)

-12.0 12.0
-9.0 9.0
-6.0 6.0
[-12.   -11.51 -11.02 -10.53 -10.04]
[-9.   -8.63 -8.27 -7.9  -7.53]
[-6.   -5.88 -5.76 -5.64 -5.52]
[[-12.   -11.51 -11.02 -10.53 -10.04]
 [-12.   -11.51 -11.02 -10.53 -10.04]
 [-12.   -11.51 -11.02 -10.53 -10.04]
 [-12.   -11.51 -11.02 -10.53 -10.04]
 [-12.   -11.51 -11.02 -10.53 -10.04]]
[[-9.   -9.   -9.   -9.   -9.  ]
 [-8.63 -8.63 -8.63 -8.63 -8.63]
 [-8.27 -8.27 -8.27 -8.27 -8.27]
 [-7.9  -7.9  -7.9  -7.9  -7.9 ]
 [-7.53 -7.53 -7.53 -7.53 -7.53]]
(50, 50)
2500

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

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

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

# 高さごとに処理
dist_arr = np.tile(np.nan, reps=(x2.size, grid_size)) # 受け皿を初期化
for i in range(x2.size):

    # 座標を作成
    X = np.stack([X0.flatten(), X1.flatten(), x2[i].repeat(grid_size)], axis=1)
    
    # マハラノビス距離を計算
    dist_vec = np.array(
        [np.sqrt((x_d - mu_d).T @ inv_sigma_dd @ (x_d - mu_d)) for x_d in X]
    )

    # 距離を格納
    dist_arr[i] = dist_vec.copy()
print(dist_arr[:5, :5].round(2))
print(dist_arr.shape)

[[5.19 5.11 5.03 4.95 4.88]
 [5.16 5.08 5.   4.92 4.85]
 [5.13 5.05 4.97 4.89 4.82]
 [5.1  5.02 4.94 4.86 4.79]
 [5.07 4.99 4.91 4.83 4.76]]
(100, 2500)

　第0・1軸(x・y軸)平面における格子点と第2軸(z軸)の点を結合した点 $\mathbf{x} = (x_0, x_1, x_2)$ を作成して、マハラノビス距離 $\Delta$ を計算します。
　第2軸の値 x2 の要素ごとに、平面上の点の距離を計算して dist_arr に格納していきます。

　3次元空間におけるマハラノビス距離のグラフを作成します。

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

# 等高線の位置を指定
dist_levels = np.linspace(start=dist_min, stop=dist_max, num=7) # 線の数を指定

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

# ラベル用の文字列を作成
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 = 3$\n'
param_label += '$\\mu = ' + mu_str + '$\n'
param_label += '$\\Sigma = ' + sigma_str + '$'

# 3Dマハラノビス距離を作図
fig, ax = plt.subplots(figsize=(9, 8), dpi=100, facecolor='white', constrained_layout=True, 
                       subplot_kw={'projection': '3d'})
cs = ax.contour(X0, X1, np.linspace(dist_min, dist_max, num=grid_size).reshape(grid_shape), 
                cmap='viridis', vmin=dist_min, vmax=dist_max, levels=dist_levels) # カラーバー表示用のダミー
fig.colorbar(cs, ax=ax, shrink=0.6, pad=0.1, label=def_label)
ax.cla() # カラーバー表示用
for i in range(x2.size):
    ax.contour(X0, X1, dist_arr[i].reshape(grid_shape), offset=x2[i], 
               cmap='viridis', vmin=dist_min, vmax=dist_max, levels=dist_levels, alpha=0.5, 
               linewidths=1) # 等高線
ax.set_xlim(xmin=x0_min, xmax=x0_max)
ax.set_ylim(ymin=x1_min, ymax=x1_max)
ax.set_zlim(zmin=x2_min, zmax=x2_max)
ax.set_xlabel('$x_0$')
ax.set_ylabel('$x_1$')
ax.set_zlabel('$x_2$')
ax.set_title(param_label, loc='left')
fig.suptitle("Mahalanobis' Distance", fontsize=20)
ax.set_box_aspect([1.0, x1_size/x0_size, x2_size/x0_size])
#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()

　第2軸の値( x2 の要素)ごとに第0・1軸(x・y軸)平面における等高線を描画することで、楕円体を描画します。上手く立体的に見えない場合は、等高線の数( dist_levels や x2 の値の数)などを調整してください。

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

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

パラメータと形状の関係

　最後は、パラメータの値を少しずつ変化させて、距離の形状の変化をグラフで確認します。
　アニメーションの作図コードについては「Probability-Distribution/code/mahalanobis_distance/3d_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のノート】 - からっぽのしょこ」を参照してください。

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

参考文献

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

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

Amazon

おわりに

　本当にこれでいいのかよく分かっていませんが、いい感じのそれっぽいグラフができたので記事にしておきました。間違っていたり合っていたりしたら教えてください。
　3次元マハラノビス距離をグラフ化できたってことは、同様に3次元ガウス分布もグラフ化できるってことですよね？山の向きは逆ですが。それもやってみよう思ったのですが、Python版では1次元版も2次元版もまだ書いてなかった(とっくに書いてるつもりだった)ので、後回しにしてそれらとセットで書こうと思います。

　アレコレ苦労の末に3変数関数の等高線図を作れたので、作図処理は別の記事で解説することにしました(この記事より前に投稿するつもりだったのですが、まだ3行しか書けてません…)。他にもこの記事の補足やついでに何記事か書き終わったり書いてる途中だったりまだ書き始めてなかったりしてます。

　さてさてこの記事で、目次ページや日記的な記事も含めてですが、投稿数が700記事です！

　2024年8月16日は、BEYOOOOONDSとSeasoningSの小林萌花さんの24歳のお誕生日です♪おめでとうございます！

　アイドルであり音大卒のピアニストでもあります。私はピアノの良し悪しを判断できる感性は持ち合わせていませんが、ピアノの他にも絵が上手だったり、意外と雑だったりボブが超似合っていたり、アイドルとしては勿論のことアイドル以外の面でも語りつくせないほど素敵な方です。
　次こそはオーケストラコラボ公演を聴きに行かねばならないが次の機会はいつ頃でしょうか。

【次の内容】

からっぽのしょこ

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

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