目的

基準となる群の重心からのマハラノビスの距離を計算することにより，あるデータポイントがその基準群に帰属する確率を計算する
R には，mahalanobis 関数がある。使い方と出力の違いに注意のこと

使用法

Mahalanobis(dat, x)

引数

dat	基準群のデータ行列
x	所属確率を計算したいデータ行列

ソース

インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/mahalanobis.R", encoding="euc-jp")

# マハラノビスの距離による基準群への帰属確率
Mahalanobis <- function(dat,                 # 基準群のデータ行列
                        x)                      # 所属確率を計算するデータ行列
{
        dat <- subset(dat, complete.cases(dat))      # 欠損値を持つケースを除く
        n <- nrow(dat)                               # ケース数
        p <- ncol(dat)                               # 変数の個数
        ss <- var(dat)*(n-1)/n                       # 分散・共分散行列
        inv.ss <- solve(ss)                  # 分散共分散行列の逆行列
        m <- colMeans(dat)                   # 各変数の平均値
        dif <- t(t(x)-m)                     # 平均値からの偏差
        d2 <- apply(dif, 1,
                function(z) z %*% inv.ss %*% z) # マハラノビスの平方距離
        P <- pchisq(d2, p, lower.tail=FALSE) # 所属確率
        return(data.frame(d2=d2, P=P))
}


使用例

dat <- matrix(c(	# データ行列（4 変数，11 ケース）
	1, 2, 5, 3,
	7, 4, 8, 5,
	5, 4, 7, 1,
	2, 3, 5, 4,
	9, 5, 4, 7,
	2, 1, 4, 2,
	5, 4, 7, 4,
	2, 3, 5, 7,
	4, 1, 8, 2,
	3, 2, 5, 8,
	3, 5, 4, 8
), ncol=4, byrow=TRUE)

Mahalanobis(dat, rbind(c(1, 2, 5, 3), c(3, 1, 4, 1)))	# (1, 2, 5, 3) と (3, 1, 4, 1) から dat の重心へのマハラノビス距離を求める

出力結果例

> Mahalanobis(dat, rbind(c(1, 2, 5, 3), c(3, 1, 4, 1)))
        d2          P	# マハラノビス平方距離と，基準群への所属確率
1 2.000707 0.73562876
2 8.637195 0.07083601