R語言判別分析

本文中分三個(gè)方法介紹判別分析，Bayes判別，距離判別，F(xiàn)isher判別。前兩種判別方法都要考慮兩個(gè)、或多個(gè)總體協(xié)方差（這里是算方差，方差是協(xié)方差的一種）相等或不等的情況，由var.equal=的邏輯參數(shù)表示，默認(rèn)是FALSE，表示認(rèn)為兩總體協(xié)方差不等。用樣本的協(xié)方差可以估計(jì)總體的協(xié)方差。在Bayes方法中我們把相等和不等的兩個(gè)結(jié)果都列了出來，距離判別里我們默認(rèn)兩總體協(xié)方差不等。事實(shí)上，一般使用時(shí)，我們都以兩總體的協(xié)方差不等作為標(biāo)準(zhǔn)來進(jìn)行后續(xù)計(jì)算。

協(xié)方差計(jì)算公式： S(xy) = N,i ∑ [(xi-x均)(yi-y均)]/(N-1) ; 這里要計(jì)算的總體協(xié)方差就是方差，為 S = N,i ∑ [(xi-x均)*(xi-x均)^T]/(N-1)
1. Bayes判別

Bayes判別是假定對(duì)研究對(duì)象已有一定的認(rèn)識(shí)，這種認(rèn)識(shí)常用先驗(yàn)概率來描述，當(dāng)取得樣本后，就可以用樣本來修正已有的先驗(yàn)概率分布，得出后驗(yàn)概率分布，現(xiàn)通過后驗(yàn)概率分布進(jìn)行各種統(tǒng)計(jì)推斷.

* 1.1 兩總體判別的Bayes判別程序*

在程序中，輸入變量TrnX1，TrnX2表示X1類，X2類訓(xùn)練樣本，其輸入格式是數(shù)據(jù)框，或矩陣（樣本按行輸入）.rate= L(1|2) / L(2|1)*p2/p1 ， 缺省值為1. TstX是待測(cè)樣本，其輸入格式是數(shù)據(jù)框，或矩陣（樣本按行輸入），或向量（一個(gè)待測(cè)樣本）.如果不輸入TstX（缺省值），則待測(cè)樣本為兩個(gè)訓(xùn)練樣本之和，即計(jì)算訓(xùn)練樣本的回代情況.輸入變量var.equal是邏輯變量，var.equal=TRUE表示認(rèn)為兩總體的協(xié)方差陣是相同的；否則（缺省值）是不同的. 函數(shù)的輸出是由”1”,”2”構(gòu)成的一維矩陣，”1”代表待測(cè)樣本屬于X1類，”2”代表待測(cè)樣本屬于X2類

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
}

* 1.2 例. 用1.1中Bayes判別程序?qū)谛牟『驼Ｈ诉M(jìn)行判別分類*

為研究舒張期血壓與血漿膽固醇對(duì)冠心病的作用，調(diào)查了50~59歲的女冠心病人15例和正常人16例。他們的舒張壓（x1）與血漿膽固醇（x2）列在下表中。試用判別分析法建立判別冠心病與正常人的判別函數(shù)

#它們的先驗(yàn)概率pi分別用15/31、16/31來估計(jì)
#載入數(shù)據(jù)
TrnX1<-matrix(
   c(9.86, 13.33, 14.66, 9.33, 12.80, 10.66, 10.66, 13.33, 13.33, 13.33, 12.00, 14.66, 13.33, 12.80, 13.33,
     5.18, 3.73, 3.89, 7.10, 5.49, 4.09, 4.45, 3.63, 5.96, 5.70, 6.19, 4.01, 4.01, 3.63, 5.96),
   ncol=2)
TrnX2<-matrix(
   c(10.66, 12.53, 13.33, 9.33, 10.66, 10.66, 9.33, 10.66, 10.66, 10.66, 10.40, 9.33, 10.66, 10.66, 11.20, 9.33,
     2.07, 4.45, 3.06, 3.94, 4.45, 4.92, 3.68, 2.77, 3.21, 5.02, 3.94, 4.92, 2.69, 2.43, 3.42, 3.63),
   ncol=2)
#輸入待測(cè)樣本TrnX
TstX<-matrix(
   c(9.06,13.00,12.66,9.00,12.12,11.66,12.11,12.63,8.33,11.12,
     5.68,3.43,2.82,6.86,5.19,2.17,4.12,3.26,2.91,4.22 ),
   ncol=2)
