多元统计分析其实是围绕实际生活中变量的多维特征产生的，主要是在以下场景有所应用:

  • 数据降维和结构简化
  • 分类与判别（聚类）
  • 变量间独立性的度量： 回归分析以及典型相关分析
  • 多元数据的统计推断：以多元正态分布为代表的统计推断问题

多元随机变量

多元统计中的多其实主要说的是维度的扩充，而不是变量个数的增减。其实在概统中已有了多元变量的分布函数的概念，但注意在概统中始终都是数量值函数，这里的多元则是多元向量值函数，从这个角度出发多元统计相当于是对一元的扩充，这点在多元正态分布的推断中体现会更加明显，统计量的抽样分布其实也是对一元的扩展，只不过因为引入多元变量之后，协方差函数会变得相对复杂，这时也就需要引入矩阵来进行简化，简单来说多元统计分析其实就是建立在矩阵基础之上的多元变量的统计分析。

分布函数

与一元分布函数类似,多元随机变量X=(X1,X2,⋯ ,Xp)X=(X_1,X_2,\cdots,X_p)的密度函数可以表示为:

f(x1,x2,⋯ ,xp)=f(x1,x2,⋯ ,xp)f(x_1,x_2,\cdots,x_p)=f(x_1,x_2,\cdots,x_p)

边缘分布的定义，对于p元随机向量，其前q个变量的边际分布函数定义为:

\begin{align*} F(x_1,x_2\cdots,x_q)&=P(X_1\leq x_1,X_2\leq x_2,\cdots,X_q\leq x_q,X_{q+1}<\infty,\cdots,X_p<\infty)\\ &=F(x_1,x_2\cdots,x_q,\infty)\\ \end{align*}

若求任意q个随机变量，只需适当进行转换即可

数字特征

在正式介绍多元统计分析之前，先对多元随机变量的数字特征进行描述:

期望

多元变量的期望其实就是对用矩阵表示多元向量的每一个元素求期望，这里取X=(X1,X2⋯ ,Xp)TX=(X_1,X_2\cdots,X_p)^T，于是可以得到:

E(X)=[E(X1)E(X2)⋮E(Xp)]=[∫−∞+∞x1f(x1)dx1∫−∞+∞x2f(x2)dx2⋮∫−∞+∞xpf(xp)dxp]\mathrm{E}(X)=\begin{bmatrix} \mathrm{E}(X_1)\\ \mathrm{E}(X_2)\\ \vdots\\ \mathrm{E}(X_p) \end{bmatrix}= \begin{bmatrix} \int_{-\infty}^{+\infty}x_1f(x_1)\mathrm{d}x_1\\ \int_{-\infty}^{+\infty}x_2f(x_2)\mathrm{d}x_2\\ \vdots\\ \int_{-\infty}^{+\infty}x_pf(x_p)\mathrm{d}x_p \end{bmatrix}

接着借助矩阵可以推倒多元随机变量期望的一些性质:

E(AXB)=AE(X)BE(AX+BY)=AE(X)+BE(Y)\mathrm{E}(AXB)=A\mathrm{E}(X)B\\ \mathrm{E}(AX+BY)=A\mathrm{E}(X)+B\mathrm{E}(Y)

其中A,BA,B为合适的常数矩阵，X,YX,Y为多元随机变量。

方差

类比一元中方差的定义，可以相似的定义多元随机变量的方差，不过需要注意的当变量变为了向量之后，多元随机向量的方差也变成了一个矩阵形式，沿用上边定义的X，可以得到:

D(X)=E[(X−E(X))(X−E(X))T]=[Var(X1)Cov(X1,X2)⋯Cov(X1,Xp)Cov(X2,X1)Var(X2)⋯Cov(X2,Xp)⋮⋮⋱⋮Cov(Xp,X1)Cov(Xp,X2)⋯Var(Xp)]\mathrm{D}(X)=E[(X-\mathrm{E}(X))(X-\mathrm{E}(X))^T]=\begin{bmatrix} \mathrm{Var}(X_1) & \mathrm{Cov}(X_1,X_2) & \cdots & \mathrm{Cov}(X_1,X_p)\\ \mathrm{Cov}(X_2,X_1) & \mathrm{Var}(X_2) & \cdots & \mathrm{Cov}(X_2,X_p)\\ \vdots & \vdots & \ddots & \vdots\\ \mathrm{Cov}(X_p,X_1) & \mathrm{Cov}(X_p,X_2) & \cdots & \mathrm{Var}(X_p) \end{bmatrix}

