R语言中的多元统计之判别分析报告

前言判别分析（discriminant analysis）是多元统计分析中较为成熟的一种分类方法，它的核心思想是“分类与判断”，即根据已知类别的样本所提供的信息，总结出分类的规律性，并建立好判别公式和判别准则，在此基础上，新的样本点将按照此准则判断其所属类型。

例如，根据一年甚至更长时间的每天的湿度差及压差，我们可以建立一个用于判别是否会下雨的模型，当我们获取到某一天（建立模型以外的数据）的湿度差及压差后，使用已建立好的模型，就可以得出这一天是否会下雨的判断。

根据判别的组数来区分，判别分析可以分为两组判别和多组判别。

接下来，我们将学习三种常见的判别分析方法，分别是：•距离判别•Bayes判别•Fisher判别一、距离判别基本理论假设存在两个总体和，另有为一个维的样本值，计算得到该样本到两个总体的距离和,如果大于，则认为样本属于总体，反之样本则属于总体；若等于，则该样本待判。

这就是距离判别法的基本思想。

在距离判别法中，最核心的问题在于距离的计算，一般情况下我们最常用的是欧式距离，但由于该方法在计算多个总体之间的距离时并不考虑方差的影响，而马氏距离不受指标量纲及指标间相关性的影响，弥补了欧式距离在这方面的缺点，其计算公式如下：，为总体之间的协方差矩阵二、距离判别的R实现（训练样本）首先我们导入数据# 读取SAS数据> library(sas7bdat)> data1 <- read.sas7bdat('disl01.sas7bdat')# 截取所需列数据，用于计算马氏距离> testdata <- data1[2:5]> head(testdata,3)X1 X2 X3 X41 -0.45 -0.41 1.09 0.452 -0.56 -0.31 1.51 0.163 0.06 0.02 1.01 0.40# 计算列均值> colM <- colMeans(testdata)> colMX1 X2 X3 X4 0.096304348 -0.006956522 2.033478261 0.431739130 # 计算矩阵的协方差> cov_test <- cov(testdata)> cov_testX1 X2 X3 X4 X1 0.068183816 0.027767053 0.14996870 -0.002566763 X2 0.027767053 0.015363865 0.05878251 0.001252367X3 0.149968696 0.058782512 1.01309874 0.028607150X4 -0.002566763 0.001252367 0.02860715 0.033912464# 样本的马氏距离计算> distance <- mahalanobis(testdata,colM,cov_test)> head(distance,5)[1] 12.726465 11.224681 1.692702 1.347885 2.369820这样，我们得到了距离判别中最关键的马氏距离值，在此基础上就可以进行进一步的判别分析了。

不过我们介绍一个R的第三方包WMDB，该包的wmd()函数可以简化我们的距离判别过程，函数将输出样本的分类判别结果、错判的样本信息以及判别分析的准确度。

> library(WMDB)> head(data1,3)A X1 X2 X3 X41 1 -0.45 -0.41 1.09 0.452 1 -0.56 -0.31 1.51 0.163 1 0.06 0.02 1.01 0.40# 提取原始数据集的A列生成样品的已知类别> testdata_group <- data1$A# 转换为因子变量，用于wmd()函数中> testdata_group <- as.factor(testdata_group)> wmd(testdata,testdata_group)1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27blong 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 1 1 1 2 1 1 1 1 2 2 228 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46blong 2 2 2 2 2 2 1 2 2 2 1 1 1 1 1 2 1 2 2 [1] "num of wrong judgement"[1] 15 16 20 22 23 24 34 38 39 40 41 42 44[1] "samples divided to"[1] 2 2 2 1 1 1 1 1 1 1 1 1 1[1] "samples actually belongs to"[1] 1 1 1 2 2 2 2 2 2 2 2 2 2Levels: 1 2[1] "percent of right judgement"[1] 0.7173913由分析结果可知，根据已知分类的训练样品建立的判别规则，重新应用于训练样品后，出现了13个错判样品，拥有71.7%的准确度。

三、距离判别的R实现（测试样本）接着，当我们获取到未分类的新样本数据时，使用wmd()函数，在训练样本的基础上进行这些数据的距离判别# 导入数据，一共10个样本> data2 <- read.sas7bdat('disldp01.sas7bdat')# 截取所需列数据> newtestdata <- data2[1:4]# 进行判别分析> wmd(testdata,testdata_group,TstX = newtestdata)1 2 3 4 5 6 7 8 9 10blong 1 1 1 1 1 1 2 2 2 1根据马氏距离判别分析得到的结果，10个待判样品中，第一类7个，第二类3个。

