R语言实战——距离判别、贝叶斯判别、Fisher判别理论详细推导与R语言实现

最新推荐文章于 2022-11-16 18:10:10 发布

三只佩奇不结义

最新推荐文章于 2022-11-16 18:10:10 发布

阅读量2.4w

点赞数 55

分类专栏： R语言随笔机器学习文章标签：数据挖掘机器学习算法判别分析

本文链接：https://blog.csdn.net/qq_41196612/article/details/105339511

版权

R语言随笔同时被 2 个专栏收录

19 篇文章 28 订阅

订阅专栏

机器学习

11 篇文章 2 订阅

订阅专栏

前言

判别分析是用以判别个体所属群体的一种统计方法，它产生于20世纪30年代，近年来，在许多现代自然科学的各个分支和技术部门中，得到广泛的应用.

例如，利用计算机对一个人是否有心脏病进行诊断时，可以取一批没有心脏病的人，测其p个指标的数据，然后再取一批已知患有心胜病的人，同样也测得p个相同指标的数据，利用这些数据建一个判别函数，并求出相应的临界值，这时对于需要进行诊断的人，也同样测其p个指标的数据，将其代入判别函数，求得判别得分，再依判别临界值，即可判断此人是属于有心脏病的那一群体，还是属于没有心脏病的那一群体。又如在考古学中，对化石及文物年代的判断；在地质学中，判断是有矿还是无矿；在量管理中，判断某种产品是合格品，还是不合格品；在植物学中，对于新发现的一种植物，判断其属于那一科.总之判别分析方法在很多学科中有着广泛的应用.

判别方法有多种，这里主要介绍的是最常用的判别方法，而且是两类群体的判别方法.

1 距离判别

在这里插入图片描述

1.1 双群体

1.1.1 理论推导

在这里插入图片描述

1.1.2 R语言实现

discriminiant.distance.R(距离判别函数)

discriminant.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) == FALSE)
    TrnX1 <- as.matrix(TrnX1)
  if (is.matrix(TrnX2) != TRUE)
    TrnX2 <- as.matrix(TrnX2)
  
  # 
  nx <- nrow(TstX)    # 需要用以预测的集合的大小，或者说测试集的大小
  # 生成长度为nx的0向量，用以存储预测的标签
  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(x)
    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
}

在程序中，输入变量Trnx1、Trnx2表示X1类、X2类训练样本，其输入格式是数据框，或矩阵（样本按行输入），输入变量TstX是待测样本，其输入格式是数据框，或矩阵（样本按行输入），或向量（一个待测样本）.如果不输入Tstx（缺省值），则待测样本为两个训练样本之和，即计算训练样本的回代情况。输入变量var.equal是逻辑变量，var.equal=TRUE表示两个总体的协方差阵相同；否则（缺省值）为不同。函数的输出是由“1”和“2”构成的的一维矩阵，“1”表示待测样本属于X1类，“2”表示待测样本属于X2类。

在上述程序中，用到Mahalanobis距离函数mahalanobis（），该函数的使用格式为mahalanobis(x, center, cov, inverted=FLASE)
其中x是由样本数据构成的向量或矩阵（p维），center为样本中心，cov为样本的协方差阵.

1.1.3 实例分析

在这里插入图片描述

R语言代码

classX1 <- data.frame(
  x1=c(6.60, 6.60, 6.10, 6.10, 8.40, 7.2, 8.40, 7.50, 7.50, 8.30, 7.80, 7.80),
  x2=c(39.00,39.00, 47.00, 47.00, 32.00, 6.0, 113.00, 52.00, 52.00,113.00,172.00,172.00),
  x3=c(1.00, 1.00, 1.00, 1.00, 2.00, 1.0, 3.50, 1.00,3.50, 0.00, 1.00, 1.50),
  x4=c(6.00, 6.00, 6.00, 6.00, 7.50, 7.0, 6.00, 6.00,7.50, 7.50, 3.50, 3.00),
  x5=c(6.00, 12.00, 6.00, 12.00, 19.00, 28.0, 18.00, 12.00,6.00, 35.00, 14.00, 15.00),
  x6=c(0.12, 0.12, 0.08, 0.08, 0.35, 0.3, 0.15, 0.16,0.16, 0.12, 0.21, 0.21),
  x7=c(20.00,20.00, 12.00, 12.00, 75.00, 30.0, 75.00, 40.00,40.00,180.00, 45.00, 45.00)
)
classX2<-data.frame(
  x1=c(8.40, 8.40, 8.40, 6.3, 7.00, 7.00, 7.00, 8.30,
       8.30, 7.2, 7.2, 7.2, 5.50, 8.40, 8.40, 7.50,
       7.50, 8.30, 8.30, 8.30, 8.30, 7.80, 7.80),
  x2=c(32.0 ,32.00, 32.00, 11.0, 8.00, 8.00, 8.00,161.00,
       161.0, 6.0, 6.0, 6.0, 6.00,113.00,113.00, 52.00,
       52.00, 97.00, 97.00,89.00,56.00,172.00,283.00),
  x3=c(1.00, 2.00, 2.50, 4.5, 4.50, 6.00, 1.50, 1.50,
       0.50, 3.5, 1.0, 1.0, 2.50, 3.50, 3.50, 1.00,
       1.00, 0.00, 2.50, 0.00, 1.50, 1.00, 1.00),
  x4=c(5.00, 9.00, 4.00, 7.5, 4.50, 7.50, 6.00, 4.00,
       2.50, 4.0, 3.0, 6.0, 3.00, 4.50, 4.50, 6.00,
       7.50, 6.00, 6.00, 6.00, 6.00, 3.50, 4.50),
  x5=c(4.00, 10.00, 10.00, 3.0, 9.00, 4.00, 1.00, 4.00,
       1.00, 12.0, 3.0, 5.0, 7.00, 6.00, 8.00, 6.00,
       8.00, 5.00, 5.00,10.00,13.00, 6.00, 6.00),
  x6=c(0.35, 0.35, 0.35, 0.2, 0.25, 0.25, 0.25, 0.08,
       0.08, 0.30, 0.3, 0.3, 0.18, 0.15, 0.15, 0.16,
       0.16, 0.15, 0.15, 0.16, 0.25, 0.21, 0.18),
  x7=c(75.00,75.00, 75.00, 15.0,30.00, 30.00, 30.00, 70.00,
       70.00, 30.0, 30.0, 30.0,18.00, 75.00, 75.00, 40.00,
       40.00,180.00,180.00,180.00,180.00,45.00,45.00)
)

