Mahalanobis'距離と確率楕円の関係を書こうと思ったら、
思いの外、理論的背景が長くなったのでここで分けておきます。

Mahalanobis' Distance

点Xと群Aのマハラノビス距離は、下記で定義される。

D_{Mahalanobis}^2 = (X-\bar{A})^t cov(A) ^{-1}(X-\bar{A})　\\
\ \\
\begin{eqnarray}
cov(A):& &variance\ covariance\ matrix \\
\bar{A}:& &mean(A)\\
^t:& &transported\ matrix\\ 
^{-1}:& &inverse\ matrix \\
\end{eqnarray}

概念

2次元の場合を考える。
前回紹介した確率楕円の考え方を押さえておく。

二次元正規分布を例にして、中心からhoge%のプロットが含まれる範囲を描画してあります。
(RコードはAppendix-1)

ポイントは、同じ色の楕円上の点は同等に扱いたい、という事です。
このままでは、分布の中心(原点)からの距離が均等ではない。
そこで、補正をかけて下図のように変換したい。

2つの補正を行います。
・"長軸"をX軸、"短軸"をY軸に一致するように回転補正
・各軸方向に分散=1になるように補正

で、補正後のプロットと中心の距離のことを、Mahalanobis'距離と言います。

勘所の鋭いヒトは、これって主成分分析じゃん、と思うはずです。
その通り、でして、Mahalanobis'距離を求めるということは、分散共分散行列に基づく主成分分析をして、その主成分でプロットをスケールする、という事と同等の操作をしています。

数学的意味

投射先の軸上の分散

座標(X, Y)をθ度回転させた座標(X',Y')は、

\begin{eqnarray}
X' &=& &cosθ*X &+ sinθ*Y \\
Y' &=& -&sinθ*X &+ cosθ*Y
\end{eqnarray}

つまり回転したい向きの方向ベクトル(cosθ, sinθ)が分かれば良い。
どうやってθを決定すれば良いかというと、分布の"長軸"なわけなので、
その軸に投射した時に分散が最大になる軸を選べば良い。

少し記号を変えて、投射前座標(X1,X2)を方向ベクトル(h1,h2)の軸に投射した座標をYとします。

Y=h_1*X_1+h_2*X_2

軸上の座標の標本平均は、

\begin{eqnarray}
\bar{Y}&=&\frac{1}{n}\sum_{i=1}^n{Y_i}=\frac{1}{n}\sum_{i=1}^n{(h_1X_1+h_2X_2)} \\
&=&\frac{1}{n}\sum_{i=1}^n{h_1X_1}+\frac{1}{n}\sum_{i=1}^n{h_2X_2} \\
&=&h_1\bar{X_1}+h_2\bar{X_2}
\end{eqnarray}

これを用いて、軸(h1,h2)上の標本分散は、

\begin{eqnarray}
s^2_{(h_1, h_2)}
&=&\frac{1}{n-1}\sum_{i=1}^n{(Y_i-\bar{Y})^2} \\
&=&\frac{1}{n-1}\sum_{i=1}^n{ \Big\{(h_1X_{i1}+h_2X_{i2})-(h_1\bar{X_1}+h_2\bar{X_2}) \Big\}^2 } \\
&=&h_1^2 \frac{1}{n-1}\sum_{i=1}^n{(X_{i1}-\bar{X_1})^2}
+2h_1h_2\frac{1}{n-1}\sum_{i=1}^n{(X_{i1}-\bar{X_1})(X_{i2}-\bar{X_2})}
+h_2^2 \frac{1}{n-1}\sum_{i=1}^n{(X_{i2}-\bar{X_2})^2} \\
&=& h_1^2\ s_{11}+2h_1h_2\ s_{12}+h_2^2\ s_{22}
\end{eqnarray}

え〜と、初項から、X1の標本分散s11、X1,X2の標本共分散s12、X2の標本分散s22の形になっているという事です。

つーことで、

S= \left( 
\begin{array}{cc}
s_{11} & s_{12} \\ s_{12} & s_{22}
\end{array} \right), 
h= \left( \begin{array}{c} h_1 \\ h_2 \end{array} \right)

分散共分散行列Sと方向ベクトルhを用いて、

s^2_{(h_1, h_2)} = h^t Sh

このhが主成分の固有ベクトルだったとすると、

\begin{equation}
(S-λE)h=0
　\Rightarrow　h^t Sh = λ = s^2_{(h_1, h_2)}
\end{equation}

なる固有値λが存在する(Eは単位行列)。
すなわち、主成分軸方向の分散は、対応する固有値に一致する。

（あとで、ここに戻ってきます）

Mahalanobis'距離

