草薙の研究ログ

英語教育関係。でも最近は統計(特にR)ネタが中心。

【R】マハラノビス距離のもとめ方と判別

Rでマハラノビス距離をもとめるためには,mahalanobis関数を使う。
勉強用に書いてみた。

マハラノビス距離自体は,やっぱり英語のwiki先生が詳しい。

Mahalanobis distance - Wikipedia, the free encyclopedia

#2変量のマハラノビス距離をもとめる
#数値例の作成
library(MASS)
round(mvrnorm(n=100,mu=c(50,50),Sigma=matrix(c(100,50,50,100),2,2)))->dat

#散布図による可視化
plot(dat1,xlab="x",ylab="y",pch=20,col=2,xlim=c(0,100),ylim=c(0,100))

#平均ベクトル(2次元上の中心)
m.dat<-colMeans(dat)

#それぞれの分散共分散行列
v.dat<-var(dat)

#マハラノビス距離の計算
mahalanobis(dat,m.dat,v.dat)

#距離に基づいて判別する
#2変量2グループの数値例を作成
round(mvrnorm(n=100,mu=c(50,50),Sigma=matrix(c(100,50,50,100),2,2)))->dat1
round(mvrnorm(n=100,mu=c(50,60),Sigma=matrix(c(100,50,50,100),2,2)))->dat2

#それぞれの平均ベクトル(2次元上の中心)
m.dat1<-colMeans(dat1)
m.dat2<-colMeans(dat2)

#それぞれの分散共分散行列
v.dat1<-var(dat1)
v.dat2<-var(dat2)

#マハラノビス距離の計算
mah11<-mahalanobis(dat1,m.dat1,v.dat1)
mah12<-mahalanobis(dat2,m.dat1,v.dat1)
mah21<-mahalanobis(dat1,m.dat2,v.dat2)
mah22<-mahalanobis(dat2,m.dat2,v.dat2)

#dat1について判別
res1<-ifelse(mah11<mah21,1,2)
#dat2について判別
res2<-ifelse(mah12>mah22,2,1)

#結果としてテーブルで出力
table.result<-table(data.frame(answer=c(rep(1,100),rep(2,100)),Result=c(res1,res2)))
table.result

#正判別率
(table.result[1,1]+table.result[2,2])/(nrow(dat1)+nrow(dat2))

#誤判別率
1-(table.result[1,1]+table.result[2,2])/(nrow(dat1)+nrow(dat2))

マハラノビス距離ってほんと捗る。効果量と等価だし。