判别分析 练习题

1. 考虑两个数据集
  \(x_1=\begin{bmatrix}3 & 7\\ 2 & 4\\ 4 & 7\end{bmatrix}\) ,\(x_2=\begin{bmatrix}6 & 9\\ 5 & 7\\ 4 & 8\end{bmatrix}\),其中\(\bar{x_1}=\begin{bmatrix}3\\ 6\end{bmatrix}\),\(\bar{x_2}=\begin{bmatrix}5\\ 8\end{bmatrix}\),\(S_p=\begin{bmatrix}1 & 1\\ 1 & 2\end{bmatrix}\)。

1.1

  计算Fisher线性判别函数。    
  首先将数据集\(x_1\)、\(x_2\)、\(\bar x_1\)、\(\bar x_2\)、和合并协方差阵\(S_p\)输入:
1
2
3
4
5
6
x1 = c(3,2,4,6,5,4)
x2 = c(7,4,7,9,7,8)
y = c(1,1,1,2,2,2)
(data1 = data.frame(x1,x2,y))
plot(data1[,-3])
text(x1,x2,y,adj=-0.5)
1
2
3
library(MASS)

ld=lda(y~x1+x2,data = data1);ld
判别函数为: y=x1+(5.5511e-17)*x2
1
2
3
4
5
x1_bar = data1[1:3,-3];x1_bar  # 1-3行, 去掉 第3标签列: 第一组数据
x2_bar = data1[4:6,-3];x2_bar # 4-6行, 去掉 第3标签列: 第二组数据

(x1_mean = apply(x1_bar, 2, mean)) # 2: 按列求均值 # 第一组 按列求均值
(x2_mean = apply(x2_bar, 2, mean)) # 第二组 按列求均值
1
2
3
Sp = matrix(c(1,1,1,2),ncol = 2)
sp_re = solve(Sp) # solve求逆矩阵
(2*cov(x1_bar)+cov(x2_bar))/4
1
a = t(x1_mean-x2_mean)%*%sp_re;a # 用公式求判别直线
1
fy <- lda(y~., data = data1);fy
1
2
3
4
5
(tab = table(predict(fy)$class, data1$y)) # 每次运行都不一样,是因为 有 距离相等 情况,此时就会随机分到一个类
# 平行线,y = kx + b ,k相同->平行,点到平行线上投影相同


# Fisher 判别曲线是 x1 轴,有两个不同类的点的投影是重合的 --导致--> 判别这两个点的类别是随机的
1
sum(diag(tab))/sum(tab) # 判对率
1
predict(fy,data.frame(x1=4,x2=6))$class