#載入兩總體的貝葉斯判別函數(shù)    注 把貝葉斯判別函數(shù)存在了wd中
source("discriminiant.bayes.R")
#### 總體協(xié)方差陣相同
discriminiant.bayes(TrnX1, TrnX2,rate=16/15,var.equal=TRUE)    #rate=L(1|2)/L(2|1)* p2/p1
discriminiant.bayes(TrnX1, TrnX2,TstX,rate=16/15,var.equal=TRUE)
#### 總體協(xié)方差陣不同
discriminiant.bayes(TrnX1, TrnX2,rate=16/15)
discriminiant.bayes(TrnX1, TrnX2,TstX,rate=16/15)

在判別函數(shù)的鑒定下，有5個(gè)被錯(cuò)判，3個(gè)樣本從冠心病組被錯(cuò)判為正常組，兩個(gè)樣本從正常組錯(cuò)判到冠心病組，分別為1,6,7,17,18

待測(cè)樣本分類結(jié)果為：

                 1 2 3 4 5 6 7 8 9 10  （樣本協(xié)方差陣相同）
           blong 2 1 2 1 1 2 1 2 2  2

                 1 2 3 4 5 6 7 8 9 10  （樣本協(xié)方差陣不同）
           blong 1 1 2 1 1 2 1 2 2  2

* 1.3 用R包里的lda函數(shù)進(jìn)行判別，得到判別函數(shù)*

#也可以csv格式導(dǎo)入數(shù)據(jù)，先將以上兩類數(shù)據(jù)存為csv格式文檔
coronary_disease <- read.csv("coronary_disease.csv")
#把分組變量變?yōu)槎ㄐ宰兞?br /> group <- factor(coronary_disease$group)
#隨機(jī)抽取20個(gè)一般樣本做訓(xùn)練樣本
#train <- sample(1:31,20)
#顯示訓(xùn)練樣本中各類的比例
#table(group[train])
#group作為分組變量，X1,X2作為判別變量，使用訓(xùn)練樣本生成判別函數(shù)，先驗(yàn)概率各為50%
library(MASS)
Z <- lda(group~.,data = coronary_disease ,prior=c(1,1)/2) #subset=train
Z
判別函數(shù)是Z = -0.6379195 * X1 + -0.8001452 * X2

* 1.4 取原數(shù)據(jù)80%做train，20%做test，查看結(jié)果*

coronary_disease <- read.csv("coronary_disease.csv")
#讀取行數(shù)
N = length(coronary_disease$group)      
#ind=1的是0.7概率出現(xiàn)的行，ind=2是0.3概率出現(xiàn)的行
ind=sample(2,N,replace=TRUE,prob=c(0.8,0.2))
#生成訓(xùn)練集(這里訓(xùn)練集和測(cè)試集隨機(jī)設(shè)置為原數(shù)據(jù)集的80%,20%)
coronary_train <- coronary_disease[ind==1,]
#生成測(cè)試集
coronary_test <- coronary_disease[ind==2,]
#固定這個(gè)28分組
set.seed(7)
#可以看到這里的train中，i個(gè)1組，j個(gè)2組，共i+j個(gè)數(shù)據(jù)，則test里有(31-i-j個(gè))數(shù)據(jù)
i <- table(coronary_train$group)[[1]]
j <- table(coronary_train$group)[[2]]
i;j

* 1.4.1 用1.1中Bayes判別程序?qū)y(cè)試集進(jìn)行判別分類*

#輸入train的組別1、2，TrnXi；輸入test
TrnX1 <- coronary_train[1:i,2:3]
TrnX2 <- coronary_train[(i+1):(i+j),2:3]
TstX <- coronary_test[,2:3]
source("discriminiant.bayes.R")
#### 總體協(xié)方差陣相同
##它們的先驗(yàn)概率pi分別用i/24、j/24來估計(jì)
discriminant_train1 <- discriminiant.bayes(TrnX1, TrnX2,rate=j/i,var.equal=TRUE) ; discriminant_train1  #rate=L(1|2)/L(2|1)* p2/p1
discriminant_test1 <- discriminiant.bayes(TrnX1, TrnX2,TstX,rate=j/i,var.equal=TRUE) ; discriminant_test1
#### 總體協(xié)方差陣不同
discriminant_train2 <- discriminiant.bayes(TrnX1, TrnX2,rate=j/i) ; discriminant_train2
discriminant_test2 <- discriminiant.bayes(TrnX1, TrnX2,TstX,rate=j/i) ; discriminant_test2
* 1.4.2 導(dǎo)出結(jié)果到csv文件中*

