(2変量を扱うので少し難しくなります).
マハラノビスの距離は、
集団においてデータが集団の中心からどのくらい離れているか.
を示す値です。(中心は集団の平均)
通常の距離は2点間の距離を表すのに対して、マハラノビスの距離は集団の中心からの距離であることが特徴です。中心に近ければ起こりやすく、中心から離れると起こりにくいと考えます。
2つの集団への帰属を判断する判別分析でも用いられます。
まず1次元の場合から考えましょう。
平均0, 分散1の確率変数の理想的な分布
乱数を1000個出力しヒストグラムを描きます。
x <- rnorm(1000,0,1)
hist(x,prob=T,breaks = "Scott")
curve(dnorm(x,0,1),add=T)
曲線はガウスの発見した理論曲線で、次の式で表されます。
\(\phi \left( z\right)
=\dfrac{1}{\sqrt{2\pi }}\exp \left( -\dfrac{z^{2}}{2}\right)\sim
N(0,1)\).
中心 0 から距離 \(|z|\)
にあるデータが出力される確率は \(\phi(z)\) .
中心極限定理では、一般の確率変数でもサンプル数が大きくなれば、標本平均を標準化した確率変数は標準正規分布に近くなることが知られています。
確率は中心からの距離によって決まります。(中心=平均のこと)
\(D^2=z^2\) であるので、\(D^2\)
は自由度1のカイ二乗分布に従います。(大学レベル)
一般の正規分布では、尺度を標準偏差にすることで統一的な扱いが可能になります。
つまり中心から標準偏差いくつ離れているかを確率的距離と考えるわけです。
\(z=\dfrac{x-\mu}{\sigma}\).
標準化によって得られるzの絶対値がxのマハラノビス距離に相当します。
\(D^2=z^2\) は自由度1のカイ二乗分布に従います。
相関係数 \(\rho=0\) 、x ,y は独立であり散布図は円状になります。
library(mvtnorm)
sigma <- matrix(c(1,0,0,1), ncol=2)
x <- rmvnorm(n=10000, mean=c(0,0), sigma=sigma)
plot(x,asp=1,col="grey",xlim=c(-3,3),ylim=c(-3,3))
平均を中心として半径3の円内にほとんどの点が集まっているように見えますが、本当にそうでしょうか?
標準正規分布の場合、各成分は独立ですから、ユークリッド平方距離は自由度2のカイ二乗分布に従います。
このように、95%や99%の中心からの境界半径は、1次元の場合とは異なります。
距離は各成分の平方和であるのでカイ二乗分布に従うからです。
hist(x[,1]^2+x[,2]^2,prob=T)
curve(dchisq(x, 2),add=T)
成分の平方和の分布は、自由度2のカイ二乗分布と一致しています。
上側5%の境界点は 5.991 なので、平方和が 5.991 を超える確率は 0.05
となります。
sum(x[,1]^2+x[,2]^2>5.991)/10000
## [1] 0.051
データフレームを使う方法もあります。
成分の平方和 M を第3列に置き、Mの値が 5.991 を超える行を
nrow関数で数えます。
data <- data.frame(X=x[,1],Y=x[,2])
plot(x,asp=1,col="grey",xlim=c(-3,3),ylim=c(-3,3))
dataE <- data.frame(X=x[,1],Y=x[,2],M=x[,1]^2+x[,2]^2)
head(dataE)
X | Y | M |
---|---|---|
0.3935669 | -0.0431349 | 0.1567555 |
-1.5520218 | 0.2375760 | 2.4652142 |
-1.0286571 | -0.7677352 | 1.6475527 |
0.2533972 | 1.6495473 | 2.7852166 |
-1.2658916 | 0.1747631 | 1.6330237 |
-1.7310656 | 0.1114805 | 3.0090161 |
par(new=T)
plot(dataE[dataE $ M >5.991,][,-3],asp=1,col="lightblue",xlim=c(-3,3),ylim=c(-3,3))
#nrow:データフレームの行数をカウント
nrow(dataE[dataE $ M >5.991,])/10000
## [1] 0.051
成分の相関が0で分散が2であるとき、 各変数の分散が2なので\(\sqrt{2}\)倍に広がっています。
library(mvtnorm)
sigma <- matrix(c(2,0,0,2), ncol=2)
x <- rmvnorm(n=10000, mean=c(0,0), sigma=sigma)
data <- data.frame(X=x[,1],Y=x[,2])
plot(x,asp=1,col="grey",xlim=c(-3,3),ylim=c(-3,3))
各成分を標準偏差で割って標準化することで、2変量標準正規分布に変換することができます。
plot(x/sqrt(2),asp=1,col="grey",xlim=c(-3,3),ylim=c(-3,3))
各成分を標準化して平方和を計算しますが、相関があるので単純にはいきません。まず、一般の2変量正規分布について確認しておきます。(高校レベルを超える).
2変量正規分布
\(\begin{pmatrix} x \\ y \end{pmatrix} \sim
N_2\left( \begin{pmatrix} \mu_x \\ \mu_y \end{pmatrix},\begin{pmatrix}
\sigma_x^2 & \sigma_{xy}\\ \sigma_{xy} & \sigma_y^2
\end{pmatrix}\right) = N_2\left( \begin{pmatrix} \mu_x \\ \mu_y
\end{pmatrix},\begin{pmatrix} \sigma_x^2 & \rho\sigma_x \sigma_y \\
\rho\sigma_x \sigma_y & \sigma_y^2
\end{pmatrix}\right)\) のとき、
\(f(x,y)=\dfrac{1}{\left(\sqrt{2\pi}\right)^2\sqrt{|\Sigma|}}\exp\left[{-\dfrac{1}{2}(X-\mu)'\
\Sigma^{-1}(X-\mu)}\right]\)
\(\Sigma^{-1}= \dfrac{1}{\sigma_x^2\sigma_y^2(1-\rho^2)} \begin{pmatrix} \sigma_y^2 & -\sigma_{xy}\\ -\sigma_{xy} & \sigma_x^2 \end{pmatrix}\)
マハラノビスの平方距離
\(\begin{split} D^2&=(X-\mu)'\ \Sigma^{-1}(X-\mu)\\ &= \dfrac{1}{\sigma_x^2\sigma_y^2(1-\rho^2)} \left[\sigma_y^2(x-\mu_x)^2-2\sigma_{xy}(x-\mu_x)(y-\mu_y)+\sigma_x^2(y-\mu_y)^2\right]\\ &=\dfrac{1}{(1-\rho^2)} \left[\dfrac{(x-\mu_x)^2}{\sigma_x^2}-2\rho\dfrac{(x-\mu_x)}{\sigma_x}\dfrac{(y-\mu_y)}{\sigma_y}+\dfrac{(y-\mu_y)^2}{\sigma_y^2}\right] \\ &=\dfrac{1}{(1-\rho^2)}\left[z_x^2-2\rho z_xz_y+z_y^2\right] \end{split}\).
相関係数の働きによって、
例えば、
\((z_x,z_y)=(1,1)\)) なら\(D^2=1/0.75\times (1-1+1)=1/0.75=0.33\).
\((z_x,z_y)=(-1,1)\)) なら\(D^2=1/0.75\times (1+1+1)=1/0.75\times
3=4.0\).
\(D^2\)を等距離にするには、y=x 上では遠く配置し、y=-x 上では近く配置しなければなりません。
sigma <- matrix(c(1,0.5,0.5,1), ncol=2)
x <- rmvnorm(n=10000, mean=c(0,0), sigma=sigma)
M <-(1/(1-0.5^2))*(x[,1]^2+x[,2]^2-2*0.5*x[,1]*x[,2])
data <- data.frame(X=x[,1],Y=x[,2])
head(M)
## [1] 1.4081261 0.4871487 5.7751150 3.7833466 0.8883745 0.5778390
dataE <- data.frame(X=x[,1],Y=x[,2],M=(1/(1-0.5^2))*(x[,1]^2+x[,2]^2-2*0.5*x[,1]*x[,2]))
head(dataE)
X | Y | M |
---|---|---|
-0.7340217 | 0.4404567 | 1.4081261 |
-0.1162716 | -0.6541410 | 0.4871487 |
-0.2572745 | -2.1978626 | 5.7751150 |
-0.5033194 | -1.8787774 | 3.7833466 |
0.9059939 | 0.2279144 | 0.8883745 |
0.7425889 | 0.2305799 | 0.5778390 |
plot(x,asp=1,col="grey",xlim=c(-3,3),ylim=c(-3,3))
par(new=T)
plot(dataE[dataE $ M <1,][,-3],asp=1,col="lightblue",xlim=c(-3,3),ylim=c(-3,3))
plot(x,asp=1,col="grey",xlim=c(-3,3),ylim=c(-3,3))
dataE <- data.frame(X=x[,1],Y=x[,2],M=(1/(1-0.5^2))*(x[,1]^2+x[,2]^2-2*0.5*x[,1]*x[,2]))
#head(dataE)
par(new=T)
plot(dataE[dataE $ M >5.991,][,-3],asp=1,col="lightblue",xlim=c(-3,3),ylim=c(-3,3))
カイ二乗分布の上側5%の境界点が 5.991 であることは次のように確かめられます。
nrow(dataE[dataE $ M >5.991,])/10000
## [1] 0.0486
カイ二乗分布に従うことは、次のように確かめられます。
hist(M,prob=T)
curve(dchisq(x, 2),add=T)
草薙の研究ログ
様によれば、
Rにはマハラノビスの平方距離を算出するmahalanobis()関数がある。
標本平均と標本分散共分散行列を用いた場合 maha と
母平均と母分散共分散行列を用いた場合 maha1 で若干値が異なっている。
#平均ベクトル(2次元上の中心)
m.dat<-colMeans(x)
#それぞれの分散共分散行列
v.dat<-var(x)
#マハラノビス距離の計算
maha <- mahalanobis(x,m.dat,v.dat)
head(maha)
## [1] 1.4461143 0.5005263 5.9072865 3.8859309 0.8652062 0.5584210
maha1 <- mahalanobis(x,c(0,0),sigma)
head(maha1)
## [1] 1.4081261 0.4871487 5.7751150 3.7833466 0.8883745 0.5778390
関数を用いて計算する。
sigma <- matrix(c(1,1.2,1.2,4), ncol=2)
x <- rmvnorm(n=10000, mean=c(4,6), sigma=sigma)
M <- mahalanobis(x,c(4,6),sigma)
dataE <- data.frame(X=x[,1],Y=x[,2],M=M)
plot(x,col="grey",xlim=c(1,7),ylim=c(-1,12))
data <- data.frame(X=x[,1],Y=x[,2])
par(new=T)
plot(dataE[dataE $ M <1,][,-3],col="lightblue",xlim=c(1,7),ylim=c(-1,12))
par(new=T)
plot(x=6,y=2,col="red",xlim=c(1,7),ylim=c(-1,12))
par(new=T)
plot(x=6,y=10,col="red",xlim=c(1,7),ylim=c(-1,12))
par(new=T)
plot(x=3,y=6,col="red",xlim=c(1,7),ylim=c(-1,12))
mahalanobis(c(6,2),c(4,6),sigma)
## [1] 20
mahalanobis(c(6,10),c(4,6),sigma)
## [1] 5
mahalanobis(c(3,6),c(4,6),sigma)
## [1] 1.5625
mahalanobis(c(3,5),c(4,6),sigma)
## [1] 1.015625
x=6,y=2 のとき、M=20 であり、自由度2のカイ二乗分布では起こり得ない点なので、空白領域に出力されている。