距离判别法与R程序实战
设mu1,mu2,cov1,cov2
分别是总体G1，G2的均值向量和协方差矩阵。一样品x到各总体的马氏距离为
D(x,G1)=mahalanobis(x,mu1,cov1),
D(x,G2)=mahalanobis(x,mu2,cov2)）
一，距离判别法：
令
g(x)=D(x,G1)-D(x,G2)
若
(i) g(x)<0,则x属于G1;
(ii) g(x)>0,则x属于G2;
(iii) g(x)=0,则无法判断。
二，距离线性判别法：
假设总体G1，G2的协方差矩阵相同，即cov1=cov2=cov。
令
{D(x,G2)}^2-{D(x,G1)}^2=2w(x)，
w(x)=t(x)a-t(mu_bar)a,称之为线性判别函数。
其中a=solve(cov)(mu1-mu2),mu_bar=(mu1+mu2)/2，t(x)表示矩阵x的转置，solve(x)表示计算矩阵x的逆矩阵。
若
(i)w(x)>0,则x属于G1;
(ii)w(x)<0,则x属于G2;
(iii)w(x)=0,则无法判断。
三、R语言实现
在R语言中，假设样本资料矩阵为x，每一列表示各变量的样本，每一行表示一个案例。那么向量的样本矩阵的计算为apply(x,2,mean),变量之间的协方差矩阵为cov(x)。因此若x1，x2分别为总体G1，G2的样本，x1的样本容量为n1，x2的样本容量为n2，那么x1，x2的样本向量均值估计为
mu1=apply(x1,2,mean),mu2=apply(x2,2,mean),
若G1，G2的协方差矩阵相同，即
cov1=cov2=cov，那么
协方差矩阵的估计为
cov=((n1-1)cov(x1)+(n1-1)cov(x1))/(n1+n2-2)
例子1，
人文发展指数是联合国开发计划署于1990年5月发表的第一份《人类发展报告》中公布的。该报告建议，目前对人生发展的衡量应当以人生的三大要素为重点，衡量人生三大要素的指示指标分别采用出生时的预期寿命、成人识字率和实际人均GDP。将以上三个指示指标的数值合成为一个复合指数，即为人文发展指数。
今从1995年世界各国人文发展指数的排列中，选取高发展水平、中等发展水平的国家各五个作为两组样品，另外选四个国家作为待判样品作距离判别分析。
高速发展水平国家：
出生时的预期寿 成人识字率 调整后的人均GDP
美国 76 99 5374
日本 79.5 99 5359
瑞士 78 99 5372
阿根廷 72.1 95.9 5242
阿联酋 73.8 77.7 5370
中等发展中国家：
出生时的预期寿命 成人识字率 调整后的人均GDP
保加利亚 71.2 93 4250
古巴 75.3 94.9 3412
巴拉圭 70 91.2 3390
格鲁吉亚 72.8 99 2300
南非 62.9 80.6 3799
待判国家：
出生时的预期寿命 成人识字率 调整后的人均GDP
中国 68.5 79.3 1950
罗马尼亚 69.9 96.9 2840
希腊 77.6 93.8 5233
哥伦比亚 69.3 90.3 5158
下面有两种实现算法：
1，距离法：
data1<-read.table("clipboard",header=T)
data2<-read.table("clipboard",header=T)
data<-read.table("clipboard",header=T)
f1<-function(x){
mu1<-apply(data1,2,mean)
mu2<-apply(data2,2,mean)
cov1<-cov(data1)
cov2<-cov(data2)
g<-mahalanobis(x,mu2,cov2)-mahalanobis(x,mu1,cov1)
if(g<0) return("中等发展中国家")
else{
if (g>0) return("高速发展水平国家")
else return("无法判断")
}
}
apply(data1,1,f1)
apply(data2,1,f1)
apply(data,1,f1)
2，距离线性判别法：
data1<-read.table("clipboard",header=T)
data2<-read.table("clipboard",header=T)
data<-read.table("clipboard",header=T)
f2<-function(x){
mu_bar<-(apply(data1,2,mean)+apply(data2,2,mean))/2
n1<-dim(data1)[1]
n2<-dim(data2)[2]
cov1<-cov(data1)
cov2<-cov(data2)
cov_matrix<-((n1-1)*cov1+(n2-1)*cov2)/(n1+n2-2)
a<-solve(cov_matrix)%*%(apply(data1,2,mean)-apply(data2,2,mean))
w<-t(x-mu_bar)%*%a
if(w<0) return("中等发展中国家")
else{
if (w>0) return("高速发展水平国家")
else return("无法判断")
}
}
apply(data1,1,f2)
apply(data2,1,f2)
apply(data,1,f2)
