吴世涵的博客
0%

【数学建模笔记 16】数学建模的多元分析

发表于 分类于
本文字数： 11k 阅读时长 ≈ 10 分钟
多元分析是多变量的统计分析方法，是数理统计中应用广泛的一个重要分支。主要包括判别分析、主成分分析、因子分析、聚类分析等。本文主要介绍了判别分析、主成分分析、因子分析、聚类分析等内容。

29. 多元分析

定义

多元分析是多变量的统计分析方法，是数理统计中应用广泛的一个重要分支。

判别分析

判别分析是一种分类方法。假定有 类判别对象 ，每一类 由 个指标的 个样本确定，即 其中 矩阵的第 行是 的第 个样本点的观测值向量。

记 ， 分别表示 ​ 的均值向量和离差矩阵，即

对待判定对象 ​，有一个一般规则，可以依据 ​ 的值，对 ​ 属于 的哪一类作出判别，称判别规则，其函数称判别函数，记 。

距离判别法

根据距离最近原则判别。 若 则 。

距离 一般用马氏距离， 个总体协方差矩阵相等时

不等时

Fisher 判别法

Fisher 判别法的思想是，将样例投影到一条或多条直线上，使同类样例的投影点尽可能接近，不同类样例的投影点尽可能远离。

对于二分类情况，若给定直线 ​​​，则第 ​ 类样本的中心在直线上的投影为 ​。

第 ​ 类中第 ​ 个样本的组内偏差为 ​，故第 ​ 类的组内偏差之和为

去掉 得 即离差矩阵。

所有类的组内偏差总和即为 要使同类样例的投影点尽可能接近，只需使 ​​ 最小化。

而组间偏差同理有 要使不同类样例的投影点尽可能远离，只需使 ​​ 最大化。

于是构造代价函数 使得 ​ 最大，也就等价于

根据拉格朗日乘子法，即 则 得 ，即 为矩阵 ​ 的特征向量。于是

于是有直线 ，和 ​​ 垂直的两类的中心线可作为两类的判别直线。

即有阈值 ​​，比较 ​ 与 ​ 的大小，得出分类。

对于多分类情况，则有

之后求解 的闭式解为 ​ 的 个最大广义特征值所对应的特征向量组成的矩阵。

主成分分析

主成分分析，是利用降维把多指标转化为几个综合指标的多元统计分析方法。

设 表示以 ​ 为样本观测值的随机变量，如果能找到 ，使得方差 最大，就表明这 个变量的最大差异。

一般说来，代表原来 ​ 个变量的主成分不止一个，但不同主成分的信息不能相互包含，即协方差为 0，在几何上即方向正交。

设 表示第 个主成分 其中 ，

使 ​ 最大，​

使 并使 ​​​​​ 最大，

​ 使 并使 ​ 最大，

……

为实现上述目标，步骤如下：

  1. 假设有 ​ 个对象，​ 个指标 ，设第 个对象的第 个指标为 ，对原来的 ​​​ 个指标进行标准化 其中 得标准化的数据矩阵 ​​，记标准化后的指标为 ；

  2. 求相关系数矩阵 ，

  3. 计算相关系数矩阵 的特征值 及对应的标准正交化特征向量 ，其中 ，组成 个新的指标变量

  4. 计算主成分贡献率 累计贡献率

因子分析

因子分析将原始变量分解为若干个因子的线性组合 是 的期望变量， 称公共因子向量， 称特殊因子向量， 称因子载荷矩阵， 是变量 在公共因子 上的载荷，反映 对 的重要程度。

通常假设 ​ 互不相关且具有单位方差，​ 互不相关且与 ​​​ 互不相关， 为对角阵，于是有 每个原始变量 的方差都可以分解成共性方差 和特殊方差 之和，其中 反映全部公共因子对 的方差贡献，​ 反映特殊因子的方差贡献。

令 ，则 是公共因子 对 总方差的贡献，称 为 的贡献率。

利用主成分分析法求因子载荷矩阵，设 为相关系数矩阵 的特征值， 为对应的正交特征向量，，则 特殊因子的方差用 的对角元估计，即 一般来说，理想的载荷结构是，每一列或每一行各载荷平方接近 0 或 1，需要对因子载荷矩阵进行旋转。

选取方差最大的正交旋转矩阵 ​​，使得因子上的载荷尽量拉开距离，从而使其接近 0 或 1

聚类分析

聚类分析，就是指将相似元素聚为一类。

层次聚类

