「統計学関連なんでもあり」の過去ログ--- 000

No.00606　判別分析を行う関数disc()の返す値　　【波音】　2006/07/09(Sun) 16:00

http://aoki2.si.gunma-u.ac.jp/R/disc.htmlにあるdisc()を用いて,

dat<-matrix(c(3.2,3.5,3.0,3.6,4.9,2.9,1.8,2.5,1.7,3.6,6.5,7.5,7.6,8.5,8.1,6.9,5.9,6.8,
8.0,7.1,7.1,5.7,7.0,6.6,5.5,4.8,4.6,4.1,4.3,2.6,3.8,4.7,2.3,2.5,2.4,1.3),ncol=2)
group<-matrix(c(rep(1,10),rep(2,8)),ncol=1)
res<-disc(dat,group)
res

とすると,マハラノビスの距離の2乗が次のように計算されます。

$distance
         D to g1    D to g2
 [1,]  2.7419508 58.9097404
(以下省略)

しかし自力(というかこの関数を用いずに)で以下のようにして計算すると,数値がどうしても合わないのですが,どうしてなのでしょうか。

x<-dat[1:10,1:2]
y<-dat[11:18,1:2]
gx<-solve(cov(x))
gy<-solve(cov(y))
dim1<-matrix(c(3.2-mean(x[,1]),8.0-mean(x[,2])),ncol=2)
dim2<-matrix(c(3.2-mean(y[,1]),8.0-mean(y[,2])),ncol=2)
dim1%*%gx%*%t(dim1) #グループ1の中心までのマハラノビスの距離の2乗
dim2%*%gy%*%t(dim2) #グループ2の中心までのマハラノビスの距離の2乗

disc()を使った場合||自力でやった場合
------------------------------------
    2.7419508     ||   2.533912
   58.9097404     ||  63.18003

No.00607　Re: 判別分析を行う関数disc()の返す値　　【青木繁伸】　06/07/09(Sun) 19:34

波音さんがやったのは，各群の分散・共分散行列が等しいことを仮定しない方法（二次の判別分析）です。

disc は，線形判別分析で，各群の分散・共分散行列が等しいことを仮定し，その推定値が変動の和を（サンプルサイズー群の数）で割ったものです。
> x<-dat[1:10,1:2]
> y<-dat[11:18,1:2]
> gxy <- solve(( cov(x)*9+cov(y)*7)/16) # 共通の分散・共分散行列の逆行列
> dim1<-matrix(c(3.2-mean(x[,1]),8.0-mean(x[,2])),ncol=2)
> dim2<-matrix(c(3.2-mean(y[,1]),8.0-mean(y[,2])),ncol=2)
> dim1%*%gxy%*%t(dim1) #グループ1の中心までのマハラノビスの距離の2乗
         [,1]
[1,] 2.741951
> dim2%*%gxy%*%t(dim2) #グループ2の中心までのマハラノビスの距離の2乗
         [,1]
[1,] 58.90974

No.00608　Re: 判別分析を行う関数disc()の返す値　　【波音】　06/07/09(Sun) 23:15

回答ありがとうございます。

なぜ私がやった場合とdisc()が返してくれる結果が違うのかは分かりました。

ところで,誤判別の確率を求めるための関数を定義したいのですが,3変数以上になった場合には分散・共分散行列はどのようにして求めればよいのでしょうか。
gohan<-function(v1,v2){
	mu1<-apply(v1,2,mean)
	mu2<-apply(v2,2,mean)
	n1<-dim(v1)[1]
	n2<-dim(v2)[1]
	n3<-dim(v1)[2]
	pool.dim<-(cov(v1)*(n1-1)+cov(v2)*(n2-1))/(n1+n2-n3) #プール後の分散・共分散行列
	dim1<-matrix(mu1-mu2,ncol=dim(v1)[2]) #グループ1の列平均からグループ2の列平均を引いたもの
	cd<-solve(pool.dim) #プール後の逆行列
	cs<-(dim1%*%cd)%*%t(dim1) #誤判別の確率
	list("プール後の分散・共分散行列"=pool.dim,"プール後の逆行列"=cd,"誤判別の確率"=cs)
}
2変数の場合はうまくいったのですが,3変数以上になるとpool.dim<-(cov(v1)*(n1-1)+cov(v2)*(n2-1))/(n1+n2-n3)の分母がこれではよろしくないと思うのですが,色々やっているうちに混乱してきてしまいました。3変数以上の場合のプール後の分散・共分散行列の計算方法を教えてください。

No.00610　Re: 判別分析を行う関数disc()の返す値　　【青木繁伸】　06/07/10(Mon) 10:12
w <- matrix(apply(sapply(1:ng, function(i) var(data[group == g.name[i],])*(num[i]-1)), 1, sum), ncol=p)
の部分が，変動・共変動行列を表すことになっておりますので，これを自由度＝「全ケース数　－　グループ数」で割ってやればプールした分散・共分散行列になります。

No.00611　Re: 判別分析を行う関数disc()の返す値　　【波音】　06/07/10(Mon) 12:22

度々の質問で申し訳ありませんが,sapply()の部分は各グループの分散・共分散行列を作って,その各要素に各グループのケース数(行数)から1を引いたnum[i]-1を掛けているわけですよね。
しかしながら,1:ngのngには何が入るのかがよく分からないのです。ここには変数の数が入るのでしょうか？

No.00612　Re: 判別分析を行う関数disc()の返す値　　【青木繁伸】　06/07/10(Mon) 12:46
ng <- length(num <- table(group))
で，グループを表す変数 group の度数分布をとって，その長さ（変数値の種類の数）を計算しています。

グループの数ですね。

No.00613　Re: 判別分析を行う関数disc()の返す値　　【波音】　06/07/10(Mon) 15:17

ありがとうございます。おかげさまでプール後の分散・共分散行列までクリアすることができましたが,あと1つだけご助言願いたいのですが。

最後に誤判別の確率を求めるときに,各グループの平均を引くとき,2グループならグループ1の列平均-グループ2の列平均となりますが,3つ以上の場合がまたもや分からずに困っています。だいぶ頭をひねっても答えが浮かばないのですがよろしくお願いします。
mean <- matrix(sapply(1:p, function(i) tapply(dat[,i], group, mean)), ncol=p) #グループごとの列平均
dim1<-matrix(mean[1,]-mean[2,],ncol=p)　#この部分がいかにしても分かりません
cs<-(dim1%*%cd)%*%t(dim1) #誤判別の確率

No.00614　Re: 判別分析を行う関数disc()の返す値　　【青木繁伸】　06/07/10(Mon) 16:27

誤判別率は，二群の組み合わせごとに計算されます

1群と2群の誤判別率，1群と3群の誤判別率，2群と3群の誤判別率というように。

No.00615　Re: 判別分析を行う関数disc()の返す値　　【波音】　06/07/10(Mon) 18:06

> 誤判別率は，二群の組み合わせごとに計算されます

そうでしたか,プログラムを書くことに夢中になってしまってついウッカリとしてしまいました。度々の回答,本当にありがとうございました。