多元统计分析03：多元正态分布的参数估计

• 2021 年 11 月 4 日
• 筆記
• 《多元统计分析》学习笔记

Chapter 3 多元正态分布的参数估计

一、随机阵的正态分布

Part 1：随机阵及其运算

从这里开始我们讨论随机阵的问题。把来自 \(p\)​ 元总体的容量为 \(n\)​​​ 的简单随机样本排成一个矩阵，就得到了样本数据阵。这是一个随机阵，其定义如下：

数据阵的每一行 \(X_{(i)}’\) 都是随机向量 \((X_1,X_2,\cdots,X_p)\) 的一个简单随机样本；

数据阵的每一列 \(\mathcal{X}_j\)​ 都是随机变量 \(X_j\)​​ 的一组简单随机样本。

拉直运算指的是将随机矩阵转化为一个长的列向量，把 \(X\) 中的第 \(2\) 列接到第 \(1\) 列的后面，再把第 \(3\) 列接到第 \(2\) 列的后面，以此类推。

如果把样本数据阵写成 \(p\) 个列向量的形式，即 \(X=\left(\mathcal{X}_1,\mathcal{X}_2,\cdots,\mathcal{X}_p\right)\) ，则拉直运算就是把矩阵的每一个列向量按列排列，组成一个 \(np\) 维向量，记为

如果要对样本进行拉直（按行拉直），可以先将数据阵转置，然后进行拉直运算，组成一个 \(np\)​ 维向量，记为

特别地，如果 \(X\) 是 \(p\) 阶对称随机阵，在 \(X\) 中只包含 \(p(p+1)/2\) 个不同的随机变量，故将其直接进行拉直运算，拉直成一个 \(p^2\) 维向量是不合适的。因此，我们专门定义了对称矩阵的拉直运算，将 \(\rm X\) 拉直成一个 \(p(p+1)/2\)​ 维向量，即

克罗内克(Kronecker)积又称为矩阵的直积，其运算法则简单来说就是用左矩阵的每一个元素去数乘右矩阵，其定义如下：

设 \(A=(a_{ij})\)​​ 是 \(n\times p\)​ 的矩阵，\(B\) 是 \(m\times q\)​ 的矩阵，定义 \(A\) 和 \(B\) 的克罗内克积为

可以看出，\(A\otimes B\) 的每个元素都是一个矩阵，该矩阵为 \(A\) 的对应元素数乘 \(B\) 得到。如果将 \(A\otimes B\) 的每个元素上的矩阵写开，将得到一个 \(mn\times pq\)​​ 维的矩阵。注意：\(A\otimes B\neq B\otimes A\) 。

Part 2：随机阵的正态分布

接下来我们考虑样本数据阵的分布。如果样本来自多元正态总体 \(N_p(\mu,\Sigma)\)​ ，那么样本数据阵 \(X\)​ 的每一列都是来自一元正态总体的简单随机样本，所以是相互独立的。

根据按行拉直运算的定义，\({\rm Vec}\left(X’\right)\) 指的是将每个样本排列在一起拉直得到的列向量，所以有

其中 \(\bold{1}_n\) 表示向量元素均为 \(1\) 的 \(n\) 维常向量，\(I_n\) 表示 \(n\) 阶单位矩阵。根据克罗内克积的定义，

这样我们就可以定义随机阵的正态分布。如果一个随机矩阵 \(X\)​ 按样本拉直后满足

就称 \(X\) 服从矩阵正态分布，记作

其中

容易验证

于是随机阵的正态分布可以等价的表示为

随机阵的正态分布具有如下性质：设 \(X\sim N_{n\times p}(M,I_n\otimes\Sigma)\) ，设 \(A\) 是 \(k\times n\) 常数矩阵，\(B\) 是 \(q\times p\) 常数矩阵，\(D\) 是 \(k\times q\) 常数矩阵，如果对 \(X\) 作线性变换得到 \({\rm Z}=AXB’+D\) ，则有

二、多元正态分布的参数估计

Part 1：基本统计量