类似地可以定义两个随机向量的协方差:

Cov(X,Y)=E[(X−E(X))(Y−E(Y))T]=E(XYT)−E(X)E(Y)T\mathrm{Cov}(X,Y)=E[(X-\mathrm{E}(X))(Y-\mathrm{E}(Y))^T]=\mathrm{E}(XY^T)-\mathrm{E}(X)\mathrm{E}(Y)^T

其中X、Y为同型的列向量。
对于方差和协方差的性质有:

  • Var(X)\mathrm{Var}(X)为非负定矩阵
  • D(X+a)=D(X)\mathrm{D}(X+a)=\mathrm{D}(X)
  • D(AX)=AD(X)AT\mathrm{D}(AX)=A\mathrm{D}(X)A^T
  • Cov(AX,BY)=ACov(X,Y)BT\mathrm{Cov}(AX,BY)=A\mathrm{Cov}(X,Y)B^T
  • Cov(X+Y,Z)=Cov(X,Z)+Cov(Y,Z)\mathrm{Cov}(X+Y,Z)=\mathrm{Cov}(X,Z)+\mathrm{Cov}(Y,Z)

这些其实就是对期望计算与矩阵运算的灵活运用

多元正态分布

多元正态分布是一种经常会用到的假设分布，这里主要介绍多元正态分布的形式以及参数估计，统计推断部分不作为重点内容。

分布函数

多元正态分布的定义有好几种形式，比较常见的一种定义方式是借助密度函数,对于随机变量序列X=(X1,X2,⋯ ,Xp)TX=(X_1,X_2,\cdots,X_p)^T，若其密度函数为（注意xx是p×1p\times 1矩阵:

f(x)=1(2π)p/2∣Σ∣1/2exp⁡{−12(x−μ)TΣ−1(x−μ)}f(x)=\frac{1}{(2\pi)^{p/2}|\Sigma|^{1/2}}\exp\left\{-\frac{1}{2}(\boldsymbol{x}-\boldsymbol{\mu})^T\Sigma^{-1}(\boldsymbol{x}-\boldsymbol{\mu})\right\}

则称X满足多元正态分布，可以证明EX=μ,DX=Σ\mathrm{E}X=\mu,DX=\Sigma,因此一般简单记为:X∼Np(μ,Σ)X\sim N_p(\mu,\Sigma),考虑到这种定义方式建立在Σ\Sigma可逆的前提之下，因此另外一种定义方法为:
若X∼Nm(0,Im)X\sim N_m(\boldsymbol{0},\boldsymbol{I_m})，亦或者说X的每一个分量都服从均值为0，方差为1的正态分布，则可以得到:

Yp×1=Ap×mXm×1+μp×1Y_{p\times 1}=A_{p\times m}X_{m\times 1}+\mu_{p\times 1}

则Y∼Np(μ,AAT)Y\sim N_p(\mu,AA^T)，这种也是一种比较常用的方法，其实就是将多元正态分布看成既定数量个标准正态分布的线性组合。
最后，这里给出以特征函数定义的多元正态分布，对于随机向量X，若其特征函数为:

Φ(t)=exp⁡{itTμ−12tTΣt}\Phi(\boldsymbol{t})=\exp\left\{i\bold{t}^T\boldsymbol{\mu}-\frac{1}{2}\bold{t}^T\Sigma\bold{t}\right\}

则X∼Np(μ,Σ)X\sim N_p(\boldsymbol{\mu},\Sigma),这种定义方法也很好地回避了Σ\Sigma不可逆的问题。
多元正态分布的几条性质:

  • 独立和不相关等价。若Σ\Sigma为单位阵，则说明X的各个分量不相关，结合密度函数可以知道各个分量独立。若各个分量独立，则可知Σ\Sigma为对角阵，因此各分量不相关。
  • 线性可加。若X∼Np(μ,Σ)X\sim N_p(\boldsymbol{\mu},\boldsymbol{\Sigma})，则AX+b∼Np(Aμ+b,AΣAT)AX+b\sim N_p(A\boldsymbol{\mu}+b,A\Sigma A^T),其中AA为合适的产生常数矩阵
  • 边缘正态，多维正态分布的任意q个变量的边际分布仍服从多元正态分布(q<p)