#在原始訓(xùn)練集，測(cè)試集中加入一列Is，表示是訓(xùn)練集還是測(cè)試集
coronary_train$Is <- c(rep("train",(i+j))) ;coronary_test$Is <- c(rep("test",(31-i-j)))
#在原始訓(xùn)練集，測(cè)試集中分別加入總體協(xié)方差相同或不同時(shí)候的判定組別Discriminant
coronary_train$Belong_VarSame <- discriminant_train1[1:(i+j)]
coronary_test$Belong_VarSame <- discriminant_test1[1:(31-i-j)]
coronary_train$Belong_VarDiff <- discriminant_train2[1:(i+j)]
coronary_test$Belong_VarDiff <- discriminant_test2[1:(31-i-j)]
#合成一張表
coronary_result <- rbind(coronary_train,coronary_test)
write.csv(coronary_result,"coronary_result.csv")
* 1.4.3 總體協(xié)方差相同或不同時(shí)，train和test的正確率*

#以train里的Var相同情況為例：
true_value_VarSame_train = coronary_result[coronary_result$Is=="train",]$Belong_VarSame
predict_value_VarSame_trian = coronary_result[coronary_result$Is=="train",]$group
#計(jì)算模型精確度
error_VarSame_train = predict_value_VarSame_trian-true_value_VarSame_train
accuracy_VarSame_train = (nrow(coronary_train)-sum(abs(error_VarSame_train)))/nrow(coronary_train) #精確度--判斷正確的數(shù)量占總數(shù)的比例
accuracy_VarSame_train

2. 距離判別（二分類與多分類）

* 2.1 兩分類的距離判別程序*

在程序中，輸入變量TrnX1，TrnX2表示X1類，X2類訓(xùn)練樣本，其輸入格式是數(shù)據(jù)框，或矩陣（樣本按行輸入）.TstX是待測(cè)樣本，其輸入格式是數(shù)據(jù)框，或矩陣（樣本按行輸入），或向量（一個(gè)待測(cè)樣本）.如果不輸入TstX（缺省值），則待測(cè)樣本為兩個(gè)訓(xùn)練樣本之和，即計(jì)算訓(xùn)練樣本的回代情況.輸入變量var.equal是邏輯變量，var.equal=TRUE表示認(rèn)為兩總體的協(xié)方差陣是相同的；否則（缺省值）是不同的. 函數(shù)的輸出是由”1”,”2”構(gòu)成的一維矩陣，”1”代表待測(cè)樣本屬于X1類，”2”代表待測(cè)樣本屬于X2類

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
}

* 2.2 多分類的距離判別程序*

distinguish.distance<-function
   (TrnX, TrnG, TstX = NULL, var.equal = FALSE){
   if ( is.factor(TrnG) == FALSE){
       mx<-nrow(TrnX); mg<-nrow(TrnG)
       TrnX<-rbind(TrnX, TrnG)
       TrnG<-factor(rep(1:2, c(mx, mg)))
   }
   if (is.null(TstX) == TRUE) TstX<-TrnX
   if (is.vector(TstX) == TRUE)  TstX<-t(as.matrix(TstX))
   else if (is.matrix(TstX) != TRUE)
      TstX<-as.matrix(TstX)
   if (is.matrix(TrnX) != TRUE) TrnX<-as.matrix(TrnX)

   nx<-nrow(TstX)
   blong<-matrix(rep(0, nx), nrow=1, dimnames=list("blong", 1:nx))
   g<-length(levels(TrnG))
   mu<-matrix(0, nrow=g, ncol=ncol(TrnX))
   for (i in 1:g)
      mu[i,]<-colMeans(TrnX[TrnG==i,])
   D<-matrix(0, nrow=g, ncol=nx)
   if (var.equal == TRUE  || var.equal == T){
      for (i in 1:g)
         D[i,]<- mahalanobis(TstX, mu[i,], var(TrnX),tol=2e-21)
   }
   else{
      for (i in 1:g)
         D[i,]<- mahalanobis(TstX, mu[i,], var(TrnX[TrnG==i,]),tol=2e-21)
   }
   for (j in 1:nx){
      dmin<-Inf
      for (i in 1:g)
          if (D[i,j]<dmin){
             dmin<-D[i,j]; blong[j]<-i
      }
   }
   blong
}