设总体 \(X=(X_1,X_2,\cdots,X_p)\)​ 服从 \(p\)​ 元正态分布 \(N_p(\mu,\Sigma)\)​ ，这里我们主要讨论参数 \(\mu\)​ 和 \(\Sigma\)​​ 的极大似然估计及其性质。设随机阵 \(X\)​ 表示一组样本容量为 \(n\)​​​​​ 的简单随机样本：

首先从样本数据阵 \(X\)​ 出发，可以定义如下相关的统计量。

样本均值向量，即对 \(X\) 的每个分量求样本均值，得到的一个 \(p\)​ 维向量：

其中，\(\bar{x}_j\)​ 表示第 \(j\) 个分量 \(X_j\)​ 的样本均值：

样本离差阵（交叉乘积阵），类比于一元总体的简单随机样本的离差平方和：

在已知样本数据阵的情况下，常用最后一个表达式计算样本离差阵。由样本离差阵的定义，易知 \(A\) 是一个 \(p\times p\)​​ 的对称矩阵，且有

样本协方差阵，其定义类似于样本方差，由样本离差阵除以自由度可得：

所以 \(S\) 也是一个 \(p\times p\) 的对称矩阵，其对角线元素 \(s_{jj}\)​ 的表达式为：

易知 \(s_{jj}\)​ 表示分量 \(X_j\)​ 的样本方差，其平方根 \(\sqrt{s_{jj}}\)​ 表示分量 \(X_j\)​​ 的样本标准差。此外 \(S\) 的非对角线元素 \(s_{ij}\ (i\neq j)\) 表示分量 \(X_i\) 和 \(X_j\)​​​ 的样本协方差。

有时我们也将样本协方差阵定义为

样本相关阵，其元素由样本相关系数构成，因此用样本协方差阵的元素即可定义：

易知 \(R\) 是一个对角线元素均为 \(1\) 的 \(p\times p\) 的对称矩阵。

Part 2：似然函数

用极大似然法求参数 \(\mu\) 和 \(\Sigma\) 的极大似然估计量，首先需要写出似然函数。似然函数就是样本 \(X\) 的联合密度函数，只不过这里的每个样本都是 \(p\) 元正态随机向量，也就是 \(n\) 个 \(p\) 元正态密度函数的乘积。

我们可以使用拉直运算，将 \({\rm Vec}(X’)\) 的联合密度函数看成参数 \(\mu\) 和 \(\Sigma\)​​ 的函数，就得到了我们所需要的似然函数，记为 \(L(\mu,\Sigma)\) ：

由此求得对数似然函数 \(l(\mu,\Sigma)\) 为：

由上式最后的部分是一个实数，所以可以利用矩阵的迹的有关性质进行变换：

于是我们可以将对数似然函数写为

Part 3：极大似然估计

求解极大似然估计，需要最大化似然函数。一种方法是我们可以对向量 \(\mu\) 和矩阵 \(\Sigma\)​​ 求导，但矩阵微商的计算比较麻烦，所以这里我们介绍一个引理。

引理：设 \(B\) 是 \(p\) 阶正定矩阵，则有 \({\rm tr}B-\ln\left|B\right|\geq p\) ，且等号成立当且仅当 \(B=I_p\) 。

由于 \(B\) 正定，所以 \(B\) 的全部特征值 \(\lambda_1,\lambda_2,\cdots,\lambda_p>0\) ，且 \(\left|B\right|=\lambda_1\lambda_2\cdots\lambda_p\) 。

利用不等式 \(\ln(1+x)\leq x\) 可得

\[\begin{aligned}
\ln|B|&=\sum_{j=1}^p\ln\lambda_j=\sum_{j=1}^p\ln(1+\lambda_j-1) \leq\sum_{j=1}^p(\lambda_j-1)={\rm tr} B-p \ .
\end{aligned}
\]

所以

\[{\rm tr} B-\ln |B|\geq p \ .
\]

由于不等式 \(\ln(1+x)\leq x\)​ 的等号成立条件是 \(x=0\)​ ，所以当且仅当 \(\lambda_1=\lambda_2=\cdots=\lambda_p=1\)​​​ 时上式等号成立，即 \(B=I_p\) 。

