【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，作者：薛毅）。