* 2.3 例2. 多分類的距離判別實(shí)例，判別冠心病和正常人*

例. 從心電圖的5個(gè)不同指標(biāo)（x1~x5）對(duì)健康人（C=1）、血管硬化癥患者（C=2）和冠心病患者（C=3）的數(shù)據(jù)進(jìn)行判別分析，使用典型判別法。

#用和張佳玉sas文檔里一樣的數(shù)據(jù)輸入
X <- data.frame(
   x1 = c(8.11,9.36,9.85,2.55,6.01,9.64,4.11,8.9,7.71,7.51,8.06,6.8,8.68,5.67,8.1,3.71,5.37,5.22,4.71,4.71,3.36,8.27),
   x2 = c(251.01,185.39,249.58,137.13,231.34,231.38,260.25,259.51,273.84,303.59,231.03,308.9,258.69,255.54,476.69,316.12,274.57,330.34,331.47,352.5,347.31,189.59),
   x3 = c(13.23,9.02,15.61,9.21,14.27,13.03,14.72,14.16,16.01,19.14,14.41,15.11,14.02,15.13,7.38,17.12,16.75,18.19,21.26,20.79,17.9,12.74),
   x4 = c(5.46,5.66,6.06,6.11,5.21,4.86,5.36,4.91,5.15,5.7,5.72,5.52,4.79,4.97,5.32,6.04,4.98,4.96,4.3,5.07,4.65,5.46),
   x5 = c(7.31,5.99,6.11,4.35,8.79,8.53,10.03,9.79,8.79,8.53,6.15,8.49,7.16,9.43,11.32,8.17,9.67,9.61,13.72,11,11.19,6.94)
)
G <- factor(c("1","1","1","1","1","1","1","1","1","1","1","2","2","2","2","2","2","3","3","3","3","3"),levels = c("1","2","3"))
source("distinguish.distance.R")
distinguish.distance(X,G)   #默認(rèn)協(xié)方差不同

#用薛毅書中例子里的數(shù)據(jù)輸入
X1 <- iris[,1:4]
G1 <- gl(3,50)
source("distinguish.distance.R")
distinguish.distance(X1,G1)    #默認(rèn)協(xié)方差不同

這里用薛毅的數(shù)據(jù)時(shí)沒有問題，換成自己的數(shù)據(jù)會(huì)一直報(bào)錯(cuò)，后來找到了原因。

distinguish.distance.R程序中有兩行代碼是mahalanobis(TstX, mu[i,], var(TrnX[TrnG==i,]))，因?yàn)槲覀兊臄?shù)據(jù)太小，出現(xiàn)如下錯(cuò)誤：
（Error in solve.default(cov, …) :
system is computationally singular: reciprocal condition number = 2.60317e-20）
mahalanobis()函數(shù)會(huì)將這個(gè)很小的值認(rèn)為是0，無法繼續(xù)后面的運(yùn)算，所以報(bào)錯(cuò)

解決辦法是，加一個(gè)參數(shù),tol=2e-21，這個(gè)值設(shè)置成比error里的值小就行。
3. Fisher判別

在程序中，輸入變量TrnX1，TrnX2表示X1類，X2類訓(xùn)練樣本，其輸入格式是數(shù)據(jù)框，或矩陣（樣本按行輸入）.TstX是待測(cè)樣本，其輸入格式是數(shù)據(jù)框，或矩陣（樣本按行輸入），或向量（一個(gè)待測(cè)樣本）.如果不輸入TstX（缺省值），則待測(cè)樣本為兩個(gè)訓(xùn)練樣本之和，即計(jì)算訓(xùn)練樣本的回代情況. 函數(shù)的輸出是由”1”,”2”構(gòu)成的一維矩陣，”1”代表待測(cè)樣本屬于X1類，”2”代表待測(cè)樣本屬于X2類

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
}

