マハラノビス距離
（2変量を扱うので少し難しくなります）.
マハラノビスの距離は、
集団においてデータが集団の中心からどのくらい離れているか.
を示す値です。（中心は集団の平均）
通常の距離は２点間の距離を表すのに対して、マハラノビスの距離は集団の中心からの距離であることが特徴です。中心に近ければ起こりやすく、中心から離れると起こりにくいと考えます。
2つの集団への帰属を判断する判別分析でも用いられます。
まず１次元の場合から考えましょう。
1 １次元.
1.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 と書きます。
\(D^2=z^2\) であるので、\(D^2\)
は自由度１のカイ二乗分布に従います。（大学レベル）
1.2 正規分布.
一般の正規分布では、尺度を標準偏差にすることで統一的な扱いが可能になります。
つまり中心から標準偏差いくつ離れているかを確率的距離と考えるわけです。
\(z=\dfrac{x-\mu}{\sigma}\).
標準化によって得られるzの絶対値がxのマハラノビス距離に相当します。
\(D^2=z^2\) は自由度１のカイ二乗分布に従います。
2 ２次元.
2.1 2変量標準正規分布.
相関係数 \(\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))
平均を中心として半径３の円内にほとんどの点が集まっているように見えますが、本当にそうでしょうか？
標準正規分布の場合、各成分は独立ですから、ユークリッド平方距離は自由度２のカイ二乗分布に従います。
- カイ二乗分布の上側５％の境界点は 5.991.
95%の点が ユークリッド平方距離 ＜ 5.991.
95%の点が ユークリッド距離 ＜ 2.364 なので、半径2.4 の円内に95%の点が集まっていることになります。
- カイ二乗分布の上側１％の境界点は 9.21.
99%の点が ユークリッド平方距離 ＜ 9.21.
99%の点が ユークリッド距離 ＜ 3.034 なので、半径3.0 の円内に99%の点が集まっていることになります。
このように、95%や99%の中心からの境界半径は、１次元の場合とは異なります。
距離は各成分の平方和であるのでカイ二乗分布に従うからです。
2.1.1 各成分の平方和の分布.
hist(x[,1]^2+x[,2]^2,prob=T)
curve(dchisq(x, 2),add=T)
成分の平方和の分布は、自由度２のカイ二乗分布と一致しています。
上側５％の境界点は 5.991 なので、平方和が 5.991 を超える確率は 0.05
となります。
sum(x[,1]^2+x[,2]^2>5.991)/10000## [1] 0.051データフレームを使う方法もあります。
成分の平方和 M を第３列に置き、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.0513 2変量正規分布.
3.1 相関のない場合.
成分の相関が０で分散が２であるとき、 各変数の分散が２なので\(\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))
3.2 相関のある場合.
各成分を標準化して平方和を計算しますが、相関があるので単純にはいきません。まず、一般の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}\).
相関係数の働きによって、
- y=x から近ければ、距離は \(2\rho z_xz_y\) 小さく計算され.
- y=x から遠ければ、距離は \(2\rho z_xz_y\) 大きく計算されます.
例えば、
\((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 上では近く配置しなければなりません。
3.2.1 \(D^2<1\) である点を薄青表示.
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.5778390dataE <- 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))
3.2.2 \(D^2>5.991\) である点を薄青表示.
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)
3.3 Rのマハラノビス関数.
草薙の研究ログ
様によれば、
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.5584210maha1 <- mahalanobis(x,c(0,0),sigma)
head(maha1)## [1] 1.4081261 0.4871487 5.7751150 3.7833466 0.8883745 0.57783904 例7_6.
関数を用いて計算する。
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] 20mahalanobis(c(6,10),c(4,6),sigma)## [1] 5mahalanobis(c(3,6),c(4,6),sigma)## [1] 1.5625mahalanobis(c(3,5),c(4,6),sigma)## [1] 1.015625x=6,y=2 のとき、M=20 であり、自由度２のカイ二乗分布では起こり得ない点なので、空白領域に出力されている。