边缘分布的均值和方差可以根据原来分布的期望方差得到（注意期望方差的实际含义，边际分布不改变两个变量的数字特征)

多元样本的数字特征

在用样本对总体进行估计时，首先要对总体进行抽样，假设抽取的第ii个样本记为Xi=(Xi1,⋯ ,Xip)X_i=(X_{i1},\cdots,X_{ip}),一般习惯将XiX_i作为行向量进行拼接，最后得到一个n×pn\times p的矩阵来代表样本矩阵，这里记作Xn×pX_{n\times p},在引入了样本矩阵之后，为了表示方便，这里给出样本方差和样本离差的定义:

\begin{equation*} \begin{aligned} \overline{X}&=（\sum_i \dfrac{X_{i1}}{n},\sum_i \dfrac{X_{i2}}{n},\cdots,\sum_i \dfrac{X_{ip}}{n}）^T\\ &=X^T \boldsymbol{1}_{n\times 1}\times \dfrac{1}{n} \end{aligned} \end{equation*}

接着定义样本离差阵S：

\begin{equation*} \begin{aligned} S&=\sum_i (X_i-\overline{X})(X_i-\overline{X})^T\\ &=X^TX-n\overline{X}^T\overline{X}\\ &=X^T(I_n-\frac{1}{n}\boldsymbol{1_{n\times 1}}\boldsymbol{1_{1\times n}})X \end{aligned} \end{equation*}

则样本协差阵为:

V=1nSV=\frac{1}{n}S

参数估计

定义了样本的离差阵之后，这里给出多元正态分布的mle:

\begin{equation*} \begin{aligned} L(\boldsymbol{\mu},\boldsymbol{\Sigma})=&\sum_{i=1}^n \log\left( \frac{1}{(2\pi)^{p/2}|\Sigma|^{1/2}}\exp\left\{-\frac{1}{2}(X_i-\boldsymbol{\mu})^T\Sigma^{-1}(X_i-\boldsymbol{\mu})\right\}\right)\\ =&-\frac{n}{2}\log|\Sigma|-\frac{1}{2}\sum_{i=1}^n(X_i-\boldsymbol{\mu})^T\Sigma^{-1}(X_i-\boldsymbol{\mu})\\ =&-\frac{n}{2}\log|\Sigma|-\frac{n}{2}\left\{(\overline{\boldsymbol{X}}-\boldsymbol{\mu})^T\Sigma^{-1}(\overline{\boldsymbol{X}}-\boldsymbol{\mu})\right\}-\frac{1}{2}\sum_{i=1}^n\left\{(X_i-\overline{\boldsymbol{X}})^T\Sigma^{-1}(X_i-\overline{\boldsymbol{X}})\right\}\\ \end{aligned} \end{equation*}

由上易知总体均值的mle为样本均值,对于总体协差阵的mle，只需考虑第一项和最后一项的取值，于是有:

\begin{equation*} \begin{aligned} L(\overline{X},\Sigma)=&-\frac{n}{2}\log|\Sigma|-\frac{1}{2}\sum_{i=1}^n\left\{(X_i-\overline{\boldsymbol{X}})^T\Sigma^{-1}(X_i-\overline{\boldsymbol{X}})\right\}\\ =&-\frac{n}{2}\log|\Sigma|-\frac{1}{2}\sum_{i=1}^n\left\{tr\left[\Sigma^{-1}(X_i-\overline{\boldsymbol{X}})(X_i-\overline{\boldsymbol{X}})^T\right]\right\}\\ =&-\frac{n}{2}\log|\Sigma|-\frac{1}{2}tr\left[\Sigma^{-1}\sum_{i=1}^n(X_i-\overline{\boldsymbol{X}})(X_i-\overline{\boldsymbol{X}})^T\right]\\ =&-\frac{n}{2}\log|\Sigma|-\frac{1}{2}tr\left[\Sigma^{-1}S\right]\\ =&-\frac{n}{2}\log|\Sigma|-\frac{n}{2}tr\left[\Sigma^{-1}V\right]\\ \end{aligned} \end{equation*}

