【R统计】多类别的判别

8.2 某医院研究心电图指标对健康人（1），硬化病患者（2）和冠心病患者（3）的鉴别能力，先训练样本如下。试用距离判别（考虑方差相同与方差不同两种情况）、Bayes判别（考虑方差同和不同两种情况，且先验概率为11/23,7/23,5/23）对数据进行分析。

原始数据 data.txt如下， 3,4,5,6列为心电图指标，2列为类别

1	1	8.11	261.01	13.23	7.36
2	1	9.36	185.39	9.02	5.99
3	1	9.85	249.58	15.61	6.11
4	1	2.55	137.13	9.21	4.35
5	1	6.01	231.34	14.27	8.79
6	1	9.64	231.38	13.03	8.53
7	1	4.11	260.25	14.72	10.02
8	1	8.90	259.91	14.16	9.79
9	1	7.71	273.84	16.01	8.79
10	1	7.51	303.59	19.14	8.53
11	1	8.06	231.03	14.41	6.15
12	2	6.80	308.90	15.11	8.49
13	2	8.68	258.69	14.02	7.16
14	2	5.67	355.54	15.03	9.43
17	2	3.71	316.32	17.12	8.17
17	2	5.37	274.57	16.75	9.67
18	2	9.89	409.42	19.47	10.49
19	3	5.22	330.34	18.19	9.61
20	3	4.71	331.47	21.26	13.72
21	3	4.71	352.50	20.79	11.00
22	3	3.36	347.31	17.90	11.19
23	3	8.27	189.56	12.74	6.94

函数distinguish.distance.R

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))
   }
   else{
      for (i in 1:g)
         D[i,]<- mahalanobis(TstX, mu[i,], var(TrnX[TrnG==i,]))
   }
   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
}

　　

函数distinguish.bayes.R

distinguish.bayes<-function
   (TrnX, TrnG, p=rep(1, length(levels(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){
         d2 <- mahalanobis(TstX, mu[i,], var(TrnX))
         D[i,] <- d2 - 2*log(p[i])
      }
   }
   else{
      for (i in 1:g){
         S<-var(TrnX[TrnG==i,])
         d2 <- mahalanobis(TstX, mu[i,], S)
         D[i,] <- d2 - 2*log(p[i])-log(det(S))
      }
   }
   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
}

　　

运行脚本

data <- read.table("data.txt");
X <- data[,3:6];
G <- factor(data[,2]);

source("distinguish.distance.R");
distinguish.distance(X,G);
#      1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
#blong 2 1 1 1 1 1 1 1 1  1  1  2  1  2  2  1  2  2  3  3  3  3
distinguish.distance(X,G,var.equal=TRUE);
#      1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
#blong 1 1 1 1 1 1 3 1 1  3  1  2  1  2  2  3  2  2  3  3  3  1

source("distinguish.bayes.R");
distinguish.bayes(X,G, p=c(rep(11/23,11),rep(7/23,7),rep(5/23,5)));
#      1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
#blong 2 1 1 1 1 1 1 1 1  1  1  2  1  2  2  1  2  2  3  3  3  3
distinguish.bayes(X,G, p=c(rep(11/23,11),rep(7/23,7),rep(5/23,5)),var.equal=TRUE);
#      1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
#blong 1 1 1 1 1 1 3 1 1  3  1  2  1  2  2  3  2  2  3  3  3  1

　

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