首先固定 \(\Sigma>0\) ，由二次型的性质知

以上不等式当且仅当 \(\mu=\bar{X}\) 时等号成立。

进一步取 \(B=\Sigma^{-1/2}\dfrac An\Sigma^{-1/2}\)​​ 正定，利用引理可得

以上不等式当且仅当 \(\Sigma=\dfrac An\)​​​​​​​​​ 时等号成立。

注意这里的第四个等号，只有在矩阵的迹运算和行列式运算中才成立，其原理是

这里用 \(\det(\cdot)\) 表示行列式运算。对于矩阵运算不具有这一性质，即

由以上的推导过程可知参数 \(\mu\) 和 \(\Sigma\) 的极大似然估计量为

似然函数的最大值为

三、参数估计的性质

Part 1：基本统计量的性质

定理：设 \(\bar{X}\) 和 \(A\) 分别为 \(p\) 元正态总体 \(N_p(\mu,\Sigma)\) 的样本均值向量和样本离差阵，样本容量为 \(n\) ，则

(1) \(\bar{X}\sim N_p\left(\mu,\dfrac1n\Sigma\right)\) ；

(2) \(A\xlongequal{d}\displaystyle\sum_{k=1}^{n-1}Z_kZ_k’\) ，其中 \(Z_1,Z_2,\cdots,Z_{n-1}\) 独立同 \(N_p(0,\Sigma)\) 分布；

(3) \(\bar{X}\) 和 \(A\) 相互独立；

(4) \(P(A>0)=1\ \iff\ n>p\) ，即 \(A\) 以概率 \(1\) 正定当且仅当 \(n>p\) 。

该定理的证明和数理统计中一元正态分布的抽样分布类似，需要构造一个正交矩阵，设为 \(\Gamma\) 且具有如下形式

\[\Gamma=\left[\begin{array}{cccc}
\gamma_{11} & \gamma_{12} & \cdots & \gamma_{1n} \\
\vdots &\vdots & & \vdots \\
\gamma_{(n-1),1} & \gamma_{(n-1),2} & \cdots & \gamma_{(n-1),n} \\
\cfrac1{\sqrt{n}} &\cfrac1{\sqrt{n}} & \cdots & \cfrac1{\sqrt{n}}
\end{array}\right]=(\gamma_{ij})_{n\times n} \ .
\]

对样本数据阵构造正交变换，令

\[{\rm Z}=\left[\begin{array}{c}
Z_1′ \\
Z_2′ \\
\vdots \\
Z_n’ \\
\end{array}\right]=\Gamma\left[\begin{array}{c}
X_{(1)}’ \\
X_{(2)}’ \\
\vdots \\
X_{(n)}’ \\
\end{array}\right]=\Gamma X \ ,
\]

即对任意的 \(k=1,2,\cdots,n\)​ 都有

\[Z_k=\left(X_{(1)},X_{(2)},\cdots,X_{(n)}\right)\left[\begin{array}{c}
\gamma_{k1} \\
\gamma_{k2} \\
\vdots \\
\gamma_{kn} \\
\end{array}\right] \ ,
\]

特别地，当 \(k=n\) 时有

\[Z_n=\frac1{\sqrt{n}}\sum_{i=1}^nX_{(i)} \ .
\]

容易证明 \(Z_k\) 是一个 \(p\)​​ 维正态随机向量，且由正交矩阵的性质知