で、求めたい値は、
1. 元座標(X1, X2)を固有ベクトルh1(h11, h12)とh2(h21,h22)を使って主成分座標(pX1, pX2)に変換
2. 主成分座標(pX1, pX2)と中心の距離を固有値λ(λ1,λ2)のルートで割った距離を算出
という事になる。
元座標(X1, X2)は、中心原点(0,0)に平行移動しておくものとします。(回転の関係で)

さて、2から攻めると、

\begin{eqnarray}
D^2_{mahalanobis} 
&=& \frac{pX_1^2+pX_2^2}{s^2}\\
&=& \frac{pX_1^2+pX_2^2}{λ}\\
&=& ( \begin{array}{cc}pX_1 & pX_2 \end{array} ) 
\Big( \begin{array}{c} \frac{pX_1}{λ_1} \\ \frac{pX_2}{λ_2} \end{array} \Big)  \\
&=& ( \begin{array}{cc}pX_1 & pX_2 \end{array} ) 
\Big( \begin{array}{cc} \frac{1}{λ_1} & 0\\ 0 & \frac{1}{λ_2} \end{array} \Big)   ( \begin{array}{cc}pX_1 & pX_2 \end{array} )^t 
\\
\end{eqnarray}

1を盛り込む。

pX_1 = h_{11} X_1+h_{12}X_2\\
pX_2 = h_{21} X_1+h_{22}X_2

より、

\begin{eqnarray}
( \begin{array}{cc}pX_1 & pX_2 \end{array} ) &=&
( \begin{array}{cc}X_1 & X_2 \end{array} )
\Big( \begin{array}{cc} h_{11} & h_{21}\\ h_{12} & h_{22} \end{array} \Big) \\
& \equiv &
( \begin{array}{cc}X_1 & X_2 \end{array} )H
\end{eqnarray}

Hは固有ベクトルマトリックス。

また、

( \begin{array}{cc}pX_1 & pX_2 \end{array} )^t =
H^t
( \begin{array}{cc}X_1 & X_2 \end{array} )^t \\

これらを用いて、

\begin{eqnarray}
D^2_{mahalanobis} 
&=& ( \begin{array}{cc}pX_1 & pX_2 \end{array} ) 
\Big( \begin{array}{cc} \frac{1}{λ_1} & 0\\ 0 & \frac{1}{λ_2} \end{array} \Big)   ( \begin{array}{cc}pX_1 & pX_2 \end{array} )^t 
\\
&=& ( \begin{array}{cc}X_1 & X_2 \end{array} )
H
\Big( \begin{array}{cc} \frac{1}{λ_1} & 0\\ 0 & \frac{1}{λ_2} \end{array} \Big)
H^t
( \begin{array}{cc}X_1 & X_2 \end{array} )^t

\end{eqnarray}

大丈夫ですか？
ページトップの数式と見比べてください。大分近づいてきました。
あとは、

H
\Big( \begin{array}{cc} \frac{1}{λ_1} & 0\\ 0 & \frac{1}{λ_2} \end{array} \Big)
H^t
 \equiv S^{-1}

が言えればOKですね？

少しトリッキーですが、

\begin{eqnarray}
\Big( \begin{array}{cc}
 \frac{1}{λ_1} & 0\\ 0 & \frac{1}{λ_2} 
\end{array} \Big)
&=&
\Big( \begin{array}{cc}
 \frac{λ_2}{λ_1λ_2} & 0\\ 0 & \frac{λ_1}{λ_1λ_2} 
\end{array} \Big) \\
&=&
\frac{1}{λ_1λ_2}
\Big( \begin{array}{cc} λ_2 & 0\\ 0 & λ_1 \end{array} \Big) \\
&=&
\Big( \begin{array}{cc} λ_1 & 0\\ 0 & λ_2 \end{array} \Big)^{-1} \\
&\equiv& \Lambda^{-1}
\end{eqnarray}

と置く。
λ大文字は、対角成分が固有値の行列です。

H \Lambda^{-1} H^t \equiv S^{-1}

が言えればOK。Sは分散共分散行列。

ふ〜。さて、(あとでここに)と言っていた式を再掲。

(S-λE)h=0

これを今風の記号に書き直すと、

(S-λ_1 E)\Big(\begin{array}{c}h_{11}\\ h_{12}\end{array}\Big)=0 \\
(S-λ_2 E)\Big(\begin{array}{c}h_{21}\\ h_{22}\end{array}\Big)=0

これを少し整理して、

S\Big(\begin{array}{c}h_{11}\\ h_{12}\end{array}\Big) = λ_1\Big(\begin{array}{c}h_{11}\\ h_{12}\end{array}\Big)\\
S\Big(\begin{array}{c}h_{21}\\ h_{22}\end{array}\Big)=
λ_2\Big(\begin{array}{c}h_{21}\\ h_{22}\end{array}\Big)

