Mahalanobis 距离
定义：马氏距离(Mahalanobis distance)是由印度统计学家马哈拉诺比斯(P. C. Mahalanobis)提出的，表示数据的协方差距离。它是一种有效的计算两个未知样本集的相似度的方法。与欧氏距离不同的是它考虑到各种特性之间的联系并且是尺度无关，即独立于测量尺度。 对于一个均值为 μ，协方差矩阵为 Σ 的多变量向量，其马氏距离公式为：
距离判别
定义：是由训练样品得出每个分类的重心坐标，然后对新样品求出它们离各个类别重心的距离远近，从而归入离得最近的类。距离判别的特点是直观、简单，适合于对自变量均为连续变量的情况下进行分类。 本实验分别讨论两总体协方差相同和不同的情况。
两总体样本 X1 和 X2 的协方差相同，即 μ1、μ2 不等 ， Σ1 = Σ2 = Σ
判别准则公式:
两总体样本 X1 和 X2 的协方差不同，即 μ1、μ2 不等 ， Σ1、Σ2、Σ 均不等
判别函数:
Bayes 判别
定义：贝叶斯判别是根据最小风险代价判决或最大似然比判决，是根据贝叶斯准则进行判别分析的一种多元统计分析法。基本思想是：设有两个总体，它们的先验概率分别为 p1 、 p2，各总体的密度函数为 f1(x) 、f2(x) ，在观测到一个样本 x 的情况下，可以证明极小化平均误判损失函数的划分区域 R1 和 R2 来作为Bayes 判别的判别准则。
判别准则：
样本协方差相同：
样本协方差不同：
Fisher 判别
定义：费歇尔（Fisher）判别是一种先进行高维向低位投影，再根据距离判别的一种方法。借助方差分析的思想构造判别函数（相当于一种投影），使组间区别最大、组内离差最小，然后代入新样本数据，将其与判别临界值比较以确定应判为至哪一总体。
构造判别准则公式：
数据准备
本实验以某中小企业的破产模型作为数据来源，选定了 4 个经济指标，包括：
X1 总负债率（现金收益/总负债）
X2 收益性指标（纯收入/总财产）
X3 短期支付能力（流动资产/流动负债）
X4 生产效率性指标（流动资产/纯销售额）
其中，数据包含 17 个破产企业（1 类）和 21 个正常运行企业（2 类），以及 8 个待判定的企业。
读取数据，并将表格数据分为训练样本和测试样本，trainX1 代表第一类训练样本，trainX2 代表第二类训练样本，testX 为测试样本。
> options(stringsAsFactors=F)
> bankruptcy <- readLines("http://labfile.oss.aliyuncs.com/courses/916/bankruptcy.csv")
> bankruptcy <- unlist(strsplit(bankruptcy, split=","))
> bankruptcy <- as.data.frame(matrix(bankruptcy, ncol=5, byrow=T))
> colnames(bankruptcy) <- bankruptcy[1,]
> bankruptcy <- bankruptcy[-1,]
> bankruptcy <- as.data.frame(sapply(bankruptcy, as.numeric))
Warning message:
In lapply(X = X, FUN = FUN, ...) : NAs introduced by coercion
>
> trainX1 <- data.frame(x1=bankruptcy[1:17, 1], x2=bankruptcy[1:17, 2], x3=bankruptcy[1:17, 3], x4=bankruptcy[1:17, 4])
> trainX2 <- data.frame(x1=bankruptcy[18:38, 1], x2=bankruptcy[18:38, 2], x3=bankruptcy[18:38, 3], x4=bankruptcy[18:38, 4])
> testX <- data.frame(x1=bankruptcy[39:46, 1], x2=bankruptcy[39:46, 2], x3=bankruptcy[39:46, 3], x4=bankruptcy[39:46, 4])
距离判别
将距离判别中的算法编写成名为 discriminiant.distance 的函数，函数中，输入变量 TrnX1,TrnX2 分别表示 X1 类和 X2 类的训练样本，格式可以是矩阵、数据框，输入变量 TstX 是测试样本，格式是矩阵、数据框或向量，var.equal 是逻辑变量，TRUE 表示协方差阵相同，默认为不同，函数输出以 “1” 和“2” 构成的一维矩阵，“1” 代表测试样本属于X1 类，“2” 代表属于 X2 类。
距离判别主要用到了 mahalanobis() 函数，格式为：
mahalanobis(data,center,cov)
data:样本数据构成的向量或矩阵
center：样本中心
cov：样本协方差
> discriminiant.distance<-function(TrnX1,TrnX2,TstX=NULL,var.equal=F){
+ if (is.null(TstX)==T)
+ TstX <- rbind(TrnX1,TrnX2)
+ if (is.vector(TstX)==T)
+ TstX<- t(as.matrix(TstX))
+ else if (is.matrix(TstX)!=T)
+ TstX<- as.matrix(TstX)
+ if (is.matrix(TrnX1)!=T)
+ TrnX1 <- as.matrix(TrnX1)
+ if (is.matrix(TrnX2)!=T)
+ TrnX2 <- as.matrix(TrnX2)
+ nx <- nrow(TstX)
+ blong <- matrix(rep(0,nx),nrow=1,byrow=T,dimnames=list("blong",1:nx))
+ mu1<-colMeans(TrnX1)
+ mu2<-colMeans(TrnX2)
+
+ if (var.equal==T || 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.distance(trainX1, trainX2)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
blong 1 1 1 1 1 1 1 1 1 1 2 2 1 1 1 2 1 2 2 2
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37
blong 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
38
blong 2
通过与原始样本对比，可以看出样本 11 12 和 16 被判错
> discriminiant.distance(trainX1, trainX2, var.equal=T)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
blong 1 1 1 1 1 1 1 1 1 1 2 2 1 1 1 1 1 2 1 2
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37
blong 2 2 2 2 1 2 1 2 2 1 2 1 2 2 2 2 2
38
blong 2
样本 11 12 19 25 27 30 32 被判错
> discriminiant.distance(trainX1, trainX2, testX)
1 2 3 4 5 6 7 8
blong 1 1 1 1 2 2 2 2
Bayes 判别
将Bayes 判别中的算法编写成名为 discriminiant.bayes 的函数，函数中，输入变量 TrnX1,TrnX2 分别表示 X1 类和 X2 类的训练样本，格式可以是矩阵、数据框，输入变量 TstX 是测试样本，格式是矩阵、数据框或向量，rate 是损失比和先验概率的乘积，默认为 1 ，var.equal 是逻辑变量，TRUE 表示协方差阵相同，默认为不同，函数输出以 “1” 和“2” 构成的一维矩阵，“1” 代表测试样本 属于X1 类，“2” 代表属于 X2 类。
> discriminiant.bayes <- function(TrnX1, TrnX2, rate=1, TstX=NULL, var.equal=F){
+ if (is.null(TstX)==T)
+ TstX <- rbind(TrnX1, TrnX2)
+ if (is.vector(TstX)==T)
+ TstX <- t(as.matrix(TstX))
+ else if (is.matrix(TstX)!=T)
+ TstX <- as.matrix(TstX)
+ if (is.matrix(TrnX1)!=T)
+ TrnX1 <- as.matrix(TrnX1)
+ if (is.matrix(TrnX2)!=T)
+ TrnX2 <- as.matrix(TrnX2)
+
+ nx <- nrow(TstX)
+ blong <- matrix(rep(0, nx), nrow=1, byrow=T, dimnames=list("blong", 1:nx))
+ mu1 <- colMeans(TrnX1)
+ mu2 <- colMeans(TrnX2)
+
+ if (var.equal==T || 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.bayes(trainX1, trainX2, rate=21/17)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
blong 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 2 1 2
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37
blong 2 2 2 2 1 2 1 1 2 2 2 2 2 2 1 2 2
38
blong 2
样本 12 19 25 27 28 35 被判错
> discriminiant.bayes(trainX1, trainX2, rate=21/17, var.equal=T)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
blong 1 1 1 1 1 1 1 1 1 1 2 2 1 1 1 1 1 2 2 2
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37
blong 2 2 2 2 1 2 2 2 2 1 2 1 2 2 2 2 2
38
blong 2
样本 11 12 25 30 32 被判错
> discriminiant.bayes(trainX1, trainX2, rate=21/17, testX, var.equal=T)
1 2 3 4 5 6 7 8
blong 1 1 1 1 2 2 2 2
Fisher 判别
Fisher (费歇儿)判别的准则要求类内方差尽量小、类间方差尽量大。将 Fisher 判别中的算法编写成名为 discriminiant.fisher 的函数，函数中，输入变量 TrnX1,TrnX2 分别表示 X1 类和 X2 类的训练样本，格式可以是矩阵、数据框，输入变量 TstX 是测试样本，格式是矩阵、数据框或向量，var.equal 是逻辑变量，TRUE 表示协方差阵相同，默认为不同，函数输出以 “1” 和“2” 构成的一维矩阵，“1” 代表测试样本 属于X1 类，“2” 代表属于 X2 类。
> discriminiant.fisher <- function(TrnX1, TrnX2, TstX=NULL){
+ if (is.null(TstX)==T)
+ TstX <- rbind(TrnX1, TrnX2)
+ if (is.vector(TstX)==T)
+ TstX <- t(as.matrix(TstX))
+ else if (is.matrix(TstX)!=T)
+ TstX <- as.matrix(TstX)
+ if (is.matrix(TrnX1)!=T)
+ TrnX1 <- as.matrix(TrnX1)
+ if (is.matrix(TrnX2)!=T)
+ TrnX2 <- as.matrix(TrnX2)
+
+ nx <- nrow(TstX)
+ blong <- matrix(rep(0, nx), nrow=1, byrow=T, 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
+ }
> discriminiant.fisher(trainX1, trainX2)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
blong 1 1 1 1 1 1 1 1 1 1 2 2 1 1 1 1 1 2 1 2
21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37
blong 2 2 2 2 1 2 1 2 2 1 2 1 2 2 2 2 2
38
blong 2
样本 11 12 19 25 27 30 32 被判错
> discriminiant.fisher(trainX1, trainX2, testX)
1 2 3 4 5 6 7 8
blong 1 1 1 1 2 2 2 2
三种判别法的对比
1 、距离判别法与 Fisher 判别法未对总体的分布提出特定的要求，而 Bayes 判别法要求总体的分布明确。
2 、在正态等协差阵的条件下，Bayes 判别法（不考虑先验概率的影响）等价于距离判别准则和 Fisher 线性判别法。
3 、当 K 个总体的均值向量共线性较高时，Fisher 判别法可用较少的判别函数进行判别。
4 、距离判别法和 Fisher 判别法的不足是没有考虑各总体出现的概率大小，也给不出预测的后验概率及错判率的估计，以及错判之后的损失。而这些不足恰是 Bayes 的优点。但是若给定的先验概率不符合客观实际时，Bayes 判别法也可能会导致错误的结论。
判对率
当不同类样本的协方差矩阵相同时，我们可以在 R 中使用 MASS 包的 lda 函数实现线性判别。利用 table 函数建立混淆矩阵，比对真实类别和预测类别。
> library(MASS)
> train <- bankruptcy[1:38,]
> head(train)
总负债率 收益性指标 短期支付能力 生产效率指标 类别
1 -0.45 -0.41 1.09 0.45 1
2 -0.56 -0.31 1.51 0.16 1
3 0.06 0.02 1.01 0.40 1
4 -0.07 -0.09 1.45 0.26 1
5 -0.10 -0.09 1.56 0.67 1
6 -0.14 -0.07 0.71 0.28 1
线性判对率
> distance.lda <- lda(类别~., data=train)
> table<-table(train$类别,predict(distance.lda,train)$class)
> table
1 2
1 15 2
2 3 18
> sum(diag(prop.table(table)))
[1] 0.8684211
Bayes 判对率
> bayes.lda = lda(类别~.,data = train,prior = c(17,21)/38)
> table <- table(train$类别, predict(bayes.lda, train)$class)
> table
1 2
1 15 2
2 3 18
> sum(diag(prop.table(table)))
[1] 0.8684211
线性判别分析预测股票涨跌情况
一元线性判别分析( LDA 技术)
运用 Bayes 定理进行分类，假设观测分成 K 类，K>=2，即定性响应变量 Y 可以取 k 个不同的无序值。设为 πk 为一个随机选择的观测来自第 k 类的先验概率 (prior)，记为 πk=P(Y=k)。设 f(X)=Pr(X=x|Y=k)为表示第 k 类观测的 X 的概率密度。
贝叶斯定理公式：
特征变量满足正态分布，即概率密度函数：
根据贝叶斯定理变换：
取对数转换，可见最大等价于下式最大：
用 LDA 预测股票涨跌
> library(ISLR)
> library(MASS)
> attach(Smarket)
由于 Smarket 数据集中数据较多我们采用 lda 函数，选取 Lag1、Lag2 变量样本进行分析，将 2005 年前的数据作为训练样本，2005 年的数据作为测试样本。
> lda.fit <- lda(Direction~Lag1+Lag2, data=Smarket, subset=Year<2005)
> lda.fit
Call:
lda(Direction ~ Lag1 + Lag2, data = Smarket, subset = Year <
2005)
Prior probabilities of groups:
Down Up
0.491984 0.508016
Group means:
Lag1 Lag2
Down 0.04279022 0.03389409
Up -0.03954635 -0.03132544
Coefficients of linear discriminants:
LD1
Lag1 -0.6420190
Lag2 -0.5135293
Group means 对每类每个变量计算平均，是用来估计参数 μ 。通过 Group means 矩阵可看出：当股票下跌时，前两天的投资回报率会趋向于正；当股票上涨时，前两天的投资回报率会趋向于负。Coefficients of linear discriminants 则是线性模型的系数，说明当 -0.642Log-0.513Log 很大时，LDA分类器预测上涨；-0.642Log-0.513Log 很小时，预测下跌。
> Smarket.2005 <- subset(Smarket, Year=2005)
> lda.pred <- predict(lda.fit, Smarket.2005)$class
> data.frame(lda.pred)[1:5,]
[1] Up Down Down Up Up
Levels: Down Up
> table <- table(lda.pred, Smarket.2005$Direction)
> sum(diag(prop.table(table)))
[1] 0.5248