\[\begin{aligned}
&{\rm E}(Z_k)=\sum_{i=1}^n\gamma_{ki}{\rm E}\left(X_{(i)}\right)=\left\{\begin{array}{ll}
0 \ , & k\neq n \ . \\
\sqrt{n}\mu \ , & k=n \ .
\end{array}\right. \\ \\
&\begin{aligned}
{\rm Cov}(Z_k,Z_l)&={\rm E}\left[\left(Z_k-{\rm E}(Z_k)\right)\left(Z_l-{\rm E}(Z_l)\right)’\right] \\ \\
&=\sum_{i=1}^n\gamma_{ki}\gamma_{li}\Sigma=\left\{\begin{array}{ll}
O \ , & k\neq l \ . \\
\Sigma \ , & k=l \ .
\end{array}\right.
\end{aligned}
\end{aligned}
\]

(1) 由已经证明的性质知

\[Z_n=\frac1{\sqrt{n}}\sum_{i=1}^nX_{(i)}=\sqrt{n}\bar{X}\sim N_p\left(\sqrt{n}\mu,\Sigma\right) \ ,
\]

从而可得

\[\bar{X}=\frac{1}{\sqrt{n}}Z_n\sim N_p\left(\mu,\frac1n\Sigma\right) \ .
\]

(2) 因为

\[\sum_{i=1}^nZ_iZ_i’={\rm Z}'{\rm Z}=X’\Gamma’\Gamma X=X’X=\sum_{i=1}^nX_{(i)}X_{(i)}’ \ ,
\]

于是有

\[\begin{aligned}
\sum_{i=1}^{n-1}Z_iZ_i’&=\sum_{i=1}^nX_{(i)}X_{(i)}’-Z_nZ_n’=\sum_{i=1}^nX_{(i)}X_{(i)}’-n\bar{X}\bar{X}’ \\ \\
&=\sum_{i=1}^n\left(X_{(i)}-\bar{X}\right)\left(X_{(i)}-\bar{X}\right)’=A \ .
\end{aligned}
\]

(3) 因为 \(A\) 是 \(Z_1,Z_2,\cdots,Z_{n-1}\) 的函数，\(\bar{X}\) 是 \(Z_n\) 的函数，而 \(Z_1,Z_2,\cdots,Z_{n-1}\)​ 和 \(Z_n\) 相互独立，故 \(A\) 和 \(\bar{X}\) 也相互独立。

(4) 根据以上证明，我们可以令 \(B=\left(Z_1,Z_2,\cdots,Z_{n-1}\right)\)​ 从而 \(A=BB’\)​ 。

如果 \(A\) 正定，则 \(A\) 的秩为 \(p\)​ ​，从而 \(B\) 的秩也为 \(p\) ，于是 \(n-1\geq p\) ，即 \(n>p\) 。

如果 \(n>p\)​ ，要证 \(A\)​ 以概率 \(1\)​ 正定，只需证 \(B\)​ 的前 \(p\)​ 个分量线性相关的概率为 \(0\)​ 。由于 \(B\)​ 是一个多元正态随机阵，所以 \(B\)​ 的前 \(p\)​ 个分量的任意线性组合服从多元正态分布。

所以对于任意不全为零的 \(\beta_1,\beta_2,\cdots,\beta_p\in\mathbb{R}\)​ ，由连续型随机变量的性质知

\[P\left(\sum_{i=1}^p\beta_iZ_i=0\right)=0
\]

进而在统计意义下 \(B\) 的前 \(p\) 个分量以概率 \(1\) 线性无关，从而 \(A\) 以概率 \(1\) 正定。

Part 2：极大似然估计的性质

无偏性：样本均值向量 \(\bar{X}\)​ 是 \(\mu\)​ 的无偏估计，样本协方差阵 \(S=\dfrac1{n-1}A\)​ 是 \(\Sigma\)​ 的无偏估计，但 \(\Sigma\)​ 的极大似然估计量 \(\hat\Sigma=\dfrac1nA\)​ 不是 \(\Sigma\)​​​ 的无偏估计。

有效性：样本均值向量和样本协方差阵 \((\bar{X},S)\) 是 \((\mu,\Sigma)\) 的一致最小方差无偏估计量，也是 \((\mu,\Sigma)\) 的充分完备统计量。

相合性：当 \(n\to\infty\)​​​ 时 \(\bar{X},\hat\Sigma\)​ 是 \(\mu,\Sigma\)​​​​​ 的强相合估计。利用 \({\rm E}(\bar{X})=\mu\)​​ 和 Kolmogorov 强大数定律可知

由于 \(Z_1,Z_2,\cdots,Z_{n-1}\) 独立同分布服从于 \(N_p(0,\Sigma)\) ，所以 \({\rm E}\left(Z_iZ_i’\right)=\Sigma\) ，再利用 Kolmogorov 强大数定律可知