path <- 'E:\\桌面文档\\学习文件\\大三\\多元统计\\统计建模R语言书籍'
setwd(path)   # 设置工作路径(work directory)
source('distanceDiscriminant.R')
# 方差相同
discriminant.distance(classX1,classX2,var.equal=TRUE)
# 方差不同
discriminant.distance(classX1,classX2)

在认为两总体协方差阵相同的情况下，将训练样本回代进行判别，有三个点判错，分别是第9号样本、第28号样本和第29号样本.

在认为两总体协方差阵不同的情况下，将训练样本回代进行判别，只有一个点判错，是第9号样本。

1.2 多群体

1.2.1 理论推导

在这里插入图片描述

1.2.2 R语言实现

distinguish.distance(多分类问题的距离判别函数)

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)))  # 1重复mx遍，2重复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用于存放预测值
  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,])
  print(mu)
  # 计算马氏距离
  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,]))
  }
  print(D)
  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
}

程序分别考虑了总体协方差阵相同和总体协方差阵不同的两种情况。输入变量TrnX表示训练样本，其输入格式是矩阵（样本按行输入），或数据框.TrnG是因子变量，表示输入训练样本的分类情况.输入变量TstX是待测样本，其输入格式是矩阵（样本按行输入），或数据框，或向量（一个待测样本）.如果不输入TstX（缺省值），则待测样本为训练样本.输入变量var.equal是逻辑变量，var.equal=TRUE表示计算时认为总体协方差阵是相同的；否则（缺省值）是不同的.函数的输出是由数字构成的的一维矩阵，数字表示相应的类.为了与前一个程序兼容，对于二分类问题，也可以按照discriminiant.distance函数的输入格式输入.

1.2.3 实例分析

在这里插入图片描述
R软件中提供了Iris数据，数据的前四列是数据的四个属性，第五列标明数据属于哪一类.

代码：

### 多群体距离判别
X <- iris[,1:4]
G <- gl(3,50)    # 与rep函数有所不同，gl是将数字先重复再接起来
path <- 'E:\\桌面文档\\学习文件\\大三\\多元统计\\统计建模R语言书籍'
setwd(path)
source('distanceDiscriminant2.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 23 24
blong 1 1 1 1 1 1 1 1 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
      25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45
blong  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
      46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
blong  1  1  1  1  1  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
      67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87
blong  2  2  2  2  3  2  3  2  2  2  2  2  2  2  2  2  2  3  2  2  2
      88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105
blong  2  2  2  2  2  2  2  2  2  2  2  2   2   3   3   3   3   3
      106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
blong   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3
      121 122 123 124 125 126 127 128 129 130 131 132 133 134 135
blong   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3
      136 137 138 139 140 141 142 143 144 145 146 147 148 149 150
blong   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3

从计算结果可以看出，只有第71号样本、第73号样本和第84号样本错判，回代的判别正确率为147/150=98%.

2 贝叶斯判别

Bayes判别是假定对研究对象已有一定的认识，这种认识常用先验概率来描述，当取得样本后，就可以用样本来修正已有的先验概率分布，得出后验概率分布，现通过后验概率分布进行各种统计推断。

2.1 双群体

2.1.1 理论推导

在这里插入图片描述

2.1.2 R语言实现

discriminant.bayes.R(双群体最小ECM贝叶斯判别函数)

discriminant.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前面的1/2乘到beta上面去了
    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
}

