【R统计】距离判别、Bayes判别和Finsher判别

习题：根据经验，今天和昨天的湿度差X1及今天的压温差X2是预报天气下雨和不下雨的两个重要因素。现有数据如下，今测得x1=8.1，x2=2.0，试问预报明天下雨还是预报明天不下雨？分别用距离判别、Bayes判别（考虑方差相同和方差不同两种情况）和Fisher判别来得到你所需要的结论。

数据：

雨天		非雨天
X1(湿度差)	X2(压温差)	X1(湿度差)	X2((压温差)
-1.9	3.2	0.2	0.2
-6.9	10.4	-0.1	7.5
5.2	2.0	0.4	14.6
5.0	2.5	2.7	8.3
7.3	0.0	2.1	0.8
6.8	12.7	-4.6	4.3
0.9	-15.4	-1.7	10.9
-12.5	-2.5	-2.6	13.1
1.5	1.3	2.6	12.8
3.8	6.8	-2.8	10.0

 

解答：

数据文件 Rain.txt:

-1.9    3.2
-6.9    10.4
5.2    2.0
5.0    2.5
7.3    0.0
6.8    12.7
0.9    -15.4
-12.5    -2.5
1.5    1.3
3.8    6.8

数据文件  Fine.txt

0.2    0.2
-0.1    7.5
0.4    14.6
2.7    8.3
2.1    0.8
-4.6    4.3
-1.7    10.9
-2.6    13.1
2.6    12.8
-2.8    10.0

函数文件 discriminiant.distance.R  

discriminiant.distance<-function
   (TrnX1, TrnX2, TstX = NULL, var.equal = FALSE){
   if (is.null(TstX) == TRUE) TstX<-rbind(TrnX1,TrnX2)
   if (is.vector(TstX) == TRUE)  TstX<-t(as.matrix(TstX))
   else if (is.matrix(TstX) != TRUE)
      TstX<-as.matrix(TstX)
   if (is.matrix(TrnX1) != TRUE) TrnX1<-as.matrix(TrnX1)
   if (is.matrix(TrnX2) != TRUE) TrnX2<-as.matrix(TrnX2)

   nx<-nrow(TstX)
   blong<-matrix(rep(0, nx), nrow=1, byrow=TRUE, 
         dimnames=list("blong", 1:nx))
   mu1<-colMeans(TrnX1); mu2<-colMeans(TrnX2) 
   if (var.equal == TRUE  || var.equal == T){
      S<-var(rbind(TrnX1,TrnX2))
      w<-mahalanobis(TstX, mu2, S)-mahalanobis(TstX, mu1, S)
   }
   else{
      S1<-var(TrnX1); S2<-var(TrnX2)
      w<-mahalanobis(TstX, mu2, S2)-mahalanobis(TstX, mu1, S1)
   }
   for (i in 1:nx){
      if (w[i]>0)
          blong[i]<-1
      else
          blong[i]<-2
   }
   blong
}

　

函数文件 discriminiant.bayes.R

discriminiant.bayes<-function
   (TrnX1, TrnX2, rate=1, TstX = NULL, var.equal = FALSE){
   if (is.null(TstX) == TRUE) TstX<-rbind(TrnX1,TrnX2)
   if (is.vector(TstX) == TRUE)  TstX<-t(as.matrix(TstX))
   else if (is.matrix(TstX) != TRUE)
      TstX<-as.matrix(TstX)
   if (is.matrix(TrnX1) != TRUE) TrnX1<-as.matrix(TrnX1)
   if (is.matrix(TrnX2) != TRUE) TrnX2<-as.matrix(TrnX2)

   nx<-nrow(TstX)
   blong<-matrix(rep(0, nx), nrow=1, byrow=TRUE, 
         dimnames=list("blong", 1:nx))
   mu1<-colMeans(TrnX1); mu2<-colMeans(TrnX2) 
   if (var.equal == TRUE  || var.equal == T){
      S<-var(rbind(TrnX1,TrnX2)); beta<-2*log(rate)
      w<-mahalanobis(TstX, mu2, S)-mahalanobis(TstX, mu1, S)
   }
   else{
      S1<-var(TrnX1); S2<-var(TrnX2)
      beta<-2*log(rate)+log(det(S1)/det(S2))
      w<-mahalanobis(TstX, mu2, S2)-mahalanobis(TstX, mu1, S1)
   }
   
   for (i in 1:nx){
      if (w[i]>beta)
          blong[i]<-1
      else
          blong[i]<-2
   }
   blong
}

　　

函数文件 discriminiant.fisher.R

discriminiant.fisher<-function(TrnX1, TrnX2, TstX = NULL){
   if (is.null(TstX) == TRUE)    TstX<-rbind(TrnX1,TrnX2)
   if (is.vector(TstX) == TRUE)  TstX<-t(as.matrix(TstX))
   else if (is.matrix(TstX) != TRUE)
      TstX<-as.matrix(TstX)
   if (is.matrix(TrnX1) != TRUE)  TrnX1<-as.matrix(TrnX1)
   if (is.matrix(TrnX2) != TRUE)  TrnX2<-as.matrix(TrnX2)

   nx<-nrow(TstX)
   blong<-matrix(rep(0, nx), nrow=1, byrow=TRUE, 
         dimnames=list("blong", 1:nx))
   n1<-nrow(TrnX1); n2<-nrow(TrnX2) 
   mu1<-colMeans(TrnX1); mu2<-colMeans(TrnX2) 
   S<-(n1-1)*var(TrnX1)+(n2-1)*var(TrnX2) 
   mu<-n1/(n1+n2)*mu1+n2/(n1+n2)*mu2
   w<-(TstX-rep(1,nx) %o% mu) %*% solve(S, mu2-mu1);
   for (i in 1:nx){
       if (w[i]<=0)
           blong[i]<-1
       else
           blong[i]<-2
   }
   blong
}

　　

运行脚本：

##原始数据
rain <- read.table("Rain.txt");
fine <- read.table("Fine.txt");
test <- c(8.1, 2.0);

#距离判别
source("discriminiant.distance.R")
discriminiant.distance(rain, fine, test, var.equal=TRUE);
# rain
discriminiant.distance(rain, fine, test);
#rain

#Bayes判别
source("discriminiant.bayes.R");
discriminiant.bayes(rain, fine, 1, test, var.equal=TRUE);
#rain
discriminiant.bayes(rain, fine, 1, test);
#rain

#Fish 判别
source("discriminiant.fisher.R");
discriminiant.fisher(rain, fine, test);
#rain

#综上，预报明天下雨

　　

输出均为 第一种类型（1），也就是下雨（rain）。

　

博文源代码和习题均来自于教材《统计建模与R软件》（ISBN：9787302143666，作者：薛毅）。