16.1
对以下样本数据进行主成分分析:
X
=
[
2
3
3
4
5
7
2
4
5
5
6
8
]
X = \left[\begin{array}{llllll}2 & 3 & 3 & 4 & 5 & 7 \\ 2 & 4 & 5 & 5 & 6 & 8\end{array}\right]
X=[223435455678]
由于手解数据不是那么“友好”所以直接用代码求解:
import numpy as np
X = np.array([[2, 3, 3, 4, 5, 7],
[2, 4, 5, 5, 6, 8]], dtype='float')
def normalize_data(data_array):
m, n = data_array.shape
for i in range(m):
data_array[i] = data_array[i] - data_array[i].mean()
data_array[i] = data_array[i] / np.sqrt(data_array[i].var())
return data_array
X = normalize_data(X)
# 利用奇异值分解进行PCA
X_prime = X.T / np.sqrt(X.shape[1] - 1)
U, Sigma, VT = np.linalg.svd(X_prime)
print(f'主成分方差贡献相对大小/主成分方差为:{Sigma}')
print(f'主成分所在的轴向/主成分投影矢量为:{[x for x in VT]}')
# 计算主成分与各原变量的相关系数/因子负荷
factor_loading = np.zeros((VT.shape[0], VT.shape[1]))
for j in range(factor_loading.shape[0]):
for i in range(factor_loading.shape[1]):
factor_loading[j, i] = Sigma[j] * VT[j, i]
print(f'主成分与各原变量的相关系数/因子负荷:\n{factor_loading}')
# 计算主成分矢量对于各样本的方差贡献率
contribution_2samples = np.zeros(VT.shape[1])
for i in range(len(contribution_2samples)):
contribution_2samples[i] = (factor_loading[:, i] ** 2).sum()
print(f'主成分矢量对于各样本的方差贡献率:\n{contribution_2samples}')
# 计算各个样本的主成分值
principle_conmponents = np.zeros((VT.shape[0], X.shape[1]))
for i in range(X.shape[1]):
principle_conmponents[:, i] = VT @ X[:, i]
print(f'PCA matrix:\n{principle_conmponents}')
主成分方差贡献相对大小/主成分方差为:[1.52983485 0.24414203]
主成分所在的轴向/主成分投影矢量为:[array([0.70710678, 0.70710678]), array([ 0.70710678, -0.70710678])]
主成分与各原变量的相关系数/因子负荷:
[[ 1.0817566 1.0817566 ]
[ 0.17263449 -0.17263449]]
主成分矢量对于各样本的方差贡献率:
[1.2 1.2]
PCA matrix:
[[-2.02792041 -0.82031104 -0.4330127 0. 0.82031104 2.46093311]
[ 0.2958696 -0.04571437 -0.4330127 0. 0.04571437 0.1371431 ]]
分析:
首先从主成分方差可以看出,第一个主成分远大于第二个,所以数据主要分布在第一个轴上(贡献率
1.52983485
0.24414203
+
1.52983485
=
86.2
%
\frac{1.52983485}{0.24414203 + 1.52983485}= 86.2\%
0.24414203+1.529834851.52983485=86.2%),或者说其实数据本身更接近一个一维分布;
从主成分投影矢量可以看出,第一个轴其实就是二维坐标系中,与原来的 x y xy xy轴呈45度的, x = y x=y x=y直线的方向,第二个就是正交的,在原来 x y xy xy坐标系里 x = − y x=-y x=−y的方向。这个特点其实从原本的数据中是能看出来的,数据确实主要分布在 x = y x=y x=y直线的附近;从SVD分析的角度,可知,多个样本主要贡献了/对应了两种模式,对于最重要的第一种模式,其对原本的两个特征贡献相同,说明这种模式贡献到了原两个特征的平分处(角平分线),也反映了数据主要方差集中在这个方向,对于第二种模式,贡献也是相同的,正负反映了其对应的贡献方向(原基组下为 x = − y x=-y x=−y)
从相关系数可以看出,第一个主成分(投影到 x = y x=y x=y直线上的分量)贡献最大,且确实对原本的变量贡献程度相等,第二个主成分是投影到 x = − y x=-y x=−y直线上的数据值,相关程度也相等,正负相关,说明沿着这个主成分轴方向x值(原第一变量)变大,y值(原第二变量)变小,确实也如此;
从主成分矢量对于样本方差贡献率可知,因为没有截断,所以是完全贡献,超过1是因为这里是样本PCA,不是总体/布居PCA;
从投影主成分轴之后的PCA matrix
可以看出,其基本都落在新的第一主成分轴上(因为第二主成分都接近0),也就是原变量标准基下的
x
=
y
x=y
x=y轴。
16.2
证明样本协方差矩阵 S S S是总体协方差矩阵 Σ \Sigma Σ的无偏估计。
证明:
若
x
1
,
x
2
,
⋯
,
x
n
x_{1}, x_{2}, \cdots, x_{n}
x1,x2,⋯,xn来自独立同分布
X
X
X,其分布满足:
E
(
X
)
=
μ
E(X)=\mu
E(X)=μ;
Cov
(
x
)
=
E
[
(
X
−
μ
)
(
X
−
μ
)
T
]
=
Σ
\operatorname{Cov}(x)=E\left[(X-\mu)(X-\mu)^{T}\right]=\Sigma
Cov(x)=E[(X−μ)(X−μ)T]=Σ。n个x为取出的样本,则样本均值为:
x
ˉ
=
1
n
(
x
1
+
x
2
+
⋯
+
x
n
)
\bar{x}=\frac{1}{n}\left(x_{1}+x_{2}+\cdots+x_{n}\right)
xˉ=n1(x1+x2+⋯+xn)
则样本均值的期望为:
E
(
x
ˉ
)
=
E
(
x
1
+
x
2
+
⋯
+
x
n
n
)
=
1
n
(
E
(
x
1
)
+
E
(
x
2
)
+
⋯
+
E
(
x
n
)
)
=
1
n
(
n
×
μ
)
=
μ
\begin{aligned} E(\bar{x}) &=E\left(\frac{x_{1}+x_{2}+\cdots+x_{n}}{n}\right) \\ &=\frac{1}{n}\left(E\left(x_{1}\right)+E\left(x_{2}\right)+\cdots+E\left(x_{n}\right)\right) \\ &=\frac{1}{n}(n \times \mu) =\mu \end{aligned}
E(xˉ)=E(nx1+x2+⋯+xn)=n1(E(x1)+E(x2)+⋯+E(xn))=n1(n×μ)=μ
其中第二个等号是因为它们都来自相同分布,因此期望都相同。
样本均值的协方差为:
Cov
(
x
ˉ
)
=
Cov
(
x
ˉ
,
x
ˉ
)
=
Cov
(
∑
i
=
1
n
1
n
x
i
,
∑
j
=
1
n
1
n
x
j
)
=
∑
i
=
1
n
∑
j
=
1
n
Cov
(
1
n
x
i
,
1
n
x
j
)
\begin{aligned} \operatorname{Cov}(\bar{x}) &=\operatorname{Cov}(\bar{x}, \bar{x}) \\ &=\operatorname{Cov}\left(\sum_{i=1}^{n} \frac{1}{n} x_{i}, \sum_{j=1}^{n} \frac{1}{n} x_{j}\right) \\ &=\sum_{i=1}^{n} \sum_{j=1}^{n} \operatorname{Cov}\left(\frac{1}{n} x_{i}, \frac{1}{n} x_{j}\right) \end{aligned}
Cov(xˉ)=Cov(xˉ,xˉ)=Cov(i=1∑nn1xi,j=1∑nn1xj)=i=1∑nj=1∑nCov(n1xi,n1xj)
最后一个等于号是因为各个样本都是独立同分布(iid),所以互不相关,所以求和号可以提出来。且因为互不相关,因此对于
i
≠
j
i \neq j
i=j:
Cov
(
1
n
x
i
,
1
n
x
j
)
=
1
n
2
Cov
(
x
i
,
x
j
)
=
0
\operatorname{Cov}\left(\frac{1}{n} x_{i}, \frac{1}{n} x_{j}\right)=\frac{1}{n^{2}} \operatorname{Cov}\left(x_{i}, x_{j}\right)=0
Cov(n1xi,n1xj)=n21Cov(xi,xj)=0
从而可以去掉一个求和号:
Cov
(
x
ˉ
)
=
∑
i
=
1
n
∑
j
=
1
n
Cov
(
1
n
x
i
,
1
n
x
j
)
=
∑
i
=
1
n
Cov
(
1
n
x
i
,
1
n
x
i
)
=
∑
i
=
1
n
1
n
2
Cov
(
x
i
,
x
i
)
=
∑
i
=
1
n
1
n
2
Cov
(
x
i
)
=
∑
i
=
1
n
1
n
2
Cov
(
X
)
=
n
×
1
n
2
Σ
=
1
n
Σ
\begin{aligned} \operatorname{Cov}(\bar{x}) &=\sum_{i=1}^{n} \sum_{j=1}^{n} \operatorname{Cov}\left(\frac{1}{n} x_{i}, \frac{1}{n} x_{j}\right) \\ &=\sum_{i=1}^{n} \operatorname{Cov}\left(\frac{1}{n} x_{i}, \frac{1}{n} x_{i}\right) \\ &=\sum_{i=1}^{n} \frac{1}{n^{2}} \operatorname{Cov}\left(x_{i}, x_{i}\right) =\sum_{i=1}^{n} \frac{1}{n^{2}} \operatorname{Cov}\left(x_{i}\right) =\sum_{i=1}^{n} \frac{1}{n^{2}} \operatorname{Cov}\left(X\right) \\ &=n \times \frac{1}{n^{2}} \Sigma =\frac{1}{n} \Sigma \end{aligned}
Cov(xˉ)=i=1∑nj=1∑nCov(n1xi,n1xj)=i=1∑nCov(n1xi,n1xi)=i=1∑nn21Cov(xi,xi)=i=1∑nn21Cov(xi)=i=1∑nn21Cov(X)=n×n21Σ=n1Σ
其中倒数第二行最后一个等号是因为它们都来自相同分布,因此协方差都相同。
令
A
=
∑
i
=
1
n
(
x
i
−
x
ˉ
)
(
x
i
−
x
ˉ
)
T
A=\sum_{i=1}^{n}\left(x_{i}-\bar{x}\right)\left(x_{i}-\bar{x}\right)^{T}
A=i=1∑n(xi−xˉ)(xi−xˉ)T
则:
E
(
A
)
=
E
[
∑
i
=
1
n
(
x
i
−
x
ˉ
)
(
x
i
−
x
ˉ
)
T
]
=
E
{
∑
i
=
1
n
[
(
x
i
−
μ
)
−
(
x
ˉ
−
μ
)
]
[
(
x
i
−
μ
)
−
(
x
ˉ
−
μ
)
]
T
}
=
E
{
∑
i
=
1
n
[
(
x
i
−
μ
)
(
x
i
−
μ
)
T
+
x
ˉ
x
ˉ
T
−
x
i
x
ˉ
T
+
x
i
μ
T
−
x
ˉ
x
i
T
+
μ
x
i
T
−
μ
μ
T
]
}
=
∑
i
=
1
n
E
[
(
x
i
−
μ
)
(
x
i
−
μ
)
T
]
+
E
[
n
x
ˉ
x
ˉ
T
−
(
∑
i
=
1
n
x
i
)
x
ˉ
T
+
(
∑
i
=
1
n
x
i
)
μ
T
−
x
ˉ
(
∑
i
=
1
n
x
i
T
)
+
μ
(
∑
i
=
1
n
x
i
T
)
−
n
μ
μ
T
]
=
∑
i
=
1
n
E
[
(
x
i
−
μ
)
(
x
i
−
μ
)
T
]
+
E
(
n
x
ˉ
x
ˉ
T
−
n
x
ˉ
x
ˉ
T
+
n
x
ˉ
μ
T
−
n
x
ˉ
x
ˉ
T
+
n
μ
x
ˉ
T
−
n
μ
μ
T
)
=
∑
i
=
1
n
E
[
(
x
i
−
μ
)
(
x
i
−
μ
)
T
]
−
n
E
(
x
ˉ
x
ˉ
T
−
x
ˉ
μ
T
−
μ
x
ˉ
T
+
μ
μ
T
)
=
∑
i
=
1
n
E
[
(
x
i
−
E
(
x
i
)
)
(
x
i
−
E
(
x
i
)
)
T
]
−
n
×
E
[
(
x
ˉ
−
μ
)
(
x
ˉ
−
μ
)
T
]
=
∑
i
=
1
n
Cov
(
x
i
)
−
n
×
E
[
(
x
ˉ
−
μ
)
(
x
ˉ
−
μ
)
T
]
=
n
×
Σ
−
n
×
E
[
(
x
ˉ
−
μ
)
(
x
ˉ
−
μ
)
T
]
=
n
×
Σ
−
n
×
E
[
(
x
ˉ
−
E
(
x
ˉ
)
)
(
x
ˉ
−
E
(
x
ˉ
)
)
T
]
=
n
×
Σ
−
n
×
Cov
(
x
ˉ
)
=
n
×
Σ
−
n
×
1
n
Σ
=
(
n
−
1
)
Σ
\begin{aligned} E(A) &=E\left[\sum_{i=1}^{n}\left(x_{i}-\bar{x}\right)\left(x_{i}-\bar{x}\right)^{T}\right] \\ &=E\left\{\sum_{i=1}^{n}\left[\left(x_{i}-\mu\right)-(\bar{x}-\mu)\right]\left[\left(x_{i}-\mu\right)-(\bar{x}-\mu)\right]^{T}\right\} \\ &=E\left\{\sum_{i=1}^{n}\left[\left(x_{i}-\mu\right)\left(x_{i}-\mu\right)^{T}+\bar{x} \bar{x}^{T}-x_{i} \bar{x}^{T}+x_{i} \mu^{T}-\bar{x} x_{i}^{T}+\mu x_{i}^{T}-\mu \mu^{T}\right]\right\} \\ &=\sum_{i=1}^{n} E\left[\left(x_{i}-\mu\right)\left(x_{i}-\mu\right)^{T}\right]+E\left[n \bar{x} \bar{x}^{T}-\left(\sum_{i=1}^{n} x_{i}\right) \bar{x}^{T}+\left(\sum_{i=1}^{n} x_{i}\right) \mu^{T}-\bar{x}\left(\sum_{i=1}^{n} x_{i}^{T}\right)+\mu\left(\sum_{i=1}^{n} x_{i}^{T}\right)-n \mu \mu^{T}\right] \\ &=\sum_{i=1}^{n} E\left[\left(x_{i}-\mu\right)\left(x_{i}-\mu\right)^{T}\right]+E\left(n \bar{x} \bar{x}^{T}-n \bar{x} \bar{x}^{T}+n \bar{x} \mu^{T}-n \bar{x} \bar{x}^{T}+n \mu \bar{x}^{T}-n \mu \mu^{T}\right) \\ &=\sum_{i=1}^{n} E\left[\left(x_{i}-\mu\right)\left(x_{i}-\mu\right)^{T}\right]-n E\left(\bar{x} \bar{x}^{T}-\bar{x} \mu^{T}-\mu \bar{x}^{T}+\mu \mu^{T}\right) \\ &=\sum_{i=1}^{n} E\left[\left(x_{i}-E\left(x_{i}\right)\right)\left(x_{i}-E\left(x_{i}\right)\right)^{T}\right]-n \times E\left[(\bar{x}-\mu)(\bar{x}-\mu)^{T}\right] \\ &=\sum_{i=1}^{n} \operatorname{Cov}\left(x_{i}\right)-n \times E\left[(\bar{x}-\mu)(\bar{x}-\mu)^{T}\right] \\ &=n \times \Sigma-n \times E\left[(\bar{x}-\mu)(\bar{x}-\mu)^{T}\right] \\ &=n \times \Sigma-n \times E\left[(\bar{x}-E(\bar{x}))(\bar{x}-E(\bar{x}))^{T}\right] \\ &=n \times \Sigma-n \times \operatorname{Cov}(\bar{x}) \\ &=n \times \Sigma-n \times \frac{1}{n} \Sigma \\ &=(n-1) \Sigma \end{aligned}
E(A)=E[i=1∑n(xi−xˉ)(xi−xˉ)T]=E{i=1∑n[(xi−μ)−(xˉ−μ)][(xi−μ)−(xˉ−μ)]T}=E{i=1∑n[(xi−μ)(xi−μ)T+xˉxˉT−xixˉT+xiμT−xˉxiT+μxiT−μμT]}=i=1∑nE[(xi−μ)(xi−μ)T]+E[nxˉxˉT−(i=1∑nxi)xˉT+(i=1∑nxi)μT−xˉ(i=1∑nxiT)+μ(i=1∑nxiT)−nμμT]=i=1∑nE[(xi−μ)(xi−μ)T]+E(nxˉxˉT−nxˉxˉT+nxˉμT−nxˉxˉT+nμxˉT−nμμT)=i=1∑nE[(xi−μ)(xi−μ)T]−nE(xˉxˉT−xˉμT−μxˉT+μμT)=i=1∑nE[(xi−E(xi))(xi−E(xi))T]−n×E[(xˉ−μ)(xˉ−μ)T]=i=1∑nCov(xi)−n×E[(xˉ−μ)(xˉ−μ)T]=n×Σ−n×E[(xˉ−μ)(xˉ−μ)T]=n×Σ−n×E[(xˉ−E(xˉ))(xˉ−E(xˉ))T]=n×Σ−n×Cov(xˉ)=n×Σ−n×n1Σ=(n−1)Σ
而样本协方差
S
S
S与
A
A
A关系为:
S
=
1
n
−
1
∑
i
=
1
n
(
x
i
−
x
ˉ
)
(
x
i
−
x
ˉ
)
T
=
1
n
−
1
A
\begin{aligned} S &=\frac{1}{n-1} \sum_{i=1}^{n}\left(x_{i}-\bar{x}\right)\left(x_{i}-\bar{x}\right)^{T} \\ &= \frac{1}{n-1} A \end{aligned}
S=n−11i=1∑n(xi−xˉ)(xi−xˉ)T=n−11A
因此,样本协协方差的期望为:
E
(
S
)
=
1
n
−
1
E
(
A
)
=
1
n
−
1
(
n
−
1
)
Σ
=
Σ
\begin{aligned} E(S) &= \frac{1}{n-1} E(A) \\ &= \frac{1}{n-1} (n-1) \Sigma = \Sigma \end{aligned}
E(S)=n−11E(A)=n−11(n−1)Σ=Σ
16.3
设X维数据规范化样本矩阵,则主成分分析等价于求解一下最优化问题:
min
L
∥
X
−
L
∥
F
s.t.
rank
(
L
)
≤
k
\begin{array}{c}\min _{L}\|X-L\|_{F} \\ \text { s.t. } \quad \operatorname{rank}(L) \leq k\end{array}
minL∥X−L∥F s.t. rank(L)≤k
这里
F
F
F是弗罗贝尼乌斯范数,
k
k
k为主成分甘薯。试问为什么?
首先PCA的求解完全可以用SVD方法进行,只不过对于原来的矩阵进行变形,但是:
min
L
∥
X
′
n
−
1
−
L
∥
F
\min _{L}\left\|\frac{X^{\prime}}{\sqrt{n-1}}-L\right\|_{F}
Lmin∥∥∥∥n−1X′−L∥∥∥∥F
与这个优化问题完全等价,剩下的就是为啥PCA或者SVD与求解这个最优化问题等价,而这个等价性,第二版书的第287页
,定理15.3
已经给出了阐述,并给出了证明,直接参阅即可。