R 使用决策树、距离判断、bayes贝叶斯、Fisher 进行分类预测

R对于分类可以使用 距离判断，bayes贝叶斯，fisher，决策树等方法， 下面分别来演示

使用决策树方法对一份儿童纠正脊柱手术数据进行分析

install.packages("rpart")

library(rpart)

fit <- rpart(Kyphosis ~ Age + Number + Start, data = kyphosis)

plot(fit)

text(fit, use.n = TRUE)

距离判断，bayes贝叶斯，fisher 通过下面2个题目来说明，答案中使用的函数在文章后面有给出

rain.data <- data.frame(X1=c(-1.9,-6.9,5.2,5.0,7.3,6.8,0.9,-12.5,1.5,3.8),X2=c(3.2,10.4,2.5,2.0,0.0,12.7,-15.4,-2.5,1.3,6.8))

nonrain.data <- data.frame(X1=c(0.2,-0.1,0.4,2.7,2.1,-4.6,-1.7,-2.6,2.6,-2.8),X2=c(0.2,7.5,14.6,8.3,0.8,4.3,10.9,13.1,12.8,10.0))

distancePredit <- discriminiant.distance(rain.data,nonrain.data,data.frame(X1=8.1,X2=2.0))

distancePredit

#距离判断明天下雨

bayesPredit <- discriminiant.bayes(rain.data,nonrain.data,TstX=data.frame(X1=8.1,X2=2.0),var.equal=TRUE)

bayesPredit

#bayes 方差相同时判断为明天下雨

bayesPredit1 <- discriminiant.bayes(rain.data,nonrain.data,TstX=data.frame(X1=8.1,X2=2.0))

bayesPredit1

#bayes 方差不相同时判断为明天下雨

fisherPredit <- discriminiant.fisher(rain.data,nonrain.data,TstX=data.frame(X1=8.1,X2=2.0))

fisherPredit

#fisher 判断为明天下雨

data <- data.frame(type=c(1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,3,3,3,3,3),

                   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.80,8.68,5.67,3.71,5.37,9.89,5.22,4.71,4.71,3.36,8.27),

                   x2=c(261.01,185.39,249.58,137.13,231.34,231.38,260.25,259.91,273.84,303.59,231.03,308.9,258.69,355.54,316.32,274.57,409.42,330.34,331.47,352.5,347.31,189.56),

                   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.03,17.12,16.75,19.47,18.19,21.26,20.19,17.9,12.74),

                   x4=c(7.36,5.99,6.11,4.35,8.79,8.53,10.02,9.79,8.79,8.53,6.15,8.49,7.16,9.43,8.17,9.67,10.49,9.61,13.72,11,11.19,6.4))

distancePredict1 <- distinguish.distance(data[,-1],as.factor(data[,1]),TstX = NULL, var.equal = TRUE)

distancePredict1

#使用距离判断（方差相同时），有5个误判，误判率23%

distancePredict2 <- distinguish.distance(data[,-1],as.factor(data[,1]),TstX = NULL, var.equal = FALSE)

distancePredict2

#使用距离判断（方差相同时），有3个误判，误判率24%

bayesPredict1 <- distinguish.bayes(data[,-1],as.factor(data[,1]),p=c(11/23,7/23,5/23),TstX = NULL, var.equal = TRUE)

bayesPredict1

#使用bayes判断（方差相同时），有5个误判，误判率23%

bayesPredict2 <- distinguish.bayes(data[,-1],as.factor(data[,1]),p=c(11/23,7/23,5/23),TstX = NULL, var.equal = FALSE)

bayesPredict2

#使用bayes判断（方差相同时），有7个误判，误判率30%

下面为使用的函数

在程序中，输入变量TrnXl、TrnX2表示X1类、X2类训练样本，其输入格式

足数据框，或矩阵（样本按行输入）rate=(L(1/2)/L(2/1)).(p2/p1) 缺省值为1 TstX足待测样

本，其输入格式是数据框，或矩阵（样本接行输入），或向量（一个待测样本）如

果不输入TstX（缺省值）．则待测样本为两个训练样本之和，即计算训练样本的回

代情况输入变量var.equal是逻辑变量，var.equal:TRUE表示认为两总体的

协方差阵足相同的；否则（缺省值）是不同的函数的输出足由“l“和；2”构成

的的一维矩阵，“1”表示待测样本属于x,类，“2”表示待测样本属于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

}

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.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

}

参数说明：

程序分别考虑了总体协方差阵相同和协方差阵不同的情况输入变量TrnX表示

训练样本，其输入格式是矩阵（样本按行输入），或数据框  TrnG是因子变量，表

示训练样本的分类情况输入变量，p是先验概率，缺省值均为1输入变量TstX

是待测样本，其输入格式是矩阵（样本按行输入），或数据框，或向量（一个待测样

本）如果不输入TstX(缺省值)、则待测样本为训练样本输入变量var.equal是

逻辑变量，var.equal=TRUE表示认为总体协方差阵是相同的；否则（缺省值）是

不同的函数的输出是由数字构成的的一维矩阵，数字表示相应的类为了与前

面两总体的判别程序兼容，对于二分类问题，也可以按照discriminiant.bayes

函数的输入格式输入

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

}

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

}