* 3.1 例3. Fisher判別分析法建立心肌梗塞病和正常人的判別函數(shù)*

為研究心肌梗塞病的危險(xiǎn)因素，某研究者考察了兩組人群（心梗組與正常組，兩組間受試對(duì)象的年齡與性別構(gòu)成接近）有關(guān)的指標(biāo)10多項(xiàng)，現(xiàn)取其中血脂方面的6項(xiàng)指標(biāo)：tc（總膽固醇）、tg（甘油三酯）、hdlc（高密度脂蛋白膽固醇）、ldlc（低密度脂蛋白膽固醇）、apoa（載脂蛋白AI）、apob（載脂蛋白B）。指標(biāo)的測(cè)定結(jié)果如下表所示，每組各取30例，試作判別分析。

classX1 <- data.frame(
    x1 = c(245,236,238,233,240,235,204,200,297,177,200,195,166,144,233,143,228,264,178,240,180,161,236,168,174,215,268,213,285,193
),
    x2 = c(157,275,354,250,149,166,365,95,240,97,172,211,217,111,107,91,223,186,131,127,211,91,106,106,141,168,185,387,154,123),
    x3 = c(38,40,38,31,35,40,38,43,38,49,43,47,33,28,42,24,34,41,49,33,27,39,36,36,28,38,28,22,39,42),
    x4 = c(168,125,126,150,170,164,90,100,207,108,116,106,86,46,156,108,136,183,98,174,106,88,104,104,103,134,203,141,210,121),
    x5 = c(1.1,1.22,0.9,1.02,1.26,1.3,1.33,1.24,1.14,1.49,1.25,1.22,1.1,0.71,0.95,0.67,1.05,1.22,1.18,0.78,0.85,0.94,0.87,0.87,0.81,0.88,0.75,0.8,1.17,1.12),
    x6 = c(1.01,1.12,1.06,0.98,1.13,1.15,0.95,0.98,1.51,1.02,1.03,0.94,0.74,0.65,0.77,0.65,0.84,0.92,1.27,0.9,0.69,0.52,0.58,0.73,0.73,0.87,0.97,0.78,1.37,1)
)
classX2 <- data.frame(
    x1 = c(174,106,173,178,198,180,134,204,168,219,189,180,177,172,166,210,166,223,136,156,201,134,195,262,194,165,183,200,171,222),
    x2 = c(140,110,82,100,112,114,60,118,80,157,158,90,227,55,217,166,217,186,72,107,117,58,93,257,171,70,249,191,309,350),
    x3 = c(47,52,53,43,53,48,36,63,52,28,43,59,75,51,33,42,33,73,67,45,45,60,51,62,42,36,44,58,52,13),
    x4 = c(120,40,103,117,123,110,84,119,90,142,115,102,64,102,86,130,86,113,46,106,147,65,141,142,114,110,88,100,51,57),
    x5 = c(0.84,1.08,0.97,0.98,0.98,1.02,0.98,1.02,1.07,1.02,0.92,1.32,1.4,1.31,1.1,1.28,1.1,1.62,1.45,0.93,1.06,1.03,1.22,1.56,1.11,1.22,1.12,1.61,1.37,0.36),
    x6 = c(0.57,0.87,0.66,0.65,0.72,0.8,0.58,0.84,0.8,0.83,0.8,0.9,0.99,0.97,0.74,1.02,0.74,0.98,0.84,0.74,0.85,0.54,0.72,0.8,0.71,0.96,0.96,0.77,0.69,1.39)
)
source("discriminiant.fisher.R")
discriminiant.fisher(classX1,classX2)

將訓(xùn)練樣本回代進(jìn)行判別，有10個(gè)點(diǎn)判錯(cuò)，分別是10、19、22、24、40、45、46、47、56、60號(hào)樣本.

#粗略導(dǎo)出以上三個(gè)例子里的結(jié)果
write.csv(discriminiant.bayes(TrnX1, TrnX2,TstX,rate=16/15,var.equal=TRUE),"disease1.csv")
write.csv(distinguish.distance(X1,G1),"disease2.csv")
write.csv(discriminiant.fisher(classX1,classX2),"disease3.csv")

R里還可以用knn來進(jìn)行判別，百度有很多博文有詳細(xì)介紹，這里不作討論。