R语言判别分析

自己整理编写的R语言常用数据分析模型的模板，原文件为Rmd格式，直接复制粘贴过来，作为个人学习笔记保存和分享。部分参考薛毅的《统计建模与R软件》和《R语言实战》

本文中分三个方法介绍判别分析，Bayes判别，距离判别，Fisher判别。前两种判别方法都要考虑两个、或多个总体协方差（这里是算方差，方差是协方差的一种）相等或不等的情况，由var.equal=的逻辑参数表示，默认是FALSE，表示认为两总体协方差不等。用样本的协方差可以估计总体的协方差。在Bayes方法中我们把相等和不等的两个结果都列了出来，距离判别里我们默认两总体协方差不等。事实上，一般使用时，我们都以两总体的协方差不等作为标准来进行后续计算。

协方差计算公式： S(xy) = N,i ∑ [(xi-x均)(yi-y均)]/(N-1) ; 这里要计算的总体协方差就是方差，为 S = N,i ∑ [(xi-x均)*(xi-x均)^T]/(N-1)

1. Bayes判别

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

* 1.1 两总体判别的Bayes判别程序*

在程序中，输入变量TrnX1，TrnX2表示X1类，X2类训练样本，其输入格式是数据框，或矩阵（样本按行输入）.rate= L(1|2) / L(2|1)*p2/p1 ， 缺省值为1. TstX是待测样本，其输入格式是数据框，或矩阵（样本按行输入），或向量（一个待测样本）.如果不输入TstX（缺省值），则待测样本为两个训练样本之和，即计算训练样本的回代情况.输入变量var.equal是逻辑变量，var.equal=TRUE表示认为两总体的协方差阵是相同的；否则（缺省值）是不同的. 函数的输出是由”1”,”2”构成的一维矩阵，”1”代表待测样本属于X1类，”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
}

* 1.2 例. 用1.1中Bayes判别程序对冠心病和正常人进行判别分类*

为研究舒张期血压与血浆胆固醇对冠心病的作用，调查了50~59岁的女冠心病人15例和正常人16例。他们的舒张压（x1）与血浆胆固醇（x2）列在下表中。试用判别分析法建立判别冠心病与正常人的判别函数

#它们的先验概率pi分别用15/31、16/31来估计
#载入数据
TrnX1<-matrix(
   c(9.86, 13.33, 14.66, 9.33, 12.80, 10.66, 10.66, 13.33, 13.33, 13.33, 12.00, 14.66, 13.33, 12.80, 13.33, 
     5.18, 3.73, 3.89, 7.10, 5.49, 4.09, 4.45, 3.63, 5.96, 5.70, 6.19, 4.01, 4.01, 3.63, 5.96),
   ncol=2)
TrnX2<-matrix(
   c(10.66, 12.53, 13.33, 9.33, 10.66, 10.66, 9.33, 10.66, 10.66, 10.66, 10.40, 9.33, 10.66, 10.66, 11.20, 9.33,
     2.07, 4.45, 3.06, 3.94, 4.45, 4.92, 3.68, 2.77, 3.21, 5.02, 3.94, 4.92, 2.69, 2.43, 3.42, 3.63),
   ncol=2)
#输入待测样本TrnX
TstX<-matrix(
   c(9.06,13.00,12.66,9.00,12.12,11.66,12.11,12.63,8.33,11.12,
     5.68,3.43,2.82,6.86,5.19,2.17,4.12,3.26,2.91,4.22 ),
   ncol=2)
#载入两总体的贝叶斯判别函数    注 把贝叶斯判别函数存在了wd中
source("discriminiant.bayes.R")
#### 总体协方差阵相同
discriminiant.bayes(TrnX1, TrnX2,rate=16/15,var.equal=TRUE)    #rate=L(1|2)/L(2|1)* p2/p1
discriminiant.bayes(TrnX1, TrnX2,TstX,rate=16/15,var.equal=TRUE)
#### 总体协方差阵不同
discriminiant.bayes(TrnX1, TrnX2,rate=16/15