discriminant_analysis
YuDalang
2016-09-04
判别分析有：距离判别、贝叶斯判别、FISHER判别
1. 数据展示
head(iris)## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosa数据的图形展示
plot(iris[,1:4],col=as.numeric(iris[,5]))#同pairs()函数相同
2. 距离判别
在下面的程序中：程序分别考虑了总体协方差阵相同和总体协方差阵不同的两种情况。输入变量TrnX表示训练样本，其输入格式是矩阵（样本按行输入）或者数据框。TrnG是因子变量，便是输入样本的分类情况。输入变量TstX是待测样本，其输入格式同TrnX如果不输入那么待测样本为训练样本。输入变量var.equal是逻辑变量，TURE时表示的是总体方差相同，否则是不同。
#编写距离判别公式:来自统计建模
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)))
}
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(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,])
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,]))
}
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
}
#用测试数据进行计算
X<-iris[,1:4]
G<-gl(3,50)
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 25 26
## 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 1 1
## 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
## blong 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72
## blong 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 2
## 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95
## blong 3 2 2 2 2 2 2 2 2 2 2 3 2 2 2 2 2 2 2 2 2 2 2
## 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113
## blong 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3
## 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130
## blong 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147
## blong 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## 148 149 150
## blong 3 3 371、73和84号样本错判，回代的判别正确率为147/150=147/150
3. Baye判别
在下面的程序中：变量p是先验概率
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(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,])
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
}
#下面用Iris数据进行测试
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 25 26
## 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 1 1
## 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
## blong 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72
## blong 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 2 3 2
## 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95
## blong 3 2 2 2 2 3 2 2 2 2 2 3 2 2 2 2 2 2 2 2 2 2 2
## 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113
## blong 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3
## 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130
## blong 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147
## blong 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## 148 149 150
## blong 3 3 3第69/71/73/78和81号样本错判，正确率为145/150=145/150
4. Fisher 线性判别
- 在R中实现线性判别分析的命令是 MASS 包中的LDA函数(linear discriminant Analysis)
lda(formula , data , ... , subset , na.action)
说明：formula用法为groups~x1+x2+…… ,groups表明总体的已知的确定的分类，x1,x2,...表示分类指标，subset指明训练样本。- 下面以 Iris 数据为例进行判别分析, 首先从150条数据中抽取90条数据建立判别函数
require("MASS")## Loading required package: MASSset.seed(1314)
train <- sample(1:150, 90)
#table是从中随机选出60个作为训练集
table(iris$Species[train])##
## setosa versicolor virginica
## 30 31 29进行线性判别分析
#LDA函数(linear discriminant Analysis)是实现线性判断的一个函数
z <- lda(Species ~ ., iris, prior = c(1,1,1)/3, subset = train)#点表示所有的分类变量
#线性判别分析的结果
print(z)## Call:
## lda(Species ~ ., data = iris, prior = c(1, 1, 1)/3, subset = train)
##
## Prior probabilities of groups:
## setosa versicolor virginica
## 0.3333333 0.3333333 0.3333333
##
## Group means:
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## setosa 5.013333 3.456667 1.440000 0.2366667
## versicolor 5.967742 2.832258 4.287097 1.3451613
## virginica 6.582759 2.986207 5.582759 1.9827586
##
## Coefficients of linear discriminants:
## LD1 LD2
## Sepal.Length 0.8095703 0.2655179
## Sepal.Width 1.4207588 -2.4118265
## Petal.Length -2.0779530 0.5599844
## Petal.Width -2.9318010 -2.3536921
##
## Proportion of trace:
## LD1 LD2
## 0.9943 0.0057预测剩余60个样品的类别号
# iris[-train, ]
# dim(iris[-train, ])
# iris[-train, ]这个语句的意思就是，取train序号之外的60行数据
test.sp<-predict(z, iris[-train, ])$class
test.sp## [1] setosa setosa setosa setosa setosa setosa
## [7] setosa setosa setosa setosa setosa setosa
## [13] setosa setosa setosa setosa setosa setosa
## [19] setosa setosa versicolor versicolor versicolor versicolor
## [25] versicolor versicolor versicolor versicolor versicolor versicolor
## [31] virginica versicolor versicolor versicolor versicolor versicolor
## [37] versicolor versicolor versicolor virginica virginica virginica
## [43] virginica virginica virginica virginica virginica virginica
## [49] virginica virginica virginica versicolor virginica virginica
## [55] virginica virginica virginica virginica virginica virginica
## Levels: setosa versicolor virginicatest.true<-iris[-train,5]
table(test.true,test.sp)## test.sp
## test.true setosa versicolor virginica
## setosa 20 0 0
## versicolor 0 18 1
## virginica 0 1 20线性判别分析的投影效果
plot(z,col=as.numeric(iris[train,5]))
str(z)## List of 10
## $ prior : Named num [1:3] 0.333 0.333 0.333
## ..- attr(*, "names")= chr [1:3] "setosa" "versicolor" "virginica"
## $ counts : Named int [1:3] 30 31 29
## ..- attr(*, "names")= chr [1:3] "setosa" "versicolor" "virginica"
## $ means : num [1:3, 1:4] 5.01 5.97 6.58 3.46 2.83 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:3] "setosa" "versicolor" "virginica"
## .. ..$ : chr [1:4] "Sepal.Length" "Sepal.Width" "Petal.Length" "Petal.Width"
## $ scaling: num [1:4, 1:2] 0.81 1.421 -2.078 -2.932 0.266 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:4] "Sepal.Length" "Sepal.Width" "Petal.Length" "Petal.Width"
## .. ..$ : chr [1:2] "LD1" "LD2"
## $ lev : chr [1:3] "setosa" "versicolor" "virginica"
## $ svd : num [1:2] 36.96 2.81
## $ N : int 90
## $ call : language lda(formula = Species ~ ., data = iris, prior = c(1, 1, 1)/3, subset = train)
## $ terms :Classes 'terms', 'formula' length 3 Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width
## .. ..- attr(*, "variables")= language list(Species, Sepal.Length, Sepal.Width, Petal.Length, Petal.Width)
## .. ..- attr(*, "factors")= int [1:5, 1:4] 0 1 0 0 0 0 0 1 0 0 ...
## .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. ..$ : chr [1:5] "Species" "Sepal.Length" "Sepal.Width" "Petal.Length" ...
## .. .. .. ..$ : chr [1:4] "Sepal.Length" "Sepal.Width" "Petal.Length" "Petal.Width"
## .. ..- attr(*, "term.labels")= chr [1:4] "Sepal.Length" "Sepal.Width" "Petal.Length" "Petal.Width"
## .. ..- attr(*, "order")= int [1:4] 1 1 1 1
## .. ..- attr(*, "intercept")= int 1
## .. ..- attr(*, "response")= int 1
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. ..- attr(*, "predvars")= language list(Species, Sepal.Length, Sepal.Width, Petal.Length, Petal.Width)
## .. ..- attr(*, "dataClasses")= Named chr [1:5] "factor" "numeric" "numeric" "numeric" ...
## .. .. ..- attr(*, "names")= chr [1:5] "Species" "Sepal.Length" "Sepal.Width" "Petal.Length" ...
## $ xlevels: Named list()
## - attr(*, "class")= chr "lda"为了与前面两个进行比较，用所有的样本进行预测
z <- lda(Species ~ ., iris, prior = c(1,1,1)/3)
print(z)## Call:
## lda(Species ~ ., data = iris, prior = c(1, 1, 1)/3)
##
## Prior probabilities of groups:
## setosa versicolor virginica
## 0.3333333 0.3333333 0.3333333
##
## Group means:
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## setosa 5.006 3.428 1.462 0.246
## versicolor 5.936 2.770 4.260 1.326
## virginica 6.588 2.974 5.552 2.026
##
## Coefficients of linear discriminants:
## LD1 LD2
## Sepal.Length 0.8293776 0.02410215
## Sepal.Width 1.5344731 2.16452123
## Petal.Length -2.2012117 -0.93192121
## Petal.Width -2.8104603 2.83918785
##
## Proportion of trace:
## LD1 LD2
## 0.9912 0.0088iris.pred<-predict(z)$class #同predict(z,iris)$class
table(iris$Species,iris.pred)## iris.pred
## setosa versicolor virginica
## setosa 50 0 0
## versicolor 0 48 2
## virginica 0 1 49准确率为：147/150
注意：
Fisher判别法的重点就是选择适当的“投影轴”，对于线性判别来说，可以先将样本点投影到一维空间（u(x)=a1x1），如果划分的效果不理想，可以考虑投影到二维空间（u(x)=a1x1+a2x2），以此类推。对于非线性判别来说，可以使用二次型判别函数，可以将样本点投影到若干种二次曲面中，实现理想的判别效果。
往往实际应用中，复杂的数据无法使用线性判别函数得到理想的判别效果，就必须使用类似于二次判别函数的非线性分类方法。
关于Fisher线性判别和二次判别方法的实现可以使用MASS包中的lda()函数和qda()函数。lda()函数和qda()函数的语法和参数如下：
lda(formula, data, …, subset, na.action) lda(x, grouping, prior = proportions, tol = 1.0e-4,method, CV = FALSE, nu, …) lda(x, grouping, …, subset, na.action) qda(formula, data, …, subset, na.action) qda(x, grouping, prior = proportions, method, CV = FALSE, nu, …) qda(x, grouping, …, subset, na.action)
formula 指定参与模型计算的变量，以公式形式给出，类似于y=x1+x2+x3；
data 用于指定需要分析的数据对象；
na.action 指定缺失值的处理方法，默认情况下，缺失值的存在使算法无法运行，当设置为“na.omit”时则会删除含有缺失值的样本；
x 指定需要处理的数据，可以是数据框形式，也可以是矩阵形式；
grouping 为每个观测样本指定所属类别；
prior 可为各个类别指定先验概率，默认情况下用各个类别的样本比例作为先验概率；
method 指定计算均值和方差的方法，默认情况下使用最大似然估计，还可以使用mve法和稳健估计法（基于t分布）；
CV 是否进行交叉验证，默认情况下不进行交叉验证。例子
# Fisher判别法的应用
#通过抽样建立训练样本和测试样本
index <- sample(2,size = nrow(iris),
replace = TRUE, prob = c(0.75,0.25))
train <- iris[index == 1,]
test <- iris[index == 2,]
#加载R包并使用Fisher线性判别法
library(MASS)
lda_res3 <- lda(Species~., data = train)
lda_pre3 <- predict(lda_res3, newdata = test[,1:4])
#生成实际与预判交叉表和预判精度
table(test$Species,lda_pre3$class)##
## setosa versicolor virginica
## setosa 16 0 0
## versicolor 0 13 0
## virginica 0 2 16sum(diag(table(test$Species,lda_pre3$class)))/sum(table(test$Species,lda_pre3$class))## [1] 0.9574468#使用Fisher二次判别法
lda_res4 <- qda(Species~., data = train)
lda_pre4 <- predict(lda_res3, newdata = test[,1:4])
#生成实际与预判交叉表和预判精度
table(test$Species,lda_pre4$class)##
## setosa versicolor virginica
## setosa 16 0 0
## versicolor 0 13 0
## virginica 0 2 16sum(diag(table(test$Species,lda_pre4$class)))/sum(table(test$Species,lda_pre4$class))## [1] 0.9574468