何がしたいかというと、

\begin{eqnarray}
S \Big(\begin{array}{cc}h_{11} & h_{21} \\ h_{12} & h_{22} \end{array} \Big)&=&
\Big(\begin{array}{cc}λ_{1}h_{11} & λ_2h_{21} \\ λ_1h_{12} & λ_2h_{22} \end{array} \Big) \\
&=&\Big(\begin{array}{cc}h_{11} & h_{21} \\ h_{12} & h_{22} \end{array} \Big)
\Big(\begin{array}{cc}λ_1 & 0 \\ 0 & λ_2 \end{array} \Big)
\end{eqnarray}

すなわち、

SH=H \Lambda

これを用いて、

\begin{eqnarray}
SH&=&H \Lambda \\
S^{-1}SH&=&S^{-1}H\Lambda \\
H&=&S^{-1}H\Lambda \\
H\Lambda^{-1}&=&S^{-1}H\Lambda \Lambda^{-1}=S^{-1}H\\
H\Lambda^{-1}H^{-1}&=&S^{-1}HH^{-1}=S^{-1}
\end{eqnarray}

おっと、ここで、

\begin{eqnarray}
H^tH&=&
\Big( \begin{array}{cc} h_{11} & h_{12} \\ h_{21} & h_{22} \end{array} \Big)
\Big( \begin{array}{cc} h_{11} & h_{21} \\ h_{12} & h_{22} \end{array} \Big)\\
&=& \Big( \begin{array}{cc} h_{11}^2+h_{12}^2 & h_{11}h_{21}+h_{12}h_{22} \\
 h_{11}h_{21}+h_{12}h_{22} & h_{21}^2+h_{22}^2 \end{array} \Big) \\
&=&\Big( \begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array} \Big) \\

\therefore H^t&=&H^{-1}
\end{eqnarray}

対角成分＝1は、単位ベクトルだから。
それ以外の成分=0は、直行ベクトルだから。

やっとなっ、

H\Lambda^{-1}H^{t}=S^{-1}

従って、

\begin{eqnarray}
D^2_{mahalanobis} 
&=& ( \begin{array}{cc}X_1 & X_2 \end{array} )
H \Lambda^{-1} H^t
( \begin{array}{cc}X_1 & X_2 \end{array} )^t\\
&=& ( \begin{array}{cc}X_1 & X_2 \end{array} )
S^{-1}
( \begin{array}{cc}X_1 & X_2 \end{array} )^t
\end{eqnarray}

っという事になりました。

Appendix

Appendix-1. 確率楕円

詳しくは前回

library("ellipse")
library("ggplot2")

# "回転前"のデータ
set.seed(13)
dat0 <- matrix(rnorm(1000), , 2)

# 回転
sita <- pi/4
rot <- matrix(c(cos(sita), -sin(sita), sin(sita), cos(sita)),2)
dat <- data.frame(dat0 * c(2, 1) %*% rot)

# 楕円を描く確率水準
level <- seq(0.1, 0.9, by=0.1)

# 楕円座標の算出
datE <- NULL
for(i in level){
  datE <-rbind(datE, cbind(i, ellipse(cov(dat), centre=apply(dat, 2, mean), level=i)))
}
datE <- data.frame(datE)
datE$i <- factor(datE$i)

# カラーコードの設定
Col <- colorRamp(c("Red","skyblue"))
color <- rgb(Col(seq(0, 1, length=length(level)))/255)
# 軸範囲
ax <- max(abs(dat))

# 描画
ggplot()+
  layer(
    data = dat,
    mapping = aes(x=X1, y=X2),
    geom = "point",
    stat = "identity",
    position = "identity",
    params = list(color="grey")
  )+
  layer(
    data = datE,
    mapping = aes(x=X1, y=X2, color=i),
    geom = "path",
    stat = "identity",
    position = "identity",
    params = list(size=1.2)
  )+
  xlim(-ax, ax)+
  ylim(-ax, ax)+
  scale_color_manual(values=color)+
  theme_bw()

っとまぁ、最近はggplot2を使うようにしていますが、普通にplotを使っても書けます。

# ggplot2を使わない描画
par(mar=c(2,2,1,1))
plot(0,0,type="n", xlim=c(-ax,ax), ylim=c(-ax,ax), mgp=c(2, 0.3, 0), tcl=-0.2)
  points(dat, col="grey", pch=16)
for(i in 1:length(level)){
  lines(datE[datE$i==level[i], 2:3], col=color[i], lwd=2)
}

ん〜。どちらが見直した時に読みやすいか、の好みで。
plotも、結構使い勝手いいと思いますよっ。

参考文献

静岡大金久保先生のページより、判別分析(マハラノビス)
書籍ですが、『多次元データ解析法(Rで学ぶデータサイエンス)』