若取A=Σ−1VA=\Sigma^{-1}V，则可得上式为(log⁡∣A∣−log⁡∣V∣−tr(A))×n2\left(\log|A|-\log|V|-tr(A)\right)\times \frac{n}{2}，考虑到∣A∣=∏λi,tr(A)=∑λi|A|=\prod \lambda_i,tr(A)=\sum \lambda_i，而函数log⁡x−x\log x-x在x=1x=1处取得最大值，因此当A=IA=I时，L(X‾,Σ)L(\overline{X},\Sigma)取得最大值，即Σ−1V=I\Sigma^{-1}V=I，所以Σ^=V\hat{\Sigma}=V，于是有多元正态分布总体参数的mle为：

\begin{equation*} \begin{aligned} \hat{\boldsymbol{\mu}}&=\overline{\boldsymbol{X}}\\ \hat{\boldsymbol{\Sigma}}&=V\\ \end{aligned} \end{equation*}

考虑到总体协差阵的估计不无偏，修正后得到新的统计量:

Σ^=Sn−1\hat{\boldsymbol{\Sigma}}=\frac{S}{n-1}

可以证明上述两个估计量是总体参数的无偏一致有效估计

抽样分布

类似地，接下来将主要讨论样本均值和样本协差阵的抽样分布，由多元正态分布的线性可加性，很容易得到样本均值X‾∼Np(μ,Σn)\boldsymbol{\overline{X}}\sim N_p(\boldsymbol{\mu},\dfrac{\Sigma}{n}),而对于SS的分布，这里先看一维的情况，在一维正态总体下，样本方差的分布为(n−1)S2σ2∼χ2(n−1)\dfrac{(n-1)S^2}{\sigma^2}\sim \chi^2(n-1)，而对于多元正态分布，则相当于把χ2\chi^2分布做了推广，对于多元正态分布来说，SS服从Wishart分布,这里先介绍一下Wishart分布。

Wishart分布

假设Xi=(Xi1,Xi2,⋯ ,Xip)T∼Np(μi,Σ)(i=1,2⋯ ,n)X_i=(X_{i1},X_{i2},\cdots,X_{ip})^T\sim N_p(\mu_i,\Sigma)(i=1,2\cdots,n)，记:

Wp×p=∑inXiXiT\boldsymbol{W}_{p\times p}=\sum_i^n X_iX_i^T

则称W\boldsymbol{W}服从自由度为n的Wishart分布，简记为W∼Wp(n,Σ,Z)\boldsymbol{W}\sim W_p(n,\Sigma,Z)，其中ZZ为:

Z=(μ1,μ2,⋯ ,μn)(μ1,μ2,⋯ ,μn)T=∑inμiμiTZ=\left(\mu_1,\mu_2,\cdots,\mu_n \right)(\mu_1,\mu_2,\cdots,\mu_n)^T=\sum_i^n\mu_i\mu_i^T

若Z=0Z=\bold{0}，则可简记为W∼Wp(n,Σ)\boldsymbol{W}\sim W_p(n,\Sigma)
Wishart分布的概率密度函数为:

