第16章 主成分分析

习题16.1

  对以下样本数据进行主成分分析：

解答：

解答思路：

1. 给出主成分分析算法
2. 自编程基于numpy库实现主成分分析算法

解答步骤：

第1步：主成分分析算法

  主成分分析方法主要有两种，可以通过相关矩阵的特征值分解或样本矩阵的奇异值分解进行。

1. 相关矩阵的特征值分解算法

  根据书中第16.2.2节的相关矩阵的特征值分解算法具体步骤：

（1）对观测数据按照如下公式进行规范化处理，得到规范化数据矩阵，仍以$X$表示

$$$$

（2）依据规范化数据矩阵，计算样本相关矩阵$R$

$$$$

（3）求样本相关矩阵$R$的$k$个特征值和对应的$k$个单位特征向量。 求解$R$的特征方程

$$$$

求方差贡献度$\displaystyle \sum_{i=1}^k \eta_i$达到预定值的主成分个数$k$。
求前$k$个特征值对应的单位特征向量

$$$$

（5）计算$k$个主成分$y_i$与原变量$x_i$的相关系数$\rho(x_i,y_i)$，以及$k$个主成分对原变量$x_i$的贡献率$v_i$。

（6）计算$n$个样本的$k$个主成分值
将规范化样本数据带入$k$个主成分式，得到$n$个样本的主成分值。第$j$个样本$x_j = (x_{1j}, x_{2j}, \cdots, x_{mj})^T$的第$i$个主成分值是

$$$$

i = 1,2,\cdots,m, \quad j=1,2,\cdots, n

$k$个主成分$y_1,y_2,\cdots,y_k$的累计方差贡献率定义为$k$个方差之和与所有方差之和的比

$$$$

$X'$每一列的均值为零。
（2）对矩阵$X'$进行截断奇异值分解，得到

$$$$

第2步：自编程基于numpy库实现主成分分析算法

import numpy as np

def pca_svd(X, k):
    """
    样本矩阵的奇异值分解的主成分分析算法
    :param X: 样本矩阵X
    :param k: 主成分个数k
    :return: 特征向量V，样本主成分矩阵Y
    """
    n_samples = X.shape[1]

    # 构造新的n×m矩阵
    T = X.T / np.sqrt(n_samples - 1)

    # 对矩阵T进行截断奇异值分解
    U, S, V = np.linalg.svd(T)
    V = V[:, :k]

    # 求k×n的样本主成分矩阵
    return V, np.dot(V.T, X)

X = np.array([[2, 3, 3, 4, 5, 7],
              [2, 4, 5, 5, 6, 8]])
X = X.astype("float64")

# 规范化变量
avg = np.average(X, axis=1)
var = np.var(X, axis=1)
for i in range(X.shape[0]):
    X[i] = (X[i, :] - avg[i]) / np.sqrt(var[i])

# 设置精度为3
np.set_printoptions(precision=3, suppress=True)
V, vnk = pca_svd(X, 2)

print("正交矩阵V：")
print(V)
print("样本主成分矩阵Y：")
print(vnk)

正交矩阵V：
[[ 0.707  0.707]
 [ 0.707 -0.707]]
样本主成分矩阵Y：
[[-2.028 -0.82  -0.433  0.     0.82   2.461]
 [ 0.296 -0.046 -0.433  0.     0.046  0.137]]

习题16.2

  证明样本协方差矩阵$S$是总体协方差矩阵$\Sigma$的无偏估计。

解答：

解答思路：

1. 给出总体协方差矩阵$\Sigma$
2. 给出样本协方差矩阵$S$
3. 给出无偏估计的定义
4. 证明样本协方差矩阵$S$是总体协方差矩阵$\Sigma$的无偏估计

解答步骤：

第1步：总体协方差矩阵$\Sigma$

  根据书中第16章的定理16.1：

定理16.1 设$x$是$m$维随机变量，$\Sigma$是$x$的协方差矩阵，$\Sigma$的特征值分别是$\lambda_1 \geqslant \lambda_2 \geqslant \cdots \geqslant \lambda_m \geqslant 0$，特征值对应的单位特征向量分别是$\alpha_1, \alpha_2, \cdots, \alpha_m$，则$x$的第$k$个主成分是

$$$$

即协方差矩阵$\Sigma$的第$k$个特征值。

  根据书中第16章的本章概要的总体主成分的性质：

主成分$y$的协方差矩阵是对角矩阵

$$$$

样本协方差矩阵$S$为

$$$$

S = [s_{ij}]{m \times m} \ \displaystyle s{ij} = \frac{1}{n - 1} \sum_{k=1}^n (x_{ik} - \bar{x}i)(x{jk} - \bar{x}_j), \quad i,j=1,2,\cdots,m \end{array} \tag{16.41}

$z_1, \cdots, z_{n-1}$独立同分布于$N_p(0, \Sigma)$，于是

$$$$

则$S$是$\Sigma$的无偏估计。在实用中普遍采用$S$来估计$\Sigma$。

  无偏估计是用样本统计量来估计总体参数的一种无偏推断。估计量的期望等于被估计参数的真实值，则称此估计量为被估计参数的无偏估计，即具有无偏性。表示在多次重复下，它们的平均数接近所估计的参数真值。

第4步：证明样本协方差矩阵$S$是总体协方差矩阵$\Sigma$的无偏估计

  由样本协方差矩阵$S$定义：

  令：

  根据无偏估计的定义，证明样本协方差矩阵$S$是总体协方差矩阵$\Sigma$的无偏估计，即证明样本协方差矩阵$S$的期望等于总体协方差矩阵$\Sigma$：

即证明：

  参考《应用多元统计分析 第五版》（王学民，2017年8月第5版）书中第51页的无偏性一节的推导

依(2.2.14)式

$$$$

  因此，样本协方差矩阵$S$是总体协方差矩阵$\Sigma$的无偏估计。

习题16.3

  设$X$为数据规范化样本矩阵，则主成分等价于求解以下最优化问题：

其中，$F$是弗罗贝尼乌斯范数，$k$是主成分个数。试问为什么？

解答：

解答思路：

1. 给出弗罗贝尼乌斯范数的定义
2. 给出矩阵的最优近似
3. 证明主成分是该最优化问题的最优解，该最优化问题的最优解是$X$的主成分

解答步骤：

第1步：弗罗贝尼乌斯范数

  根据书中第15章的定义15.4的弗罗贝尼乌斯范数：

  定义15.4（弗罗贝尼乌斯范数） 设矩阵$A \in R^{m \times m}$，$A = [a_{ij}]_{m \times m}$，定义矩阵$A$的弗罗贝尼乌斯范数为

$$$$

则

$$$$

则

$$$$

X' = \frac{1}{\sqrt{n-1}}X^T

\begin{array}{cl} \min \limits_L &\displaystyle \left |\frac{1}{\sqrt{n-1}} X^T - \frac{1}{\sqrt{n-1}} L^T \right |_F \ \text{s.t.} & \text{rank}(L) \leqslant k \end{array}

X' = \frac{1}{\sqrt{n-1}}X^T, L' = \frac{1}{\sqrt{n-1}}L^T

\begin{array}{cl} \min \limits_L & |X' - L' |_F \ \text{s.t.} & \text{rank}(L) \leqslant k \ & \displaystyle X' = \frac{1}{\sqrt{n-1}}X^T \ & \displaystyle L' = \frac{1}{\sqrt{n-1}}L^T \end{array}

X' = U' \Sigma' V'^T

A \sim X' \ S \sim L'