基于R语言的判别分析

本文主要介绍了基于R语言实现距离判别、Bayes判别、Fisher判别的基本思路以及给出了具体的操作过程。

1.数据

这里总共有个20个电视品牌的数据，销售状态G1中的1表示畅销，2表示滞销；销售状态G2中的1表示畅销，2表示平销，3表示滞销。接下来只对销售状态G1做判别分析，也就是两总体的判别分析，对于多总体判别分析，思路完全一样，有兴趣的同学可以自己尝试。

2.描述性统计分析

code:

##   step1：对样本做描述性统计分析


par(mfrow = c(1,3))

##  Q和C的散点图
plot(data1[,1],data1[,2],xlab = "质量评分Q",ylab = "功能评分C",main = "Q和C的散点图");
text(data1[,1],data1[,2],data1[,4],adj = -0.8)

##  Q和P的散点图
plot(data1[,1],data1[,3],xlab = "质量评分Q",ylab = "销售价格P",main = "Q和P的散点图");
text(data1[,1],data1[,3],data1[,4],adj = -0.8)

##  C和P的散点图

plot(data1[,2],data1[,3],xlab = "功能评分C",ylab = "销售价格P",main = "C和P的散点图");
text(data1[,2],data1[,3],data1[,4],adj = -0.8)

运行结果：

这里不对描述性结果做解释.

接下来进入正题，依次为距离判别，Bayes判别和Fisher判别，必要细节在程序中都给出了相应注释，这里只对结果做简单解释，对程序不做过多的解释。

##   step2： 选择方法：  =====  距离判别法  =====



#  First：编写判别程序   【注】：该函数只适用于数据框或者矩阵

discriminant_distance <- function(trnx1 , trnx2 , tst = NULL , var.equal = FALSE){
  nx <- nrow(tst)
  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){
    n1 <- nrow(trnx1)
    n2 <- nrow(trnx2)
    S <- (n1 + n2 -1)/(n1 + n2 -2)*var(rbind(trnx1 , trnx2))
    w <- mahalanobis(tst , mu2 , S) - mahalanobis(tst , mu1 , S)
  }else{
    S1 <- var(trnx1)
    S2 <- var(trnx2)
    w <- mahalanobis(tst , mu2 , S2) - mahalanobis(tst , mu1 , S1)
  }
  for(i in 1:nx){
    if(w[i]>0) blong[i] <- 1
    else blong[i] <- 2
  }
  blong
}



#  Second：构造训练样本和测试样本

##  1. 取状态为1的前3个样本，状态为2的后两个样本作为测试样本

trnx1 <- data1[1:13,1:3]
trnx2 <- data1[14:20,1:3]

tst <- rbind(data1[1:3,1:3],data1[19:20,1:3])


#（1）  在协方差矩阵相同的情况下做判别

discriminant_distance(trnx1,trnx2,tst,var.equal = T)

#（2）  在协方差矩阵不同的情况下做判别

discriminant_distance(trnx1,trnx2,tst)


#  2. 用训练集作为测试集

tst1 <- data1[1:20,1:3]

#（1）  在协方差矩阵相同的情况下做判别

discriminant_distance(trnx1,trnx2,tst1,var.equal = T)

#（2）  在协方差矩阵不同的情况下做判别


discriminant_distance(trnx1,trnx2,tst1,var.equal = FALSE)




##   step3： 选择方法：  =====  Bayes判别法  =====


data1 <- read.csv("E:/多元统计分析/tvdata.csv")
names(data1) <- c("Q","C","P","G1","G2")


##  编写判别程序并用样本进行判别   【注】：该程序只适用于数据框或者矩阵

#（1）  在协方差矩阵相同的情况下做判别

a <- data1[1:13,1:3]
b <- data1[14:20,1:3]

n1 <- nrow(a)
n2 <- nrow(b)

mu1 <- colMeans(a)
mu2 <- colMeans(b)

sig1 <- var(a)
sig2 <- var(b)
sig <- ((n1-1)*sig1+(n2-1)*sig2)/(n1+n2-2)

beta <- log2(n2/n1)

ab <- rbind(a,b)
nx <- nrow(ab)

blong <- matrix(rep(0,nx),nrow = 1,byrow = TRUE,dimnames = list('blong',1:nx))

for(i in 1:nx){
  x <- ab[i,]
  wx <- t(t(x)-(1/2)*(mu1+mu2)) %*% solve(sig) %*% (mu1-mu2)
  if(wx > beta){
    blong[i] <- 1
  }else{
    blong[i] <- 2
  }
}

blong




#（2）  在协方差矩阵不相同的情况下做判别

p1 <- n1/(n1+n2)
p2 <- n2/(n1+n2)

##  计算k值

k <- log2(p2/p1)+(1/2)*log2(det(sig1)/det(sig2))+
  (1/2)*(t(mu1)%*%solve(sig1)%*%mu1 - t(mu2)%*%solve(sig2)%*%mu2)


##  bayes判别

ab <- rbind(a,b)
ab <- as.matrix(ab)
nx <- nrow(ab)

blong <- matrix(rep(0,nx),nrow = 1,byrow = TRUE,dimnames = list('blong',1:nx))

for(i in 1:nx){
  x <- ab[i,]
  
  wx <- (-1/2)*t(x) %*% (solve(sig1)-solve(sig2)) %*% x + 
    (t(mu1)%*%solve(sig1) - t(mu2)%*%solve(sig2))%*%x
  
  if(wx > k){
    blong[i] <- 1
  }else{
    blong[i] <- 2
  }
}

blong



##   step4： 选择方法：  =====  Fisher判别法  =====



data1 <- read.csv("E:/多元统计分析/tvdata.csv")
names(data1) <- c("Q","C","P","G1","G2")


library(MASS)


##   (1)  协方差矩阵相等时为一次判别函数


model1 <- lda(data1$G1~data1$Q+data1$C+data1$P)
z <- predict(model1)
newG1 <- z$class
cbind(data1$G1,z$x,newG1)



##   (2)  协方差矩阵不相等时为二次判别函数

model2 <- qda(data1$G1~data1$Q+data1$C+data1$P)
Z <- predict(model2)
new_G1 <- Z$class 
cbind(data1$G1,new_G1)

ok，在距离判别中，我们使用了马氏距离，这里的马氏不是马尔科夫，也不是闵可夫斯基，而是马哈拉诺比斯，具体计算公式在任何一本多元统计教材都可以找到。

以上就是基于R语言的判别分析的三种主要判别方法。