基本思想：距离相近的样品先聚为一类，距离远的后聚成类，步骤如：

  1. 将每个样品独自聚成一类，构造 个类；
  2. 根据选定的距离公式，计算样品两两间距离，构造距离矩阵；
  3. 把距离最近的两类归为一类，聚成 类；
  4. 计算新类与当前各类的距离，再选距离最近的两类归为一类，聚成 类；
  5. 按步骤 4 重复迭代，直到所有样品聚为一类。

类和类之间的距离有各种定义，如最短距离法，对于类 ​，有 除此之外，还有最长距离法、中间距离法等。

K 均值聚类

动态聚类法的思想是先粗略分一下类，然后按某种最优原则进行修正，直到类分得较合理为止，K 均值法是动态聚类的一种方法。

假定样本可分为 ，并选定 ​ 个初始聚类中心，然后根据最小距离原则将每个样本分配到某一类中，之后计算新的聚类中心，并调整聚类情况，直到收敛。

步骤如下：

  1. 将样本集任意划分为 类，记 ，计算聚类中心 ， 其中 为当前 类中的样本数，以及计算误差平方和

  2. 再令 ​，按最小距离原则重新聚类，并重新计算聚类中心；

  3. 若连续两次迭代 不变，算法终止，否则转步骤 2。

最佳簇数 k 的确定

方法有簇内离差平方和拐点法和轮廓系数法。

  • 簇内离差平方和拐点法

计算不同 值下计算簇内离差平方和，然后画图找到拐点。当斜率由大突然变小，且之后斜率变化缓慢，则认为突然变换的点就是寻找的目标点。 再增加聚类效果不再有大的变化。

  • 轮廓系数法

思想为使簇内样本密集，簇间样本分散。

对于每个样本点 ​ 有：

  • 簇内不相似度 ​：样本点 ​ 与其所在簇内其他样本的平均距离；
  • 簇间不相似度 ：样本点 与其他簇样本的平均距离的最小值。

于是有样本点 ​ 的轮廓系数 总轮廓系数为所有样本点轮廓系数的平均值。

总轮廓系数小于 0 说明聚类效果不佳；接近 1 说明簇内样本平均距离非常小，簇间最近距离非常大，聚类效果非常理想。

Python 代码

距离判别法

蠓虫分为 Af 和 Apf，测得 Af、Apf 的触角长度和翅膀长度数据。

Af: ​，​，，，​；

Apf: ，，，，；

试判别 属于哪一类，代码如下：

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
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
#! /usr/bin/env python
# -*- coding: utf-8 -*-
# @ author: Koorye
# @ date: 2021-7-30
# @ function: 距离判别法

# %%

import numpy as np
from sklearn.neighbors import KNeighborsClassifier

# %%

# 源数据
x0 = np.array([[1.24, 1.27],
              [1.36, 1.74],
              [1.38, 1.64],
              [1.38, 1.82],
              [1.38, 1.90],
              [1.14, 1.78],
              [1.18, 1.96],
              [1.20, 1.86],
              [1.26, 2.00],
              [1.28, 2.00]])

# 前 5 个属于 1 类 af，后 5 个属于 2 类 apf
label = np.array([1 for i in range(5)] + [2 for i in range(5)])

# 待判别数据
x = np.array([[1.24, 1.80]])

# %%

# 画图
import matplotlib.pyplot as plt

plt.scatter(x0[:5,0], x0[:5,1],label='Af')
plt.scatter(x0[5:,0], x0[5:,1],label='Apf')
plt.scatter(x[0][0],x[0][1], label='Unknown')
plt.legend()

# %%

# 协方差矩阵
v = np.cov(x0.T)

# 模型
knn = KNeighborsClassifier(2,
                           metric='mahalanobis',
                           metric_params={'V': v})
knn.fit(x0, label)
knn.predict(x)

输出如下：

1
array([2])

从图中大致可以看出未知样本接近 Apf 类，而计算结果为 2，即 Apf 类，与直觉相吻合。

Fisher 判别法

对上例使用 Fisher 判别法求解，代码如下：

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
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
#! /usr/bin/env python
# -*- coding: utf-8 -*-
# @ author: Koorye
# @ date: 2021-7-30
# @ function: Fisher 判别法

# %%

import numpy as np
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis

# %%

# 源数据
x0 = np.array([[1.24, 1.27],
              [1.36, 1.74],
              [1.38, 1.64],
              [1.38, 1.82],
              [1.38, 1.90],
              [1.14, 1.78],
              [1.18, 1.96],
              [1.20, 1.86],
              [1.26, 2.00],
              [1.28, 2.00]])

# 前 5 个属于 1 类 af，后 5 个属于 2 类 apf
label = np.array([1 for i in range(5)] + [2 for i in range(5)])

# 待判别数据
x = np.array([[1.24, 1.80]])

# %%

# 协方差矩阵
v = np.cov(x0.T)

# 模型