距离判别方法简单实用，它只要求知道总体的数字特征，而不涉及总体的分布，当总体均值和协方差未知时，就用样本的均值和协方差矩阵来估计，因此距离判别没有考虑到每个总体出现的机会大小，即先验概率，没有考虑到错判的损失。

因此，我们进一步学习贝叶斯判别法。

一、贝叶斯判别基本理论贝叶斯判别法的前提是假定我们已经对所要分析的数据有所了解（比如数据服从什么分别，各个类别的先验概率等），根据各个类别的先验概率求得新样本属于某类的后验概率。

该算法应用到经典的贝叶斯公式，该公式为：假设有两个总体和，分别具有概率密度函数和，并且根据以往的统计分析，两个总体各自出现的先验概率为和，当一个样本发生时，求该样本属于某一类的概率，计算公式为：这样，我们得到了该样本属于两类总体的概率，分别为和，属于哪一类总体的概率值大，我们则将样本划分到该类中。

二、贝叶斯判别的R实现在R中，我们使用klaR包中的NaiveBayes()函数实现贝叶斯判别分析，函数调用公式如下：> NaiveBayes(formula, data, ..., subset, na.action = na.pass)# formula指定参与模型计算的变量，以公式形式给出，类似于y=x1+x2+x3# na.action指定缺失值的处理方法，默认情况下不将缺失值纳入模型计算，也不会发生报错信息，当设为“na.omit”时则会删除含有缺失值的样本# 数据准备，使用R内置数据集iris# 通过抽样建立训练样本(70%)和测试样本(30%)> index <- sample(2,size = nrow(iris),replace = TRUE,prob = c(0.7,0.3))> train_data <- iris[index == 1,]> test_data <- iris[index == 2,]# 载入所用包> library(klaR)# 构建贝叶斯模型> Bayes_model <- NaiveBayes(Species ~ ., data = train_data)# 进行预测> Bayes_model_pre <- predict(Bayes_model, newdata = test_data[,1:4]) # 生成实际与预判交叉表> table(test_data$Species,Bayes_model_pre$class)setosa versicolor virginicasetosa 20 0 0versicolor 0 17 0virginica 0 3 7从上表生成的交叉表中，我们可以看到在该模型中错判了3个。

# 生成预判精度> sum(diag(table(test_data$Species,Bayes_model_pre$class)))+ / sum(table(test_data$Species,Bayes_model_pre$class))[1] 0.9361702三、Fisher判别基本理论Fisher判别法的基本思想是“投影”，将组维的数据向低维空间投影，使其投影的组与组之间的方差尽可能的大，组内的方差尽可能的小。

因此，Fisher判别法的重点就是选择适当的“投影轴”。

判别函数为，接下来我们以两类总体举例。

首先我们将样本点投影到一维空间，旋转坐标轴至总体单位尽可能分开的方向，此时分类变量被简化为一个，判别函数；如果划分的效果不理想，可以考虑投影到二维空间（）,以此类推。

上图为二维空间的Fisher判别，从图中可以看到，无论我们把总体和投影到还是轴，都不能很好的把两类总体区分出来。

为此，我们需要寻找一条合适的投影线，使得两类总体向该线投影后的区分程度达到最大，线性判别函数即为该投影线的表达形式（这里我们仅介绍Fisher判别的基本原理，不涉及参数的具体推导和求解，这些都可用R程序求得）。

四、Fisher判别的R实现在R中，我们使用MASS包中的lda()函数实现Fisher判别分析，函数调用公式如下：> lda(formula, data, ..., subset, na.action)# formula:指定参与模型计算的变量，以公式形式给出，类似于y=x1+x2+x3 # na.action:指定缺失值的处理方法，默认情况下，缺失值的存在使算法无法运行，当设置为“na.omit”时则会删除含有缺失值的样本# 数据准备，使用R内置数据集iris# 通过抽样建立训练样本(70%)和测试样本(30%)> index <- sample(2,size = nrow(iris),replace = TRUE, prob = c(0.7,0.3)) > train_data <- iris[index == 1,]> test_data <- iris[index == 2,]# 载入所用包> library(MASS)# 构建Fisher判别模型> fisher_model <- lda(Species~., data = train_data)# 进行预测> fisher_model_pre <- predict(fisher_model, newdata = test_data[,1:4]) # 生成实际与预判交叉表> table(test_data$Species,fisher_model_pre$class)setosa versicolor virginicasetosa 20 0 0versicolor 0 14 1virginica 0 0 18# 生成预判精度> sum(diag(table(test_data$Species,fisher_model_pre$class)))+ / sum(table(test_data$Species,fisher_model_pre$class))[1] 0.9811321五、Fisher判别进阶——非线性判别在判别分析的实际应用中，对复杂的数据使用线性判别可能无法得到理想的效果。