マハラノビスの距離で判別分析

今回は、あらかじめ合否のデータリストがある状態で、新たなデータを入れたときにそのデータの合否がそれぞれの要因によってどう分類されるか？を調べます。

そのための変数は、比重、色素の二種類をもつデータです。

> data = read.csv("7-1.csv")
> data
   hijyu sikiso gouhi
1  0.815     53     0
2  0.851     53     1
3  0.831     44     0
4  0.833    160     1
5  0.854     17     1
6  0.844     33     0
7  0.834     58     0
8  0.827     64     0
9  0.833     28     0
10 0.839     48     0
11 0.830     30     0
12 0.836     76     1
13 0.819     48     0
14 0.843     29     0
15 0.839    111     1
16 0.828     48     0
17 0.832     92     1
18 0.856     46     1
19 0.850    213     1
20 0.844     53     1

mahalanobisの関数は群馬大学で公開されていたので、それをインストールします

source("http://aoki2.si.gunma-u.ac.jp/R/src/mahalanobis.R", encoding="euc-jp")

合格・不合格のデータをそれぞれ分類しておきます。

data.goukaku = data[data$gouhi > 0,]
data.hugoukaku = data[data$gouhi == 0,]
okmatrix = matrix(c(data.goukaku$hijyu,data.goukaku$sikiso), ncol=2)
ngmatrix = matrix(c(data.hugoukaku$hijyu, data.hugoukaku$sikiso), ncol=2)

それぞれのデータを合否をもとに散布図に描きます。

plot(data.goukaku$hijyu, data.goukaku$sikiso, col="red")
par(new=T)
plot(data.hugoukaku$hijyu, data.hugoukaku$sikiso, col="blue")

これらのデータに対してそれぞれマハラノビスの距離を求めてみると以下のようになります。

Mahalanobis(
    okmatrix, 
    matrix(c(data$hijyu, data$sikiso), ncol=2)
)
Mahalanobis(
    ngmatrix, 
    matrix(c(data$hijyu, data$sikiso), ncol=2)
)

> Mahalanobis(
+     okmatrix, 
+     matrix(c(data$hijyu, data$sikiso), ncol=2)
+ )
           d2            P
1  14.6234953 0.0006676492
2   0.8262871 0.6615673028
3   4.1478479 0.1256916070
4   2.2011573 0.3326785313
5   2.2138264 0.3305778156
6   1.1171826 0.5720143069
7   2.3268822 0.3124093037
8   5.2037846 0.0741331648
9   4.1470580 0.1257412584
10  1.2996507 0.5221369561
11  5.4203512 0.0665251223
12  1.1888033 0.5518926786
13 11.5050621 0.0031747352
14  1.3830584 0.5008096498
15  0.3394467 0.8438982336
16  5.4350812 0.0660369650
17  2.0956566 0.3506985433
18  2.0345910 0.3615714803
19  6.6208690 0.0365003118
20  0.4793626 0.7868786014
> Mahalanobis(
+     ngmatrix, 
+     matrix(c(data$hijyu, data$sikiso), ncol=2)
+ )
             d2            P
1  3.563685e+00 1.683277e-01
2  1.012448e+01 6.331349e-03
3  4.555432e-04 9.997723e-01
4  1.315911e+02 2.662953e-29
5  8.285814e+00 1.587663e-02
6  2.264394e+00 3.223244e-01
7  2.552024e+00 2.791483e-01
8  3.085350e+00 2.138084e-01
9  2.104592e+00 3.491352e-01
10 1.663588e+00 4.352677e-01
11 2.086140e+00 3.523713e-01
12 1.225280e+01 2.184435e-03
13 2.164813e+00 3.387792e-01
14 2.342371e+00 3.099993e-01
15 5.088193e+01 8.935729e-12
16 1.725864e-01 9.173252e-01
17 2.262384e+01 1.222634e-05
18 1.170327e+01 2.875192e-03
19 3.206605e+02 2.341328e-70
20 5.226196e+00 7.330710e-02

ここに新たなデータとして、比重=0.837, 色素=65のデータを入れてみます。

test = c(0.837, 65)

> test = c(0.837, 65)
> Mahalanobis(
+     okmatrix, 
+     matrix(test, ncol=2)
+ )
        d2         P
1 1.212049 0.5455152
> 
> Mahalanobis(
+     ngmatrix, 
+     matrix(test, ncol=2)
+ )
        d2          P
1 6.439559 0.03996387

すると、ngmatrixと比較したときの方が距離が短くなるので、結果、このデータは不合格のグループに属する可能性が高い、と判定できることになります。