model = LinearDiscriminantAnalysis()
model.fit(x0, label)
model.predict(x)

输出如下：

1
array([2])

结果与距离判别法一致。

主成分分析

对下列数据进行主成分分析

序号 身高 x1 胸围 x2 体重 x3
1 149.5 69.5 38.5
2 162.5 77.0 55.5
3 162.7 78.5 50.8
4 162.2 87.5 65.5
5 156.5 74.5 49.0

代码如下：

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
27
28
29
30
31
32
33
34
35
#! /usr/bin/env python
# -*- coding: utf-8 -*-
# @ author: Koorye
# @ date: 2021-7-30
# @ function: 主成分分析

# %%

import numpy as np
import pandas as pd
from sklearn.decomposition import PCA

# %%

# 源数据
df = pd.DataFrame({
    'x1': [149.5, 162.5, 162.7, 162.2, 156.5],
    'x2': [69.5, 77, 78.5, 87.5, 74.5],
    'x3': [38.5, 55.5, 50.8, 65.5, 49]
})

# 模型
model = PCA().fit(np.array(df))

# %%

print('特征值:', model.explained_variance_)
print('贡献率:', model.explained_variance_ratio_)
print('各主成分的系数:', model.components_)

# %%

pca_df = pd.DataFrame(model.transform(np.array(df)))
pca_df.columns = ['F1', 'F2', 'F3']
pca_df

输出如下：

1
2
3
4
5
特征值: [161.47423448   9.60080441   2.28996111]
贡献率: [0.93141196 0.05537914 0.0132089 ]
各主成分的系数: [[-0.39412803 -0.5037512  -0.76869878]
 [-0.90779807  0.34389304  0.24008381]
 [ 0.14340765  0.79244703 -0.59284226]]
F1 F2 F3
0 17.867546 2.409312 0.343559
1 -4.102132 -2.731441 -1.927107
2 -1.323700 -3.525555 2.076603
3 -16.960269 3.552614 0.422142
4 4.518556 0.295070 -0.915196

层次聚类

对下列 7 种岩石聚类

1 2 3 4 5 6 7
Cu 2.9909 3.2044 2.8392 2.5315 2.5897 2.9600 3.1184
W 0.3111 0.5348 0.5696 0.4528 0.3010 3.0480 2.8395
Mo 0.5324 0.7718 0.7614 0.4893 0.2735 1.4997 1.9350

代码如下：

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
27
28
29
30
#! /usr/bin/env python
# -*- coding: utf-8 -*-
# @ author: Koorye
# @ date: 2021-7-30
# @ function: 层次聚类

# %%

import numpy as np
import pandas as pd
import scipy.cluster.hierarchy as sch

# %%

# 源数据
df = pd.DataFrame({
	'Cu': [2.9909,3.2044,2.8392,2.5315,2.5897,2.9600,3.1184],
	'W':[.3111,.5348,.5696,.4528,.3010,3.0480,2.8395],
	'Mo': [.5324,.7718,.7614,.4893,.2735,1.4997,1.9350],
})

# 计算两两距离
dist = sch.distance.pdist(df)

# 转化为距离矩阵
dist_mat = sch.distance.squareform(dist)

# 聚类并画图
z = sch.linkage(dist)
sch.dendrogram(z)

输出如下：

K 均值聚类

对鸢尾花数据集进行 K 均值聚类

ID sepal length (cm) sepal width (cm)
0 5.1 3.5
1 4.9 3.0
2 4.7 3.2
3 4.6 3.1
4 5.0 3.6
... ... ...
145 6.7 3.0
146 6.3 2.5
147 6.5 3.0
148 6.2 3.4
149 5.9 3.0

代码如下：

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
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
#! /usr/bin/env python
# -*- coding: utf-8 -*-
# @ author: Koorye
# @ date: 2021-7-30
# @ function: K 均值聚类

# %%

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score

# %%

# 源数据
df = pd.DataFrame(load_iris()['data'], columns=load_iris()['feature_names'])

# 计算 k 取不同值时的轮廓系数
score_list = []
for i in range(2,10):
    model = KMeans(i)
    model.fit(df.iloc[:,:2])
    score_list.append(silhouette_score(df, model.labels_))

plt.plot([i for i in range(2,10)], score_list)

# %%

model = KMeans(3)
model.fit(df.iloc[:,:2])
df2 = df.iloc[:,:2].copy()
df2['label'] = model.labels_

from plotnine import *

(
    ggplot(df2,aes('sepal length (cm)', 'sepal width (cm)', color='label'))
    + geom_point()
    + theme_matplotlib()
)

输出如下：

可以看到，分 3 类是轮廓系数最大，聚类效果如图。