在程序中，输入变量TrnX1、Trnx2表示X1类、X2类训练样本，其输入格式是数据框，或矩阵（样本按行输入）. $rate=\frac{L(1|2)}{L(2|1)}\cdot \frac{p_2}{p_1}$ ，缺省值为1.TstX是待测样本，其输入格式是数据框，或矩阵（样本按行输入），或向量（一个待测样本）.如果不输入TstX（缺省值），则待测样本为两个训练样本之和，即计算训练样本的回代情况.输入变量var.equal是逻辑变量，var.equal=TRUE表示认为两总体的协方差阵是相同的；否则（缺省值）是不同的.函数的输出是由“1”和“2”构成的的一维矩阵，“1”表示待测样本属于X1类，“2”表示待测样本属于X2类。`

2.1.3 实例分析

在这里插入图片描述
代码：

TrnX1<-matrix(
  c(24.8, 24.1, 26.6, 23.5, 25.5, 27.4,
    -2.0, -2.4, -3.0, -1.9, -2.1, -3.1),
  ncol=2)
TrnX2<-matrix(
  c(22.1, 21.6, 22.0, 22.8, 22.7, 21.5, 22.1, 21.4,
    -0.7, -1.4, -0.8, -1.6, -1.5, -1.0, -1.2, -1.3),
  ncol=2)
path <- 'E:\\桌面文档\\学习文件\\大三\\多元统计\\统计建模R语言书籍'
setwd(path)
source('bayesDiscriminant.R')
discriminant.bayes(TrnX1,TrnX2,rate=8/6,var.equal=TRUE)

结果：

      1 2 3 4 5 6 7 8 9 10 11 12 13 14
blong 1 1 1 2 1 1 2 2 2  2  2  2  2  2

第4号样本被错判。

2.2 多群体

2.2.1 理论推导

在这里插入图片描述

2.2.2 R语言实现

ditinguish.bayes.R(多群体ECM最小准则贝叶斯分类函数)

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(0,nrow=1,ncol=nx,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
}

程序分别考虑了总体协方差阵相同和协方差阵不同的情况.输入变量Trnx表示训练样本，其输入格式是矩阵（样本按行输入），或数据框.TrnG是因子变量，表示训练样本的分类情况.输入变量p是先验概率，缺省值均为1.输入变量TstX是待测样本，其输入格式是矩阵（样本按行输入），或数据框，或向量（一个待测样本）.如果不输入TstX（缺省值），则待测样本为训练样本.输入变量var.equal是逻辑变量，var.equal=TRUE表示认为总体协方差阵是相同的；否则（缺省值）是不同的.函数的输出是由数字构成的的一维矩阵，数字表示相应的类.为了与前面两总体的判别程序兼容，对于二分类问题，也可以按照discriminiant.bayes函数的输入格式输入.

2.2.3 实例分析

在这里插入图片描述
代码：

## 多类别
X <- iris[,1:4]
G <- gl(3,50)
path <- 'E:\\桌面文档\\学习文件\\大三\\多元统计\\统计建模R语言书籍'
setwd(path)
source('bayesDiscriminant2.R')
blong <- distinguish.bayes(X,G)

结果：

      1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
blong 1 1 1 1 1 1 1 1 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
      25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45
blong  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
      46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
blong  1  1  1  1  1  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
      67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87
blong  2  2  3  2  3  2  3  2  2  2  2  3  2  2  2  2  2  3  2  2  2
      88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105
blong  2  2  2  2  2  2  2  2  2  2  2  2   2   3   3   3   3   3
      106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
blong   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3
      121 122 123 124 125 126 127 128 129 130 131 132 133 134 135
blong   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3
      136 137 138 139 140 141 142 143 144 145 146 147 148 149 150
blong   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3

> table(blong,G)
     G
blong  1  2  3
    1 50  0  0
    2  0 45  0
    3  0  5 50

从计算结果可以看出，只有第69、71、73、78、84号样本错判，回代的判别正确率为145/150=96.67%.

3 Fisher判别

Fisher（费歇）判别是按类内方差尽量小，类间方差尽量大的准则来求判别函数的.在这里仅讨论两个总体的判别方法.

3.1 理论推导

在这里插入图片描述

3.2 R语言实现

3.2.1 自己实现

方法1

discriminant.fisher.R(双群体Fisher判别函数)

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

在程序中，输入变量TrnX1、Trnx2表示X1类、X2类训练样本，其输入格式是数据框，或矩阵（样本按行输入）.TstX是待测样本，其输入格式是数据框，或矩阵（样本按行输入），或向量（一个待测样本）.如果不输入TstX（缺省值），则待测样本为两个训练样本之和，即计算训练样本的回代情况.函数的输出是由“1”和“2”构成的的一维矩阵，“1”表示待测样本属于X1类，“2”表示待测样本属于X2类。

有一个需要注意的地方，rep(1,nx) %o% mu不可改为mu，两个算出来的结果是不一样的：

(1)情况1：
> rep(1,nrow(classX1)) %o% mu1
            x1       x2       x3 x4    x5        x6   x7
 [1,] 7.358333 73.66667 1.458333  6 15.25 0.1716667 49.5
 [2,] 7.358333 73.66667 1.458333  6 15.25 0.1716667 49.5
 [3,] 7.358333 73.66667 1.458333  6 15.25 0.1716667 49.5
 [4,] 7.358333 73.66667 1.458333  6 15.25 0.1716667 49.5
 [5,] 7.358333 73.66667 1.458333  6 15.25 0.1716667 49.5
 [6,] 7.358333 73.66667 1.458333  6 15.25 0.1716667 49.5
 [7,] 7.358333 73.66667 1.458333  6 15.25 0.1716667 49.5
 [8,] 7.358333 73.66667 1.458333  6 15.25 0.1716667 49.5
 [9,] 7.358333 73.66667 1.458333  6 15.25 0.1716667 49.5
[10,] 7.358333 73.66667 1.458333  6 15.25 0.1716667 49.5
[11,] 7.358333 73.66667 1.458333  6 15.25 0.1716667 49.5
[12,] 7.358333 73.66667 1.458333  6 15.25 0.1716667 49.5

此时
> classX1 - rep(1,nrow(classX1)) %o% mu1
           x1        x2          x3   x4    x5          x6    x7
1  -0.7583333 -34.66667 -0.45833333  0.0 -9.25 -0.05166667 -29.5
2  -0.7583333 -34.66667 -0.45833333  0.0 -3.25 -0.05166667 -29.5
3  -1.2583333 -26.66667 -0.45833333  0.0 -9.25 -0.09166667 -37.5
4  -1.2583333 -26.66667 -0.45833333  0.0 -3.25 -0.09166667 -37.5
5   1.0416667 -41.66667  0.54166667  1.5  3.75  0.17833333  25.5
6  -0.1583333 -67.66667 -0.45833333  1.0 12.75  0.12833333 -19.5
7   1.0416667  39.33333  2.04166667  0.0  2.75 -0.02166667  25.5
8   0.1416667 -21.66667 -0.45833333  0.0 -3.25 -0.01166667  -9.5
9   0.1416667 -21.66667  2.04166667  1.5 -9.25 -0.01166667  -9.5
10  0.9416667  39.33333 -1.45833333  1.5 19.75 -0.05166667 130.5
11  0.4416667  98.33333 -0.45833333 -2.5 -1.25  0.03833333  -4.5
12  0.4416667  98.33333  0.04166667 -3.0 -0.25  0.03833333  -4.5

(2)情况2：
> classX1 -  mu1
            x1        x2          x3         x4         x5           x6         x7
1   -0.7583333  38.82833  -5.0000000 -67.666667 -43.500000 -15.13000000  18.541667
2  -67.0666667 -10.50000 -14.2500000   4.541667   4.641667  -0.05166667  14.000000
3    4.6416667  39.64167   0.8283333   0.000000 -67.666667 -49.42000000  -3.250000
4    0.1000000 -26.66667 -48.5000000  -9.250000  10.541667  -7.27833333  11.828333
5   -6.8500000  30.54167  -5.3583333   7.328333  13.000000 -73.31666667  25.500000
6    7.0283333   0.00000 -72.6666667 -42.500000  12.750000  -1.15833333  22.641667
7  -41.1000000  97.75000   2.0416667  -1.358333  17.828333  -5.85000000   1.333333
8    0.1416667  51.82833  -5.0000000 -67.666667 -37.500000 -15.09000000  38.541667
9  -66.1666667   2.50000 -11.7500000   6.041667  -1.358333  -0.01166667  34.000000
10   6.8416667 105.64167  -0.1716667   1.500000 -38.666667 -49.38000000 164.750000
11   1.8000000  98.33333 -48.5000000 -11.750000  12.541667  -7.14833333  44.828333
12  -7.4500000 170.54167  -5.8583333   2.828333   9.000000 -73.45666667  -4.500000

换成mu之后结果明显有问题。

方法2

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

  # 方法2
  n1 <- nrow(TrnX1); n2 <- nrow(TrnX2)
  mu1 <- colMeans(TrnX1); mu2 <- colMeans(TrnX2)   # 计算均值向量
  S <- ( (n1-1)*var(TrnX1) + (n2-1)*var(TrnX2) )/(n1+n2-2)    # 混合样本方差
  mu <- (mu1 + mu2)/2   # 这是一个行向量
  W <- (mu2-mu1) %*% solve(S) %*% t(TstX - rep(1,nx) %o% mu)
  # 类别判别
  for (i in 1:nx){
    if (W[i] <= 0)
      blong[i] <- 1
    else
      blong[i] <- 2
  }
  blong
}

3.2.2 借助`MASS`包

MASS包是R语言自带的一个包，不用安装。

classX1$target <- 1
classX2$target <- 2
data <- rbind(classX1,classX2)   # 按行合并数据
library(MASS)
attach(data)
ld <- lda(target ~ x2+x2+x3+x4+x5+x6+x7,prior=c(0.5,0.5))
ld
lp <- predict(ld)
table(lp$class,data$target)
lp$class
detach(data)

3.3 实例分析

使用Fisher判别求解1.1.3中的实例
代码：

# Fisher判别
classX1 <- data.frame(
  x1=c(6.60, 6.60, 6.10, 6.10, 8.40, 7.2, 8.40, 7.50, 7.50, 8.30, 7.80, 7.80),
  x2=c(39.00,39.00, 47.00, 47.00, 32.00, 6.0, 113.00, 52.00, 52.00,113.00,172.00,172.00),
  x3=c(1.00, 1.00, 1.00, 1.00, 2.00, 1.0, 3.50, 1.00,3.50, 0.00, 1.00, 1.50),
  x4=c(6.00, 6.00, 6.00, 6.00, 7.50, 7.0, 6.00, 6.00,7.50, 7.50, 3.50, 3.00),
  x5=c(6.00, 12.00, 6.00, 12.00, 19.00, 28.0, 18.00, 12.00,6.00, 35.00, 14.00, 15.00),
  x6=c(0.12, 0.12, 0.08, 0.08, 0.35, 0.3, 0.15, 0.16,0.16, 0.12, 0.21, 0.21),
  x7=c(20.00,20.00, 12.00, 12.00, 75.00, 30.0, 75.00, 40.00,40.00,180.00, 45.00, 45.00)
)
classX2<-data.frame(
  x1=c(8.40, 8.40, 8.40, 6.3, 7.00, 7.00, 7.00, 8.30,
       8.30, 7.2, 7.2, 7.2, 5.50, 8.40, 8.40, 7.50,
       7.50, 8.30, 8.30, 8.30, 8.30, 7.80, 7.80),
  x2=c(32.0 ,32.00, 32.00, 11.0, 8.00, 8.00, 8.00,161.00,
       161.0, 6.0, 6.0, 6.0, 6.00,113.00,113.00, 52.00,
       52.00, 97.00, 97.00,89.00,56.00,172.00,283.00),
  x3=c(1.00, 2.00, 2.50, 4.5, 4.50, 6.00, 1.50, 1.50,
       0.50, 3.5, 1.0, 1.0, 2.50, 3.50, 3.50, 1.00,
       1.00, 0.00, 2.50, 0.00, 1.50, 1.00, 1.00),
  x4=c(5.00, 9.00, 4.00, 7.5, 4.50, 7.50, 6.00, 4.00,
       2.50, 4.0, 3.0, 6.0, 3.00, 4.50, 4.50, 6.00,
       7.50, 6.00, 6.00, 6.00, 6.00, 3.50, 4.50),
  x5=c(4.00, 10.00, 10.00, 3.0, 9.00, 4.00, 1.00, 4.00,
       1.00, 12.0, 3.0, 5.0, 7.00, 6.00, 8.00, 6.00,
       8.00, 5.00, 5.00,10.00,13.00, 6.00, 6.00),
  x6=c(0.35, 0.35, 0.35, 0.2, 0.25, 0.25, 0.25, 0.08,
       0.08, 0.30, 0.3, 0.3, 0.18, 0.15, 0.15, 0.16,
       0.16, 0.15, 0.15, 0.16, 0.25, 0.21, 0.18),
  x7=c(75.00,75.00, 75.00, 15.0,30.00, 30.00, 30.00, 70.00,
       70.00, 30.0, 30.0, 30.0,18.00, 75.00, 75.00, 40.00,
       40.00,180.00,180.00,180.00,180.00,45.00,45.00)
)

path <- 'E:\\桌面文档\\学习文件\\大三\\多元统计\\统计建模R语言书籍'
setwd(path)
source('fisherDiscriminant.R')
discriminant.fisher(classX1,classX2)

结果：

      1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
blong 1 1 1 1 1 1 1 1 1  1  1  1  2  2  2  2  2  2  2  2  2  2  2  2
      25 26 27 28 29 30 31 32 33 34 35
blong  2  2  2  1  1  2  2  2  2  2  2

将训练样本回代进行判别，有两个点判错，分别是第28，29号样本。对于多类的Fisher判别，其基本原理是相同的，最终也是可以转化为距离判别，这里就不介绍了.

F附录——MASS包的使用

F1 双群体

F1.1 线性判别

使用的数据为

G	x1	x2	G1
1	-1.9	3.2	雨
1	-6.9	0.4	雨
1	5.2	2	雨
1	5	2.5	雨
1	7.3	0	雨
1	6.8	12.7	雨
1	0.9	-5.4	雨
1	-12.5	-2.5	雨
1	1.5	1.3	雨
1	3.8	6.8	雨
2	0.2	6.2	晴
2	-0.1	7.5	晴
2	0.4	14.6	晴
2	2.7	8.3	晴
2	2.1	0.8	晴
2	-4.6	4.3	晴
2	-1.7	10.9	晴
2	-2.6	13.1	晴
2	2.6	12.8	晴
2	-2.8	10	晴

数据可视化：

d6.1 <- read.table("6_1.txt", header = TRUE) # 第一行作为列名
attach(d6.1)
plot(x1,x2)  # 绘制散点图
text(x1,x2,G,adj=-0.5)  # 添加文本标记
detach(d6.1)

在这里插入图片描述
判别分析代码：

library(MASS)
attach(d6.1)
ld <- lda(G~x1+x2)   # 做线性判别分析
ld

Z <- predict(ld)    # 预测
Z
newG <- Z$class
cbind(G,Z$x,newG)
tab <- table(G, newG)  # 生成原始类别与预测类别的列联表
sum(diag(prop.table(tab)))   # prop.table将列联表对应元素全部转为相应的概率
predict(ld, data.frame(x1=8.1,x2=2.0)) # 预测一组新数据对应的输出
detach(d6.1)

结果：

> ld
Call:
lda(G ~ x1 + x2)

Prior probabilities of groups:
  1   2 
0.5 0.5 

Group means:
     x1   x2
1  0.92 2.10
2 -0.38 8.85

Coefficients of linear discriminants:
          LD1
x1 -0.1035305
x2  0.2247957

> Z
$class
 [1] 1 1 1 1 1 2 1 1 1 1 2 2 2 2 1 2 2 2 2 2
Levels: 1 2

$posterior
            1          2
1  0.61625867 0.38374133
2  0.65888888 0.34111112
3  0.89412853 0.10587147
4  0.87143887 0.12856113
5  0.96214822 0.03785178
6  0.17275662 0.82724338
7  0.98442237 0.01557763
8  0.68514398 0.31485602
9  0.85330562 0.14669438
10 0.52789262 0.47210738
11 0.43015877 0.56984123
12 0.30676827 0.69323173
13 0.03336323 0.96663677
14 0.34672296 0.65327704
15 0.88585263 0.11414737
16 0.40213732 0.59786268
17 0.08694507 0.91305493
18 0.03480991 0.96519009
19 0.08934413 0.91065587
20 0.09926372 0.90073628

$x
           LD1
1  -0.28674901
2  -0.39852439
3  -1.29157053
4  -1.15846657
5  -1.95857603
6   0.94809469
7  -2.50987753
8  -0.47066104
9  -1.06586461
10 -0.06760842
11  0.17022402
12  0.49351760
13  2.03780185
14  0.38346871
15 -1.24038077
16  0.24005867
17  1.42347182
18  2.01119984
19  1.40540244
20  1.33503926

> cbind(G,Z$x,newG)
   G         LD1 newG
1  1 -0.28674901    1
2  1 -0.39852439    1
3  1 -1.29157053    1
4  1 -1.15846657    1
5  1 -1.95857603    1
6  1  0.94809469    2
7  1 -2.50987753    1
8  1 -0.47066104    1
9  1 -1.06586461    1
10 1 -0.06760842    1
11 2  0.17022402    2
12 2  0.49351760    2
13 2  2.03780185    2
14 2  0.38346871    2
15 2 -1.24038077    1
16 2  0.24005867    2
17 2  1.42347182    2
18 2  2.01119984    2
19 2  1.40540244    2
20 2  1.33503926    2

> sum(diag(prop.table(tab)))   # prop.table将列联表对应元素全部转为相应的概率
[1] 0.9

> predict(ld, data.frame(x1=8.1,x2=2.0)) # 预测一组新数据对应的输出
$class
[1] 1
Levels: 1 2

$posterior
          1          2
1 0.9327428 0.06725717

$x
        LD1
1 -1.591809

F1.2 二次判别

此时各群体协方差阵不同，使用的数据为：

G	Q	C	P
1	8.3	4	29
1	9.5	7	68
1	8	5	39
1	7.4	7	50
1	8.8	6.5	55
1	9	7.5	58
1	7	6	75
1	9.2	8	82
1	8	7	67
1	7.6	9	90
1	7.2	8.5	86
1	6.4	7	53
1	7.3	5	48
2	6	2	20
2	6.4	4	39
2	6.8	5	48
2	5.2	3	29
2	5.8	3.5	32
2	5.5	4	34
2	6	4.5	36

数据可视化为：

d6.2 <- read.table("6_2.txt",header = TRUE)
attach(d6.2)
plot(Q,C)
text(Q,C,G,adj=-0.5)

plot(Q,P)
text(Q,P,G,adj=-0.5)

plot(C,P)
text(C,P,G,adj=-0.5)

在这里插入图片描述

判别分析代码：

library(MASS)
qd <- qda(G~Q+C+P)
qd

Z <- predict(qd) # 预测
newG <- Z$class
cbind(G,newG)  # 按列合并

predict(qd, data.frame(Q=8,C=7.5,P=65))   # 预测一组新数据

> qd
Call:
qda(G ~ Q + C + P)

Prior probabilities of groups:
   1    2 
0.65 0.35 

Group means:
         Q        C        P
1 7.976923 6.730769 61.53846
2 5.957143 3.714286 34.00000

> cbind(G,newG)  # 按列合并
      G newG
 [1,] 1    1
 [2,] 1    1
 [3,] 1    1
 [4,] 1    1
 [5,] 1    1
 [6,] 1    1
 [7,] 1    1
 [8,] 1    1
 [9,] 1    1
[10,] 1    1
[11,] 1    1
[12,] 1    1
[13,] 1    1
[14,] 2    2
[15,] 2    2
[16,] 2    2
[17,] 2    2
[18,] 2    2
[19,] 2    2
[20,] 2    2

> predict(qd, data.frame(Q=8,C=7.5,P=65))   # 预测一组新数据
$class
[1] 1
Levels: 1 2

$posterior
          1            2
1 0.9998462 0.0001537705

F2 多总体

F2.1 线性判别

此时协方差阵相等，使用的数据为

G	Q	C	P
1	8.3	4	29
1	9.5	7	68
1	8	5	39
1	7.4	7	50
1	8.8	6.5	55
2	9	7.5	58
2	7	6	75
2	9.2	8	82
2	8	7	67
2	7.6	9	90
2	7.2	8.5	86
2	6.4	7	53
2	7.3	5	48
3	6	2	20
3	6.4	4	39
3	6.8	5	48
3	5.2	3	29
3	5.8	3.5	32
3	5.5	4	34
3	6	4.5	36

数据可视化为

d6.3 <- read.table("6_3.txt",header = TRUE)
attach(d6.3)
plot(Q,C)
text(Q,C,G,adj=-0.5)

plot(Q,P)
text(Q,P,G,adj=-0.5)

plot(C,P)
text(C,P,G,adj=-0.5)

在这里插入图片描述

线性判别分析为

ld <- lda(G~Q+C+P)
ld

Z <- predict(ld)
newG <- Z$class
cbind(G,Z$x,newG)

tab <- table(G,newG)
diag(prop.table(tab,1))

sum(diag(prop.table(tab)))

plot(Z$x, main="陈炜")
text(Z$x[,1], Z$x[,2],G, adj=-0.5,cex=0.75)

predict(ld,data.frame(Q=8,C=7.5,P=65))

结果：

> ld
Call:
lda(G ~ Q + C + P)

Prior probabilities of groups:
   1    2    3 
0.25 0.40 0.35 

Group means:
         Q        C      P
1 8.400000 5.900000 48.200
2 7.712500 7.250000 69.875
3 5.957143 3.714286 34.000

Coefficients of linear discriminants:
          LD1         LD2
Q -0.81173396  0.88406311
C -0.63090549  0.20134565
P  0.01579385 -0.08775636

Proportion of trace:
   LD1    LD2 
0.7403 0.2597 

> cbind(G,Z$x,newG)
   G        LD1          LD2 newG
1  1 -0.1409984  2.582951755    1
2  1 -2.3918356  0.825366275    1
3  1 -0.3704452  1.641514840    1
4  1 -0.9714835  0.548448277    1
5  1 -1.7134891  1.246681993    1
6  2 -2.4593598  1.361571174    1
7  2  0.3789617 -2.200431689    2
8  2 -2.5581070 -0.467096091    2
9  2 -1.1900285 -0.412972027    2
10 2 -1.7638874 -2.382302324    2
11 2 -1.1869165 -2.485574940    2
12 2 -0.1123680 -0.598883922    2
13 2  0.3399132  0.232863397    3
14 3  2.8456561  0.936722573    3
15 3  1.5592346  0.025668216    3
16 3  0.7457802 -0.209168159    3
17 3  3.0062824 -0.358989534    3
18 3  2.2511708  0.008852067    3
19 3  2.2108260 -0.331206768    3
20 3  1.5210939  0.035984885    3

> diag(prop.table(tab,1))
   1    2    3 
1.00 0.75 1.00 
> 
> sum(diag(prop.table(tab)))
[1] 0.9

> predict(ld,data.frame(Q=8,C=7.5,P=65))
$class
[1] 2
Levels: 1 2 3

$posterior
          1        2           3
1 0.2114514 0.786773 0.001775594

$x
        LD1        LD2
1 -1.537069 -0.1367865

结果的可视化：

plot(Z$x)
text(Z$x[,1], Z$x[,2],G, adj=-0.5,cex=0.75)

在这里插入图片描述

F2.2 二次判别

使用F2.1中的数据，判别分析的代码为：

qd <- qda(G~Q+C+P)
Z <- predict(qd)
newG <- Z$class
cbind(G,newG)
sum(diag(prop.table(tab)))
predict(qd,data.frame(Q=8,C=7.5,P=65))

结果：

> cbind(G,newG)
      G newG
 [1,] 1    1
 [2,] 1    1
 [3,] 1    1
 [4,] 1    1
 [5,] 1    1
 [6,] 2    2
 [7,] 2    2
 [8,] 2    2
 [9,] 2    2
[10,] 2    2
[11,] 2    2
[12,] 2    2
[13,] 2    3
[14,] 3    3
[15,] 3    3
[16,] 3    3
[17,] 3    3
[18,] 3    3
[19,] 3    3
[20,] 3    3
> sum(diag(prop.table(tab)))
[1] 0.9
> predict(qd,data.frame(Q=8,C=7.5,P=65))
$class
[1] 2
Levels: 1 2 3

$posterior
            1         2            3
1 0.008221165 0.9915392 0.0002396287

F2.3 贝叶斯判别

使用F2.1中的数据，由于在协方差阵相同的情况下，贝叶斯判别的判别函数具有与距离判别判别函数相似的形式，差别在于多了先验概率和误判代价，因此仍然可以使用线性判别分析的函数做贝叶斯判别。不考虑误判代价，则判别分析的代码为：

ld1 <- lda(G~Q+C+P,prior=rep(1,3)/3)   # 先验概率相等的贝叶斯模型
ld1

ld2 <- lda(G~Q+C+P,prior=c(5,8,7)/20)   # 先验概率不等的贝叶斯模型
ld2

Z1 <- predict(ld1)
cbind(G,Z2$x,newG=Z$class)
table(G,Z1$class)
round(Z1$post,3)   # 保留三位有效数字

predict(ld1,data.frame(Q=8,C=7.5,P=65))
predict(ld2,data.frame(Q=8,C=7.5,P=65))

结果：

> ld1
Call:
lda(G ~ Q + C + P, prior = rep(1, 3)/3)

Prior probabilities of groups:
        1         2         3 
0.3333333 0.3333333 0.3333333 

Group means:
         Q        C      P
1 8.400000 5.900000 48.200
2 7.712500 7.250000 69.875
3 5.957143 3.714286 34.000

Coefficients of linear discriminants:
          LD1         LD2
Q -0.92307369  0.76708185
C -0.65222524  0.11482179
P  0.02743244 -0.08484154

Proportion of trace:
   LD1    LD2 
0.7259 0.2741 

> ld2
Call:
lda(G ~ Q + C + P, prior = c(5, 8, 7)/20)

Prior probabilities of groups:
   1    2    3 
0.25 0.40 0.35 

Group means:
         Q        C      P
1 8.400000 5.900000 48.200
2 7.712500 7.250000 69.875
3 5.957143 3.714286 34.000

Coefficients of linear discriminants:
          LD1         LD2
Q -0.81173396  0.88406311
C -0.63090549  0.20134565
P  0.01579385 -0.08775636

Proportion of trace:
   LD1    LD2 
0.7403 0.2597 

> table(G,Z1$class)
   
G   1 2 3
  1 5 0 0
  2 1 6 1
  3 0 0 7
  
> round(Z1$post,3)   # 保留三位有效数字
       1     2     3
1  0.983 0.006 0.012
2  0.794 0.206 0.000
3  0.937 0.043 0.020
4  0.654 0.337 0.009
5  0.905 0.094 0.000
6  0.928 0.072 0.000
7  0.003 0.863 0.133
8  0.177 0.822 0.000
9  0.185 0.811 0.005
10 0.003 0.997 0.000
11 0.002 0.997 0.001
12 0.111 0.780 0.109
13 0.292 0.325 0.383
14 0.001 0.000 0.999
15 0.012 0.023 0.965
16 0.079 0.243 0.678
17 0.000 0.000 1.000
18 0.001 0.003 0.996
19 0.001 0.004 0.995
20 0.014 0.025 0.961

> predict(ld1,data.frame(Q=8,C=7.5,P=65))
$class
[1] 2
Levels: 1 2 3

$posterior
         1         2           3
1 0.300164 0.6980356 0.001800378

$x
        LD1        LD2
1 -1.426693 -0.5046594

> predict(ld2,data.frame(Q=8,C=7.5,P=65))
$class
[1] 2
Levels: 1 2 3

$posterior
          1        2           3
1 0.2114514 0.786773 0.001775594

$x
        LD1        LD2
1 -1.537069 -0.1367865

案例分析

# ==================判别分析===================== #
# 两群体同协方差
path <- "E:\\桌面文档\\学习文件\\大三\\多元统计\\实验\\第一次上机\\第一次上机课的数据\\第一次上机课的数据\\判别分析\\陈炜\\"
setwd(path)  # 设置工作路径

# ==========================================案例1========================================#
# =======线性判别分析============= #
d6.1 <- read.table("6_1.txt", header = TRUE) # 第一行作为列名
attach(d6.1)
plot(x1,x2)  # 绘制散点图
text(x1,x2,G,adj=-0.5)  # 添加文本标记
detach(d6.1)

library(MASS)
attach(d6.1)
ld <- lda(G~x1+x2)   # 做线性判别分析
ld

Z <- predict(ld)    # 预测
Z
newG <- Z$class
cbind(G,Z$x,newG)
tab <- table(G, newG)  # 生成原始类别与预测类别的列联表
sum(diag(prop.table(tab)))   # prop.table将列联表对应元素全部转为相应的概率
predict(ld, data.frame(x1=8.1,x2=2.0)) # 预测一组新数据对应的输出
detach(d6.1)

# =========================================案例2===============================#
# ========二次判别分析============ #
d6.2 <- read.table("6_2.txt",header = TRUE)
attach(d6.2)
plot(Q,C)
text(Q,C,G,adj=-0.5)

plot(Q,P)
text(Q,P,G,adj=-0.5)

plot(C,P)
text(C,P,G,adj=-0.5)

library(MASS)
qd <- qda(G~Q+C+P)
qd

Z <- predict(qd) # 预测
newG <- Z$class
cbind(G,newG)  # 按列合并

predict(qd, data.frame(Q=8,C=7.5,P=65))   # 预测一组新数据

# =========一次线性判别==================#
ld <- lda(G~.,data = d6.2)   # 假定协方差相等
ld

W <- predict(ld)  # 预测值
cbind(G,W$x,newG=W$class)
predict(ld,data.frame(Q=8,C=7.5,P=65))

detach(d6.2)

# =========================================案例3===============================#
### 多总体距离判别
## (等协方差的时候仍然是线性判别)
d6.3 <- read.table("6_3.txt",header = TRUE)
attach(d6.3)
plot(Q,C)
text(Q,C,G,adj=-0.5)

plot(Q,P)
text(Q,P,G,adj=-0.5)

plot(C,P)
text(C,P,G,adj=-0.5)

# -------------------------------------以下部分还没放进word
ld <- lda(G~Q+C+P)
ld

Z <- predict(ld)
newG <- Z$class
cbind(G,Z$x,newG)

tab <- table(G,newG)
diag(prop.table(tab,1))

sum(diag(prop.table(tab)))

plot(Z$x)
text(Z$x[,1], Z$x[,2],G, adj=-0.5,cex=0.75)

predict(ld,data.frame(Q=8,C=7.5,P=65))

## 二次判别，异方差
qd <- qda(G~Q+C+P)
Z <- predict(qd)
newG <- Z$class
cbind(G,newG)
sum(diag(prop.table(tab)))
predict(qd,data.frame(Q=8,C=7.5,P=65))

## 贝叶斯判别
ld1 <- lda(G~Q+C+P,prior=rep(1,3)/3)   # 先验概率相等的贝叶斯模型
ld1

ld2 <- lda(G~Q+C+P,prior=c(5,8,7)/20)   # 先验概率不等的贝叶斯模型
ld2

Z1 <- predict(ld1)
cbind(G,Z2$x,newG=Z$class)
table(G,Z1$class)
round(Z1$post,3)   # 保留三位有效数字

predict(ld1,data.frame(Q=8,C=7.5,P=65))
predict(ld2,data.frame(Q=8,C=7.5,P=65))

detach(d6.3)

### 案例分析
Case5 <- read.table('6_case.txt',header = T)
head(Case5,30)    # 显示前30行数据

attach(Case5)
plot(Case5[,2:5], gap=0, main="陈炜")    # 绘制两两散点图
ld <- lda(G~.,data=Case5)
ld

plot(ld,main="陈炜")

Zld <- predict(ld)
data.frame(Case5$G,Zld$class,round(Zld$x,3))
addmargins(table(Case5$G,Zld$class))

qd <- qda(G~.,data=Case5)  # 二次判别
qd

Zqd <- predict(qd)
addmargins(table(Case5$G,Zqd$class))

detach(Case5)

###########################################主成分分析
path <- "E:\\桌面文档\\学习文件\\大三\\多元统计\\实验\\第一次上机\\第一次上机课的数据\\第一次上机课的数据\\主成分分析\\陈炜\\"
setwd(path)
# (1)计算相关矩阵
X <- read.table("d7_2.txt",header=T)
cor(X)

# (2)求相关矩阵的特征值和主成分负荷
PCA <- princomp(X,cor=T)  # cor=T的意思是使用相关矩阵
PCA

PCA$loadings   # 载荷矩阵
summary(PCA,loading=TRUE) # 显示载荷矩阵

# (3)确定主成分
screeplot(PCA,type='lines',main="陈炜")   # 使用直线绘制碎石图

# (4)主成分得分
PCA$scores[,1:2]   # 其实就是每一个原始变量对应的主成分值


## 案例
Case8 <- read.table('case8.txt', header = T)
PC <- princomp(Case8,cor=T)
summary(PC)
m <- 2   # 选两个主成分(只有前两个特征值大于1)
PC$loadings[1:m]
princomp.rank(PC,m)  # 主成分排序
princomp.rank(PC,m,plot=T)

三只佩奇不结义

关注

55
点赞
踩
426

收藏

觉得还不错? 一键收藏
打赏
7
评论
R语言实战——距离判别、贝叶斯判别、Fisher判别理论详细推导与R语言实现

前言判别分析是用以判别个体所属群体的一种统计方法，它产生于20世纪30年代，近年来，在许多现代自然科学的各个分支和技术部门中，得到广泛的应用.例如，利用计算机对一个人是否有心脏病进行诊断时，可以取一批没有心脏病的人，测其p个指标的数据，然后再取一批已知患有心胜病的人，同样也测得p个相同指标的数据，利用这些数据建一个判别函数，并求出相应的临界值，这时对于需要进行诊断的人，也同样测其p个指标的数据...
复制链接

扫一扫