f(w)={∣w∣12(n−p−1)exp⁡{−12tr⁡(Σ−1w)}2np2πp(p−1)4∣Σ∣n2∏i=1pp(n−i−12),∣w∣>00, 其它 f(w)= \begin{cases}\dfrac{|w|^{\frac{1}{2}(n-p-1)}\exp \left\{-\dfrac{1}{2} \operatorname{tr}\left(\Sigma^{-1} w\right)\right\}}{2^{\frac{n p}{2}} \pi^{\frac{p(p-1)}{4}}|\Sigma|^{\frac{n}{2}} \prod_{i=1}^p p\left(\frac{n-i-1}{2}\right)}, & |w|>0 \\ 0, & \text { 其它 }\end{cases}

接下来不加证明的给出一些常用的性质:

  • 若总体X∼Np(μ,Σ)X\sim N_p(\mu,\Sigma)，则样本协方差阵SS服从Wp(n−1,Σ)W_p(n-1,\Sigma)
  • 可加性,若W1∼Wp(n1,Σ1),W2∼Wp(n2,Σ2)\boldsymbol{W}_1\sim W_p(n_1,\Sigma_1),\boldsymbol{W}_2\sim W_p(n_2,\Sigma_2)，且W1,W2\boldsymbol{W}_1,\boldsymbol{W}_2相互独立，则W1+W2∼Wp(n1+n2,Σ1+Σ2)\boldsymbol{W}_1+\boldsymbol{W}_2\sim W_p(n_1+n_2,\Sigma_1+\Sigma_2)
  • Wp×p∼Wp(n,Σ)W_{p\times p}\sim W_p(n,\Sigma)，CC为非奇异阵，则CWCT∼Wp(n,CΣCT)CWC^T\sim W_p(n,C\Sigma C^T)

除了Wishart分布，在假设检验中还有几种常用的分布，这里也做以说明：
在方差已知的情况下，可以构建(X‾−μ0)Σ−1(X‾−μ0)T\left(\overline{X}-\mu_0\right)\Sigma^{-1}\left(\overline{X}-\mu_0\right)^T作为枢轴量，可证明该枢轴量服从χ2(p)\chi^2(p)分布。而如果Σ\Sigma未知，那对总体均值的检验就会存在困难。类比一维情况下，我们引入了t统计量，t统计量其实就相当于是用样本方差(s2n−1)(\frac{s^2}{n-1})来估计总体方差σ2n\frac{\sigma^2}{n}，在多元正态分布中，其实也有类似的想法，前边提到过Sn−1\frac{S}{n-1}是总体方差Σ\Sigma的无偏估计，因此也可以很自然地来构造一个具有下列形式的统计量:

\begin{align} T^2&=\sqrt{n}\left(\overline{X}-\mu_0\right)^T\left(\frac{S}{n-1}\right)^{-1}\sqrt{n}\left(\overline{X}-\mu_0\right)\\ &=(n-1)\sqrt{n}\left(\overline{X}-\mu_0\right)^TS^{-1}\sqrt{n}\left(\overline{X}-\mu_0\right)\\ \end{align}

可以发现上述统计量其实是一维情况下的t统计量的平方，该统计量在多元正态分布假设检验中常用到，一般记为T2T^2分布，接下来给出该分布的定义:

T2分布

若X∼Np(μ,Σ)X\sim N_p(\mu,\Sigma)，S∼Wp(n,Σ)S\sim W_p(n,\Sigma),且有X,SX,S相互独立，则有:

T2=nXTS−1X∼T2(p,n)T^2=nX^TS^{-1}X\sim T^2(p,n)

该分布成为HotellingT2T^2分布,该统计量首先由Harold Hotelling给出，因此得名。
引入了该分布之后，其实可以从另外一个角度来理解刚才提到的检验统计量:

(n−1)(n(X‾−μ0))TS−1(n(X‾−μ0))(n-1)\left(\sqrt{n}\left(\overline{X}-\mu_0\right)\right)^TS^{-1}\left(\sqrt{n}\left(\overline{X}-\mu_0\right)\right)

根据前面的内容，我们易知X‾−μ0∼Np(0,Σn)\bold{\overline{X}-\mu_0}\sim N_p(0,\dfrac{\Sigma}{n})，而样本协差阵S∼Wp(n−1,Σ)S\sim W_p(n-1,\Sigma)，因此n相当于对X‾−μ0\bold{\overline{X}-\mu_0}的方差作修正，以使得X‾−μ0\bold{\overline{X}-\mu_0}与SS满足T2T^2分布的条件。
最后，关于T2T^2分布的密度函数，联想一元回归中存在F=t2F=t^2的分布，其实可以从T2T^2的分布类比到F分布，结果如下:

F=n−p+1npT2∼F(p,n−p+1)F=\frac{n-p+1}{np}T^2\sim F(p,n-p+1)

其中T2∼T2(p,n)T^2\sim T^2(p,n)。因此就可以利用T2T^2来作假设检验

wilks分布

wilks分布是在检验扩展方差分析的概念时提出的，常被用作多个正态总体均值的假设检验。当变量扩展到多元之后，变量的方差变成了一个协方差矩阵，矩阵本身的大小不可比为我们度量总体方差的大小带来了困难，进而就有学者提出了用协方差阵的行列式或者迹方法来度量协差阵的分散程度，接着就用到了wilks分布，该分布可以看成两个独立wishart分布行列式值的比:
若A1∼Wp(n1,Σ)A_1\sim W_p(n_1,\Sigma),A2∼Wp(n2,Σ)A_2\sim W_p(n_2,\Sigma),n1>pn_1>p且A1,A2A_1,A_2相互独立，则称下列分布服从Wilks分布:

Λ=∣A1∣∣A1+A2∣\Lambda=\dfrac{\left|A_1\right|}{\left|A_1+A_2\right|}

简记为Λ∼Λ(p,n1,n2)\Lambda \sim \Lambda\left(p, n_1, n_2\right) 其中n1,n2n_1, n_2 为自由度。实际应用中，常将该统计量转化为T2T^2统计量来寻找密度函数。

假设检验

关于假设检验部分的内容，这里只强调一点，严格来说，多元正态分布的假设检验可以同时拆成p个变量各自的假设检验，但是这种检验方法其实并没有考虑到变量之间的相关关系，这也是为什么要展开多元分析的原因。其实从刚才提到的几个检验统计量中也能看出来，虽然随机变量本身是多维的，但是构造的统计量都是一维的，这其实某种程度也保证了检验结果的一致性。多元统计分析相当于是对一维情况的一个扩充，但是注意扩充的是样本的维数，而不是统计量的维数。

主成分分析

主成分分析即Principle Component Analysis，是一种广泛应用的方法，在多元统计分析中，会用到PCA的几种场景是:

  • 作为一种降维方法，捕捉关键信息
  • 编制指数
  • 消除原始回归自变量之间的相关性，然后对一些主成分进行回归（称为主成分回归）
  • 作综合评价，产生客观权重（选择若干个正向/负向指标，选出若干个主成分，以方差贡献率为权重计算评分，需要注意的是指标必须都为正向或者负向，进而得到的加权分数有比较明确的经济意义）
  • 用于处理现代更加复杂的数据，如可将其用于分析具有无穷维特征的函数性数据

模型

PCA（Principal Component Analysis）是一种常用的线性降维方法。核心思想是对特征空间进行重构，PCA 通过线性变换，将 N 维空间的原始数据变换到一个较低的RR维空间(R<N(R<N)：

  • 在降维过程中，不可避免的要造成信息损失。如原来在高维空间可分的点，在低维空间可能变成一个点，变得不可分。因此，要在降维过程中尽量减少这种损失（最小化重构损失）。
  • 为使样本投影到低维空间后尽可能分散，它们的方差要尽可能大（最大化投影方差）。

这就构成了 PCA 的基本思想,其实本质上来说这两种方法是等价的。

最大化投影方差

对于最大化投影方差有:

\begin{align} Argmin&\sum_{i=1}^N w^Tx_i x_i^Tw\\ &=w^T\sum_{i=1}^N x_i x_i^Tw\\ &=w^T\Sigma w\\ s.t.&w^Tw=1\\ x_i:\text{中心化后样本向量} \end{align}

使用拉格朗日乘子法求解上述问题，得到：

Σw=λw\Sigma w=\lambda w

即投影方向为协方差阵的特征向量，当选择 q 个方向时，注意到优化目标最后可以替换成∑qλi\sum^q \lambda_i，因此选择 q 个最大的特征值对应的特征向量作为投影方向，即可实现降维。

特征向量之间线性无关保证了降维后的特征之间也是线性无关的，这是 PCA 降维的一个重要性质。

最小化重构损失

最大化投影方差从另外一个角度理解可以看成是最小化重构损失,即我们希望降维之后损失的信息尽可能少，假设 PCA 方法计算出来的所有特征向量有 p 个(μ1,μ2,…,μp\mu_1,\mu_2,\dots,\mu_p)，我们选取前 q 个向量作为我们的降维向量，最小化损失信息有:

Argmin∑i=1N∣∣∑k=1p(xi⋅μk)μk−∑k=1q(xi⋅μk)μk∣∣2=∑i=1N∣∣∑k=q+1p(xi⋅μk)μk∣∣2\begin{aligned} Argmin\quad &\sum_{i=1}^N ||\sum^p_{k=1} (x_i\cdot \mu_k)\mu_k-\sum^q_{k=1} (x_i\cdot \mu_k)\mu_k||^2\\ &=\sum_{i=1}^N ||\sum^p_{k=q+1} (x_i\cdot \mu_k)\mu_k||^2\\ \end{aligned}

一般为简便运算会先将样本中心化。因为所有的向量是线性无关的，所以可以把上述平方进行拆解，最后可以得到最小化目标函数变成了:

Argmin∑k=q+1pμkSμk\begin{aligned} Argmin\quad & \sum^p_{k=q+1} \mu_k S\mu_k\\ \end{aligned}

其实就是找到最小的 p-q 个特征向量，这也证明了最小化重构损失其实就是最大化投影方差
具体过程为选取一组 N 个 R 维的正交基组成的矩阵 P，然后令 P 左乘数据集 X 得到变换后的数据集的 X’，进而实现了数据集的维数由 N 变换为 R（R<NR<N）

机器学习中也常常会用到PCA作降维，机器学习中侧重于将PCA作为一种降维方法，而对于统计学来说，我们可能会关注主成分个数的选择问题，接下来将对几种经常用的方法进行叙述:

主成分个数的选择

关于PCA中最佳的主成分个数的选择，这里介绍几种多元统计分析中常用的方法：

  • 绝对值筛选方法，根据特征值的大小决定保留的主成分个数，一般保留特征值大于1的主成分，这种方法常用的一种可视化方法是绘制碎石图(screeplot)
  • 相对值筛选法，按照方差贡献率进行筛选，一般保留累计方差贡献率大于90%的主成分
  • Kaiser’s rule：计算特征值的统计量，过滤掉小于该统计量的所有特征值。常用的是计算特征值的平均值，但是普遍认为这样会导致筛除较少的主成分，倾向于对平均值乘一个小于1的因子再截断，从这视角出发，可以认为Kaiser’s rule是相对值筛选法的一种改进
  • 借助统计检验，建立假设检验H0:λk+1=⋯=λpH_0:\lambda_{k+1}=\cdots=\lambda_p（检验特征值是否存在显著差异），其实就是检验最后几个主成分是否存在显著差异（注意不是直接检验是否为0，严格来说p维变量的p个主成分不存在为0的特征值，为0说明总体的维数小于p)

因子分析

因子分析方法与PCA方法比较类似，它的核心思想是希望能找到n个样本或者p个变量的若干个共同因子(通常小于样本数或者变量数)，这些共同因子能够解释原始数据中的大部分信息，从而达到降维的目的,从这个角度出发，因子分析和PCA有许多共通之处，但是其实因子分析更关心的则是变量之间的相关结构，因子分析模型一般可以写作(以R型因子分析为例):

{X1=a11F1+a12F2+…+a1mFm+ε1X2=a21F1+a22F2+…+a2mFm+ε2⋯Xp=ap1F1+ap2F2+…+apmFm+εp\left\{\begin{array}{c} X_1=a_{11} F_1+a_{12} F_2+\ldots+a_{1 m} F_m+\varepsilon_1 \\ X_2=a_{21} F_1+a_{22} F_2+\ldots+a_{2 m} F_m+\varepsilon_2 \\ \quad \cdots \\ X_p=a_{p 1} F_1+a_{p 2} F_2+\ldots+a_{p m} F_m+\varepsilon_p \end{array}\right.

用矩阵表示

[X1X2⋯Xp]=[a11a12…a1ma21a22…a2m⋯ap1ap2…apm][F1F2⋯Fm]+[ε1ε2⋯εp]\left[\begin{array}{l} X_1 \\ X_2 \\ \cdots \\ X_p \end{array}\right]=\left[\begin{array}{ll} a_{11} & a_{12} \ldots a_{1 m} \\ a_{21} & a_{22} \ldots a_{2 m} \\ \cdots & \\ a_{p 1} & a_{p 2} \ldots a_{p m} \end{array}\right]\left[\begin{array}{l} F_1 \\ F_2 \\ \cdots \\ F_m \end{array}\right]+\left[\begin{array}{l} \varepsilon_1 \\ \varepsilon_2 \\ \cdots \\ \varepsilon_p \end{array}\right]

简记为

X(p×1)=AF(p×m)(m×1)+ε(p×1)\underset{(p \times 1)}{X}=\underset{(p \times m)(m \times 1)}{AF}+\underset{(p \times 1)}{\varepsilon}

对模型的一些基本假设:

  • m≤pm\le p
  • F叫做共同因子，D(F)=ImD(F)=I_m，即m个因子不相关且方差均为1
  • ε\varepsilon叫做特殊因子，D(ε)=diag(σ12,⋯ ,σp2)D(\varepsilon)=diag(\sigma_1^2,\cdots,\sigma_p^2)，即p个随机变量不相关，且方差不相同
  • Cov(F,ε)=0\mathrm{Cov}(F,\varepsilon)=\bold{0}

其中矩阵A叫做因子载荷矩阵,关于这个矩阵的名字来源可结合每个元素的实际含义来看：

\begin{equation*} \begin{aligned} X_i&=a_{i1} F_1+a_{i2} F_2+\ldots+a_{i m} F_m+\varepsilon_i\\ \mathrm{Cov}(X_i,F_j)&=\sum_i a_{ij}\mathrm{E}(F_iF_j)+\mathrm{E}(\varepsilon_iF_j)\\ &=a_{ij}\\ \mathrm{Var}(X_i)&=\sum_j a_{ij}^2+\sigma_i^2\\ \end{aligned} \end{equation*}

即aija_{ij}表示第i个变量与第j个因子的相关系数。
∑jaij2\sum_j a^2_{ij}表示共同因子对于XiX_i方差的贡献度，也叫做变量共同度；σi2\sigma_i^2则表示第i个变量的剩余方差，即由特定因子产生的方差。
最后，定义∑iaij2\sum_i a^2_{ij}表示第j个因子对X的所有变量的总影响，也叫做公共因子的方差贡献，可以用来衡量公共因子的相对重要性。

关于因子模型，需要注意以下两点:

  • 因子模型不受量纲的影响
  • 因子模型的因子载荷矩阵不唯一（因子载荷阵和因子同时作正交变换可得到新的因子和新的因子载荷矩阵，详细见高惠璇老师的应用多元分析）

模型估计

因子模型的模型估计主要是对因子载荷阵的估计，其中常用的一种方法是主成分方法，接下来将主要对该方法加以叙述:
对于X=AF+εX=AF+\varepsilon，不考虑特殊因子的方差我们有:

Var(X)=AAT\mathrm{Var}(X)=AA^T

在主成分分析法下，主要是对样本协方差阵作特征值分解(考虑到样本协差阵为n阶实对称矩阵，因此一定可以对角化):

Σ=PΛPT\Sigma=P\Lambda P^T

将Λ\Lambda拆分成diag{λ1,⋯ ,λp}diag\{\sqrt{\lambda_1},\cdots,\sqrt{\lambda_p}\}可以得到A的一个估计:

A^=PΛ12=(λ1e1,⋯ ,λpep)\hat{A}=P\Lambda^{\frac{1}{2}}=\left(\sqrt{\lambda_1}e_1,\cdots,\sqrt{\lambda_p}e_p\right)

可以发现因子分解得到的因子载荷真相当于在PCA主成分求解的系数向量上乘上对应特征值的12\frac{1}{2}次方
关于mle和主因子解法，可以参考高惠璇老师的书。

因子的正交旋转

假设存在正交矩阵P,则有:

X=AF+ε=APTPF+εX=AF+\varepsilon=AP^TPF+\varepsilon

这也是因子载荷矩阵不唯一的主要原因。

注意这里需要严谨证明新的因子模型满足因子模型的所有假设

因为因子载荷矩阵的估计结果不唯一，那么我们就希望通过正交变换得到一个比较“理想”的因子载荷阵，使得因子载荷阵每一列的取值向着尽可能两极化的方向取值，这也就能看出那些变量与哪些因子的相关性比较大，进而达到结构简化的目的。最后优化目标变成了最小化因子载荷阵每一列取值方差的和

对于列取值差异较大的情况，可以先对每行因子除上其共同度

因因子的旋转过程过于繁琐，这里不再赘述。

因子得分

在得出了因子载荷阵后，我们其实就用因子把样本线性表示出来。如果还想进一步求解因子的值，或者说用变量把因子线性表示，这里就需要对因子得分进行估计。

对F的估计与一般的参数估计是有区别的，F本身是一个随机变量
常用的估计方法有加权最小二乘法和回归法两种方法

gls

当因子载荷矩阵的估计结果得到之后，因子模型可以看做是一个以因子向量F为权重的线性回归，因为每个特殊因子的方差不同所以需要先消除异方差再作ols估计即可。

前提是特殊因子的方差是已知的

回归法

Fj=bj1X1+⋯+bjpXp+εj(j=1,⋯ ,m)F_j=b_{j 1} X_1+\cdots+b_{j p} X_p+\varepsilon_j \quad(j=1, \cdots, m)

回归法其实就是以F为因变量作回归，然后考虑到F是未知的，所以对回归模型两边乘上XiX_i来建立方程，接着再进行估计即可。该方法由汤普森提出

注意这里是对n个回归方程进行估计，建立估计方程的过程有点类似于时序中的矩估计方法。