1.2

  应用Bayes法则,在相同先验概率和相同代价下将观测值
  \(x_0'=(2,7)\)分类到总体\(G_1\)或\(G_2\)。
  将数据汇总成数据框的形式:
block6.1.2,comment
1
ld=lda(y~x1+x2,prior=c(1,1)/2,data = data1);ld
  应用Bayes准则并分类:
block6.1.3,comment
1
predict(ld,data.frame(x1=2,x2=7))$class
1
predict(ld,data.frame(x1=2,x2=7))$posterior
1
predict(ld,data.frame(x1=2,x2=7))$x
2.   设\(n_1=11\)个和\(n_2=12\)个观测值分别取自两个随机变量\(X_1\)和\(X_2\)。 假定这两个变量服从二元正态分布,且有相同的协方差矩阵:
\[\bar{x_1}=\begin{bmatrix}-1\\ -1\end{bmatrix},\bar{x_2}=\begin{bmatrix}2\\ 1\end{bmatrix},S_p=\begin{bmatrix}7.3 & -1.1\\ -1.1 & 4.8\end{bmatrix}\]   首先将数据集\(\bar x_1\)、\(\bar x_2\)、和合并协方差阵\(S_p\)输入
read6.2,comment
1
2
3
x1_bar = rep(-1,2)
x2_bar = c(2,1)
S_p = matrix(c(7.3,-1.1,-1.1,4.8),ncol=2)

2.1

  构造样本的Fisher线性判别函数。
  设判别函数为\(Y=a_2'X=(\bar x_1-\bar x_2)'S_p^{-1}X\)
block6.2.1,comment
1
ldarray=t(x1_bar-x2_bar)%*%S_p;ldarray
因此,判别函数为: \(y = -19.7*x1 -6.3*x2\)

2.2

  将观测值
  \(x_0'=(0,1)\)分配到总体\(G_1\)或\(G_2\)(假定有等代价和等先验概率)。
  可以直接判别界点然后将\(x_0'\)带入判别:
block6.2.2,comment
1
y1=(x1_bar-x2_bar)%*%x1_bar;y1
1
y2=(x1_bar-x2_bar)%*%x2_bar;y2
1
y0=(y1+y2)/2;y0
1
y=ldarray%*%array(c(0,1),dim = c(2,1));y
y<y0 那么判为G2 3.   以舒张期血压和血浆胆固醇含量预测被检验者是否患冠心病。测得15名冠心病患者和16名健康人的舒张压
  \(X_1\)(mmHg)及血浆胆固醇含量\(X_2\)(mg/dL),结果见下表。   首先读取数据:
read6.3,comment
1
2
3
4
library(openxlsx)

d6.3 = read.xlsx('../Res/mvexer5.xlsx', "E6.3");d6.3
d6.3 = d6.3[-1,];d6.3 # 去掉索引为1行: 1 <U+7F16><U+53F7>(Ak) X1A X2A <U+7F16><U+53F7>(Bk) X1B X2B
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
# 切割出来
d3.A = d6.3[,2:3]
d3.B = d6.3[,5:6]

# 重定义列名
colnames(d3.A) <- c("x1", "x2")
colnames(d3.B) <- c("x1", "x2")

# 对应增加 label 列
d3.A$label = 1
d3.B$label = 2

# 按行 合并两类
d3 = rbind(d3.A, d3.B);head(d3);tail(d3)
#attach(d3)

3.1

  对每一组数据用不同的符号作两变量的散点图,观察它们在平面上的散布情况,并判断对该组数据做判别分析是否合适。
plot6.3,comment
1
plot(d3[,-3],col=d3$label) # 去掉标签列,绘散点图
红点与黑点分割较为明显,适合做判别分析

3.2

  分别建立距离判别(等方差阵和不等方差阵)、Fisher判别和Bayes判别分析模型,计算各自的判别符合率,以此确定哪种判别方法最恰当。
  假设两类总体都服从正态分布。
(a)距离判别(等方差阵),与Fisher判别、先验概率为频率的Bayes判别等价
block6.3.1,comment
1
2
3
4
5
d3 <- apply(d3, 2, function(x)as.numeric(as.character(x))) # 转换为 数值型

d3 <- as.data.frame(d3) # 转为 data.frame

ld=lda(label~x1+x2, data=d3);ld
1
2
3
4
lp = predict(ld) # 在原训练集上 预测
G1 = lp$class # 判别结果

(tab1 = table(d3$label, G1)) # 判别矩阵
1
sum(diag(prop.table(tab1))) # 判对率
(b)距离判别(不等方差阵)
block6.3.2,comment
1
qd=qda(label~x1+x2, data=d3);qd
1
2
3
qp=predict(qd)
G2=qp$class
tab2=table(d3$label,G2);tab2 # 判别矩阵
1
sum(diag(tab2))/sum(tab2) # 判对率
(c)先验概率相同的Bayes判别
block6.3.3,comment
1
ld3=lda(label~x1+x2,data=d3,prior=c(1,1)/2);ld3
1
2
G3=predict(ld3)$class
tab3=table(d3$label,G3);tab3
1
sum(diag(tab2))/sum(tab3) # 判对率

3.3

  绘制线性判别函数图
block6.3.4,comment
1
2
3
4
5
6
attach(d3)
a1=ld$scaling[1,];a2=ld$scaling[2,]
a=1/a2*(a1*mean(x1)+a2*mean(x2))
b=-a1/a2
plot(ld,type='both')
abline(a,b)
4   对于A股市场2009年陷入财务困境的上市公司(ST公司),我们收集了8间ST公司陷入财务困境前的一年(2008年)的财务数据,同时对于财务良好的公司(非ST公司),收集了同一时期8间非ST公司对应的财务数据。数据涉及四个变量:资产负债率
  \(x_1\)、流动资产周转率\(x_2\)、总资产报酬率\(x_3\)和营业收入增长率\(x_4\)。类别变量G中2代表ST公司,1代表非ST公司。
read6.4,comment
1
2
3
4
detach(d3)
library(openxlsx)
d6.4 = read.xlsx("../Res/mvexer5.xlsx", 'E6.4');d6.4
attach(d6.4)
1
dim(d6.4)
1
2
duplicated(d6.4) # 去重  TRUE 的位置 就说明重复了
d4_uni <- unique(d6.4) # 去除相同的行
1
dim(d4_uni);head(d4_uni)

4.1

  分别建立线性判别函数、非线性判别函数和Bayes判别分析模型,计算各自的判别符合率,确定哪种判别方法最恰当。
(a)线性判别
block7.4.1,comment
1
2
3
4
library(MASS)
ld = lda(G ~ x1 + x2 + x3 + x4, data=d6.4);ld # 建立 线性判别模型
lp = predict(ld) # 根据线性判别模型 预测所属类别
G1 = lp$class # 预测的所属类别 结果
1
disTable = table(G, G1);disTable # 判别矩阵
1
rightRate = sum(diag(disTable))/ sum(disTable);rightRate # 判对率
(b)非线性判别
block6.4.2,comment
1
2
3
qd = qda(G ~ x1 + x2 + x3 + x4);qd # 建立 非线性判别模型 (二次判别)
qp = predict(qd) # 根据 非线性判别模型 预测所属类别
G2 = qp$class # 预测的所属类别 结果
1
disTable2 = table(G, G2);disTable2 # 判别矩阵
1
rightRate2 = sum(diag(disTable2))/ sum(disTable2);rightRate2 # 判对率
(c)Bayes判别
  假设两类总体均服从正态分布。则先验概率为频率的Bayes与线性判别的判别结果相同,下面建立先验概率相同的Bayes判别分析模型。
block6.4.3,comment
1
2
3
4
# 注意:先验概率个数 是分类 个数,这里分2类,所有2
ld2 = lda(G ~ x1 + x2 + x3 + x4, prior = c(1, 1)/2,data=d6.4);ld2 # 建立 Bayes判别模型 (先验概率相等)
lp2 = predict(ld2) # 根据 Bayes判别模型 预测所属类别
G3 = lp2$class # 预测的所属类别 结果
1
disTable3 = table(G, G3);disTable3 # 判别矩阵
1
rightRate3 = sum(diag(disTable3))/ sum(disTable3);rightRate3 # 判对率

4.2

  某公司2008年财务数据为:
  \(x_1=78.356,x_2=0.8895,x_3=1.800,x_4=14.1022\)。试判定2009年该公司是否会陷入财务困境。
block6.4.4,comment
1
predict(ld, data.frame(x1=78.356,x2=0.8895,x3=1.800,x4=14.1022))
判为2(陷入危机)的概率为87%,可信,2009年该公司会陷入财务困境。
1
2
# 使用完毕,解绑
detach(d6.4)
5.   植物分类之判别分析:费歇(Fisher)于1936年发表的鸢尾花(Iris)数据被广泛的作为判别分析的经典例子。数据是对3种鸢尾花(g):刚毛鸢尾花(第一组)、变色鸢尾花(第二组)和弗吉尼亚鸢尾花(第三组)各抽取一个容量为50的样本,测量其花萼长(sepallen,\(x_1\))、花萼宽(sepalwid,\(x_2\))、花瓣长(petallen,\(x_3\))、花瓣宽(petalwid,\(x_4\)),单位为mm,数据见下表。   首先读取数据:
read6.5,comment
1
2
3
4
5
6
library(openxlsx)
d1 = read.xlsx("../Res/mvexer5.xlsx", 'E6.5', cols=c(2,3,4,5));d1$label = 1;d1 # 第一组
d2 = read.xlsx("../Res/mvexer5.xlsx", 'E6.5', cols=c(7,8,9,10));d2$label = 2;d2 # 第二组
d3 = read.xlsx("../Res/mvexer5.xlsx", 'E6.5', cols=c(12,13,14,15));d3$label = 3;d3 # 第三组

d6.5 = rbind(d1,d2,d3);d6.5 # 行合并
1
2
3
library(openxlsx)
d5_csv <- read.csv('../Res/mvexer5-E6.5.csv',header=T)
head(d5_csv);dim(d5_csv)
1
2
3
4
5
6
7
8
9
10
11
12
13
d5iris <- d5_csv[-1,]
G1 <- d5iris[,2:5];G1
G2 <- d5iris[,7:10];G2
G3 <- d5iris[,12:15];G3

colnames(G1) <- c("x1", "x2", "x3", "x4")
colnames(G2) <- c("x1", "x2", "x3", "x4")
colnames(G3) <- c("x1", "x2", "x3", "x4")

diris = rbind(G1, G2, G3)
head(diris)
diris <- apply(diris, 2, function(x)as.numeric(as.character(x))) # 转换为 数值型
head(diris)
1
2
G <- rep(1:3, c(50, 50, 50))
dim(diris)
1
2
3
4
Giris <- cbind(diris, G);tail(Giris)
Giris <- as.data.frame(Giris)
library(MASS)
firis <- lda(G~x1+x2+x3+x4, data=Giris)

5.1

  非线性判别
block6.5.1,comment
1

5.2

  线性判别与Bayes判别
block7.5.2,comment
1
2
ld = lda(label~., data = d6.5);ld # 线性判别
plot(ld,col=d6.5$label)
1
2
3
lp = predict(ld) # 根据线性判别模型 预测所属类别
lpC = lp$class # 预测的所属类别结果
tab = table(d6.5$label,lpC);tab # 判别矩阵
1
sum(diag(tab))/sum(tab) # 判对率
Q&A 补充 参考