多元正态分布的极大似然估计

1. 一元正态分布的密度函数

一元正态分布的密度函数表示为：

f ( x ) = 1 ( 2 π ) σ e − ( x − μ ) 2 2 σ 2 f(x) = \frac{1}{\sqrt {(2 \pi)} \sigma} e^{- \frac{(x - \mu)^2}{2 \sigma^2}} f(x)=(2π) σ1e−2σ2(x−μ)2其中， σ > 0 \sigma >0 σ>0。由于 x x x、 μ \mu μ均为一维的数值， ( x − μ ) T (x - \mu)^T (x−μ)T与 ( x − μ ) (x - \mu) (x−μ)是等价的，所以上述密度函数又可以写成如下形式：

f ( x ) = 1 ( 2 π ) 1 2 ( σ 2 ) 1 2 e − 1 2 ( x − μ ) T ( σ 2 ) − 1 ( x − μ ) f(x) = \frac{1}{(2 \pi)^{\frac{1}{2}} (\sigma^2)^\frac{1}{2} } e^{- \frac{1}{2} (x - \mu)^T (\sigma^2)^{- 1} (x - \mu)} f(x)=(2π)21(σ2)211e−21(x−μ)T(σ2)−1(x−μ)将上式推广，就得到多元正态分布的定义。

2. 多元正态分布的密度函数

设 K K K维随机向量 x = [ x 1 . . . x K ] x = \left[ \begin{matrix} x_1 \\ ... \\ x_K \end{matrix} \right] x= x1...xK 的密度函数为：

f μ , Σ ( x ) = f μ , Σ ( x 1 , . . . , x K ) = 1 ( 2 π ) K 2 ⋅ 1 ∣ Σ ∣ 1 2 ⋅ e − 1 2 ( x − μ ) T Σ − 1 ( x − μ ) f_{\mu, \Sigma}(x) = f_{\mu, \Sigma}(x_1,...,x_K) = \frac{1}{(2 \pi)^{\frac{K}{2}}} \cdot \frac{1}{|\Sigma|^{\frac{1}{2}}} \cdot e^{- \frac{1}{2} (x - \mu)^T \Sigma^{-1} (x - \mu)} fμ,Σ(x)=fμ,Σ(x1,...,xK)=(2π)2K1⋅∣Σ∣211⋅e−21(x−μ)TΣ−1(x−μ)其中， K K K表示向量 x x x的维度。均值向量 μ \mu μ是 K K K维向量，协方差矩阵 Σ \Sigma Σ是一个 K × K K \times K K×K的对称正定阵，则称 x x x服从 K K K元正态分布，也称 x x x为 K K K维正态随机向量，简记为： x x x~ N K ( μ , Σ ) N_K (\mu, \Sigma) NK(μ,Σ)。显然当 K = 1 K = 1 K=1时，即为一元正态分布的密度函数。

注意，当 ∣ Σ ∣ = 0 |\Sigma| = 0 ∣Σ∣=0时， Σ − 1 \Sigma^{-1} Σ−1不存在， x x x也就不存在通常意义下的密度函数，然而可以形式的给出一个表达式。有些问题可以利用这一形式对 ∣ Σ ∣ ≠ 0 |\Sigma| \neq 0 ∣Σ∣=0及 ∣ Σ ∣ = 0 |\Sigma| = 0 ∣Σ∣=0的情况给出一个统一的处理。

3. 多元正态分布的极大似然估计

对于 N N N个样本点： { x 1 , . . . , x N } \{ x^1,...,x^N\} {x1,...,xN}，其似然函数为：

L ( μ , Σ ) = f μ , Σ ( x 1 ) f μ , Σ ( x 2 ) . . . f μ , Σ ( x N ) L(\mu, \Sigma) = f_{\mu,\Sigma}(x^1) f_{\mu,\Sigma}(x^2)...f_{\mu,\Sigma}(x^N) L(μ,Σ)=fμ,Σ(x1)fμ,Σ(x2)...fμ,Σ(xN) = ( 2 π ) − K N 2 ⋅ ∣ Σ ∣ − N 2 ⋅ e − 1 2 ∑ n = 1 N ( x n − μ ) T Σ − 1 ( x n − μ ) = (2 \pi)^{- \frac{KN}{2}} \cdot |\Sigma|^{- \frac{N}{2}} \cdot e^{- \frac{1}{2} \sum_{n = 1}^N (x^n - \mu)^T \Sigma^{-1} (x^n - \mu)} =(2π)−2KN⋅∣Σ∣−2N⋅e−21∑n=1N(xn−μ)TΣ−1(xn−μ)

对数似然函数：

ln ⁡ L ( μ , Σ ) = − K N 2 ln ⁡ ( 2 π ) − N 2 ln ⁡ ∣ Σ ∣ − 1 2 ∑ n = 1 N ( x n − μ ) T Σ − 1 ( x n − μ ) \ln L(\mu, \Sigma) = - \frac{KN}{2} \ln (2 \pi) - \frac{N}{2} \ln |\Sigma| - \frac{1}{2} \sum_{n = 1}^N (x^n - \mu)^T \Sigma^{-1} (x^n - \mu) lnL(μ,Σ)=−2KNln(2π)−2Nln∣Σ∣−21n=1∑N(xn−μ)TΣ−1(xn−μ) = C − N 2 ln ⁡ ∣ Σ ∣ − 1 2 ∑ n = 1 N ( x n − μ ) T Σ − 1 ( x n − μ ) = C - \frac{N}{2} \ln |\Sigma| - \frac{1}{2} \sum_{n = 1}^N (x^n - \mu)^T \Sigma^{-1} (x^n - \mu) =C−2Nln∣Σ∣−21n=1∑N(xn−μ)TΣ−1(xn−μ)

其中， C = − K N 2 ln ⁡ ( 2 π ) C = - \frac{KN}{2} \ln (2 \pi) C=−2KNln(2π)为一个常数。

【矩阵代数】

一元微积分中，导数（标量对标量的导数）与微分有联系： d f = f ′ ( x ) d x df = f^{'} (x) dx df=f′(x)dx。

多元微积分中，梯度（标量对向量的导数）也与微分有联系： d f = ∑ i = 1 n ∂ f ∂ x i d x i = ∂ f ∂ x T d x df = \sum_{i= 1}^n \frac{\partial f}{\partial x_i} dx_i = \frac{\partial f}{\partial \boldsymbol{x}}^T d \boldsymbol{x} df=∑i=1n∂xi∂fdxi=∂x∂fTdx。这里第一个等号是全微分公式，第二个等号表达了梯度与微分的联系：全微分 d f df df是 n × 1 n \times 1 n×1梯度向量 ∂ f ∂ x \frac{\partial f}{\partial \boldsymbol{x}} ∂x∂f与 n × 1 n \times 1 n×1微分向量 d x dx dx的内积。

据此，矩阵导数与微分也可建立联系： d f = ∑ i = 1 m ∑ j = 1 n ∂ f ∂ X i j d X i j = tr ( ∂ f ∂ X T d X ) df = \sum_{i = 1}^m \sum_{j = 1}^n \frac{\partial f}{\partial X_{ij}} dX_{ij} = \text{tr} \left (\frac{\partial f}{\partial X}^T dX \right) df=∑i=1m∑j=1n∂Xij∂fdXij=tr(∂X∂fTdX)。其中， tr \text{tr} tr代表迹（trace），是方阵对角线元素之和，满足性质：对尺寸相同的矩阵 A A A、 B B B有 tr ( A T B ) = ∑ i , j A i j B i j \text{tr} (A^T B) = \sum_{i,j}A_{ij}B_{ij} tr(ATB)=∑i,jAijBij，即 tr ( A T B ) \text{tr} (A^T B) tr(ATB)是矩阵 A A A、 B B B的内积。与梯度相似，这里第一个等号是全微分公式，第二个等号表达了矩阵导数矩阵与微分的联系：全微分 d f df df是 m × n m \times n m×n导数 ∂ f ∂ X \frac{\partial f}{\partial X} ∂X∂f与 m × n m \times n m×n微分矩阵 d X dX dX的内积。

矩阵微分的运算法则

加减法： d ( X ± Y ) = d X ± d Y d(X \pm Y) = dX \pm dY d(X±Y)=dX±dY
矩阵乘法： d ( X Y ) = d X Y + X d Y d(XY) = dX Y + X dY d(XY)=dXY+XdY
转置： d ( X T ) = ( d X ) T d(X^T) = (dX)^T d(XT)=(dX)T
迹： d tr ( X ) = tr ( d X ) d \text{tr} (X) = \text{tr} (dX) dtr(X)=tr(dX)
逆： d X − 1 = − X − 1 d X X − 1 dX^{-1} = -X^{-1}dX X^{-1} dX−1=−X−1dXX−1。此式可在 X X − 1 = I XX^{-1} = I XX−1=I两侧求微分来证明
行列式： d ∣ X ∣ = tr ( X ∗ d X ) d|X| = \text{tr}(X^* dX) d∣X∣=tr(X∗dX) ，其中 X ∗ X^* X∗表示 X X X的伴随矩阵，在 X X X可逆时又可以写作 d ∣ X ∣ = ∣ X ∣ tr ( X − 1 d X ) d|X|= |X|\text{tr}(X^{-1}dX) d∣X∣=∣X∣tr(X−1dX)。此式可用Laplace展开来证明，详见张贤达《矩阵分析与应用》第279页

通过矩阵导数与微分的联系： d f = tr ( ∂ f ∂ X T d X ) df = \text{tr} \left (\frac{\partial f}{\partial X}^T dX \right) df=tr(∂X∂fTdX)，在求出左侧的微分 d f df df后，可以利用如下一些迹技巧（trace trick）写成右侧的形式并得到导数：

矩阵求导的运算法则

∂ X T A X ∂ X = ( A + A T ) X \frac{\partial X^T A X}{\partial X} = (A + A^T)X ∂X∂XTAX=(A+AT)X，当 A A A为实对称矩阵时， ∂ X T A X ∂ x = 2 A X \frac{\partial X^T A X}{\partial x} = 2AX ∂x∂XTAX=2AX
当 A A A为实对称矩阵时， ∂ ( X T A X ) ∂ A = X X T \frac{\partial (X^T A X)}{\partial A} = XX^T ∂A∂(XTAX)=XXT， ∂ ln ⁡ ∣ A ∣ ∂ A = A − 1 \frac{\partial \ln |A|}{\partial A} = A^{-1} ∂A∂ln∣A∣=A−1
∂ ( X − 1 ) ∂ t = − X − 1 ∂ X ∂ t X − 1 \frac{\partial (X^{-1}) }{\partial t} =- X^{-1} \frac{\partial X}{\partial t} X^{-1} ∂t∂(X−1)=−X−1∂t∂XX−1

对数似然函数分别对 μ \mu μ、 Σ \Sigma Σ求偏导

由上，对数似然函数：

ln ⁡ L ( μ , Σ ) = C − N 2 ln ⁡ ∣ Σ ∣ − 1 2 ∑ n = 1 N ( x n − μ ) T Σ − 1 ( x n − μ ) \ln L(\mu, \Sigma) = C - \frac{N}{2} \ln |\Sigma| - \frac{1}{2} \sum_{n = 1}^N (x^n - \mu)^T \Sigma^{-1} (x^n - \mu) lnL(μ,Σ)=C−2Nln∣Σ∣−21n=1∑N(xn−μ)TΣ−1(xn−μ)

ln ⁡ L ( μ , Σ ) \ln L(\mu, \Sigma) lnL(μ,Σ)对 μ \mu μ求偏导，并令偏导为 0 0 0，即求解 ∂ [ ∑ n = 1 N ( x n − μ ) T Σ − 1 ( x n − μ ) ] ∂ μ = 0 \frac{\partial [\sum_{n = 1}^N (x^n - \mu)^T \Sigma^{-1} (x^n - \mu)]}{\partial \mu} = 0 ∂μ∂[∑n=1N(xn−μ)TΣ−1(xn−μ)]=0，记为 ∂ l 1 ∂ μ = 0 \frac{\partial l_1}{\partial \mu} = 0 ∂μ∂l1=0

将 l 1 = ∑ n = 1 N ( x n − μ ) T Σ − 1 ( x n − μ ) l_1 = \sum_{n = 1}^N (x^n - \mu)^T \Sigma^{-1} (x^n - \mu) l1=∑n=1N(xn−μ)TΣ−1(xn−μ)展开： ∑ n = 1 N [ ( x n ) T Σ − 1 x n − 2 ( x n ) T Σ − 1 μ + μ T Σ − 1 μ ] \sum_{n = 1}^N [ (x^n)^T \Sigma^{-1} x^n - 2(x^n)^T \Sigma^{-1} \mu + \mu^T \Sigma^{-1} \mu ] n=1∑N[(xn)TΣ−1xn−2(xn)TΣ−1μ+μTΣ−1μ] = ∑ n = 1 N ( x n ) T Σ − 1 x n − 2 ∑ n = 1 N ( x n ) T Σ − 1 μ + N μ T Σ − 1 μ = \sum_{n = 1}^N (x^n)^T \Sigma^{-1} x^n - 2 \sum_{n = 1}^N (x^n)^T \Sigma^{-1} \mu + N \mu^T \Sigma^{-1} \mu =n=1∑N(xn)TΣ−1xn−2n=1∑N(xn)TΣ−1μ+NμTΣ−1μ
对第2项求微分：
d ( − 2 ∑ n = 1 N ( x n ) T Σ − 1 μ ) = − 2 ∑ n = 1 N ( x n ) T Σ − 1 d μ = − 2 ∑ n = 1 N t r ( ( x n ) T Σ − 1 d μ ) d(- 2 \sum_{n = 1}^N (x^n)^T \Sigma^{-1} \mu) = -2 \sum_{n = 1}^N (x^n)^T \Sigma^{-1} d \mu = -2 \sum_{n = 1}^N tr((x^n)^T \Sigma^{-1} d \mu) d(−2∑n=1N(xn)TΣ−1μ)=−2∑n=1N(xn)TΣ−1dμ=−2∑n=1Ntr((xn)TΣ−1dμ)
所以第2项对 μ \mu μ的偏导为： ∂ [ − 2 ∑ n = 1 N ( x n ) T Σ − 1 μ ] ∂ μ = − 2 ∑ n = 1 N ( ( x n ) T Σ − 1 ) T = − 2 ∑ n = 1 N Σ − 1 x n \frac{\partial [- 2 \sum_{n = 1}^N (x^n)^T \Sigma^{-1} \mu]}{\partial \mu} = -2 \sum_{n = 1}^N ((x^n)^T \Sigma^{-1})^T = -2 \sum_{n = 1}^N \Sigma^{-1} x^n ∂μ∂[−2∑n=1N(xn)TΣ−1μ]=−2n=1∑N((xn)TΣ−1)T=−2n=1∑NΣ−1xn
第3项对 μ \mu μ的偏导为： ∂ ( N μ T Σ − 1 μ ) ∂ μ = 2 N Σ − 1 μ \frac{\partial (N \mu^T \Sigma^{-1} \mu)}{\partial \mu} = 2N \Sigma^{-1} \mu ∂μ∂(NμTΣ−1μ)=2NΣ−1μ
l 1 l_1 l1对 μ \mu μ的偏导： ∂ l 1 ∂ μ = − 2 ∑ n = 1 N Σ − 1 x n + 2 N Σ − 1 μ \frac{\partial l_1}{\partial \mu} = -2 \sum_{n = 1}^N \Sigma^{-1} x^n + 2N \Sigma^{-1} \mu ∂μ∂l1=−2n=1∑NΣ−1xn+2NΣ−1μ令其等于 0 0 0，解得极大似然估计为： μ ^ = 1 N ∑ n = 1 N x n = x ‾ \hat \mu = \frac{1}{N} \sum_{n = 1}^N x^n = \overline{x} μ^=N1n=1∑Nxn=x

ln ⁡ L ( μ , Σ ) \ln L(\mu, \Sigma) lnL(μ,Σ)对 Σ \Sigma Σ求偏导，并令偏导为 0 0 0，即求解 ∂ [ ln ⁡ ∣ Σ ∣ + 1 N ∑ n = 1 N ( x n − μ ) T Σ − 1 ( x n − μ ) ] ∂ Σ = 0 \frac{\partial [\ln |\Sigma| + \frac{1}{N} \sum_{n = 1}^N (x^n - \mu)^T \Sigma^{-1} (x^n - \mu)]}{\partial \Sigma} = 0 ∂Σ∂[ln∣Σ∣+N1∑n=1N(xn−μ)TΣ−1(xn−μ)]=0，记为 ∂ l 2 ∂ Σ = 0 \frac{\partial l_2}{\partial \Sigma} = 0 ∂Σ∂l2=0

首先求微分，使用矩阵乘法、行列式、逆等运算法则
第一项： d [ ln ⁡ ∣ Σ ∣ ] = ∣ Σ ∣ − 1 d ∣ Σ ∣ = tr ( Σ − 1 d Σ ) d [\ln |\Sigma|] = |\Sigma|^{-1} d |\Sigma| = \text{tr} (\Sigma^{-1} d \Sigma) d[ln∣Σ∣]=∣Σ∣−1d∣Σ∣=tr(Σ−1dΣ)
第二项：
d [ 1 N ∑ n = 1 N ( x n − μ ) T Σ − 1 ( x n − μ ) ] d [\frac{1}{N} \sum_{n = 1}^N (x^n - \mu)^T \Sigma^{-1} (x^n - \mu)] d[N1∑n=1N(xn−μ)TΣ−1(xn−μ)]
= 1 N ∑ n = 1 N ( x n − μ ) T d Σ − 1 ( x n − μ ) = \frac{1}{N} \sum_{n = 1}^N (x^n - \mu)^T d \Sigma^{-1} (x^n - \mu) =N1∑n=1N(xn−μ)TdΣ−1(xn−μ)
= − 1 N ∑ n = 1 N ( x n − μ ) T Σ − 1 d Σ Σ − 1 ( x n − μ ) = - \frac{1}{N} \sum_{n = 1}^N (x^n - \mu)^T \Sigma^{-1} d \Sigma \Sigma^{-1} (x^n - \mu) =−N1∑n=1N(xn−μ)TΣ−1dΣΣ−1(xn−μ)
套上迹，作交换
第一项不变 tr ( Σ − 1 d Σ ) \text{tr} (\Sigma^{-1} d \Sigma) tr(Σ−1dΣ)
第二项：
tr ( − 1 N ∑ n = 1 N ( x n − μ ) T Σ − 1 d Σ Σ − 1 ( x n − μ ) ) \text{tr} \left (- \frac{1}{N} \sum_{n = 1}^N (x^n - \mu)^T \Sigma^{-1} d \Sigma \Sigma^{-1} (x^n - \mu) \right) tr(−N1∑n=1N(xn−μ)TΣ−1dΣΣ−1(xn−μ))
= − 1 N ∑ n = 1 N tr ( ( x n − μ ) T Σ − 1 d Σ Σ − 1 ( x n − μ ) ) = - \frac{1}{N} \sum_{n = 1}^N \text{tr} ((x^n - \mu)^T \Sigma^{-1} d \Sigma \Sigma^{-1} (x^n - \mu)) =−N1∑n=1Ntr((xn−μ)TΣ−1dΣΣ−1(xn−μ))
= − 1 N ∑ n = 1 N tr ( Σ − 1 ( x n − μ ) ( x n − μ ) T Σ − 1 d Σ ) = - \frac{1}{N} \sum_{n = 1}^N \text{tr} (\Sigma^{-1} (x^n - \mu) (x^n - \mu)^T \Sigma^{-1} d \Sigma) =−N1∑n=1Ntr(Σ−1(xn−μ)(xn−μ)TΣ−1dΣ)
= tr ( − Σ − 1 1 N ∑ n = 1 N ( x n − μ ) ( x n − μ ) T Σ − 1 d Σ ) = \text{tr} \left (- \Sigma^{-1} \frac{1}{N} \sum_{n = 1}^N (x^n - \mu) (x^n - \mu)^T \Sigma^{-1} d \Sigma \right) =tr(−Σ−1N1∑n=1N(xn−μ)(xn−μ)TΣ−1dΣ)
其中，第一个等号先交换了 tr \text{tr} tr与 ∑ \sum ∑，第二个等号将 d Σ d \Sigma dΣ右边式子交换到左边，第三个等号再一次交换 tr \text{tr} tr与 ∑ \sum ∑。
所以， d l 2 = tr ( ( Σ − 1 − Σ − 1 1 N ∑ n = 1 N ( x n − μ ) ( x n − μ ) T Σ − 1 ) d Σ ) d l_2 = \text{tr} \left( (\Sigma^{-1} - \Sigma^{-1} \frac{1}{N} \sum_{n = 1}^N (x^n - \mu) (x^n - \mu)^T \Sigma^{-1}) d \Sigma \right) dl2=tr((Σ−1−Σ−1N1∑n=1N(xn−μ)(xn−μ)TΣ−1)dΣ)，对照导数与微分的联系有： ∂ l 2 ∂ Σ = ( Σ − 1 − Σ − 1 1 N ∑ n = 1 N ( x n − μ ) ( x n − μ ) T Σ − 1 ) T \frac{\partial l_2}{\partial \Sigma} = (\Sigma^{-1} - \Sigma^{-1} \frac{1}{N} \sum_{n = 1}^N (x^n - \mu) (x^n - \mu)^T \Sigma^{-1})^T ∂Σ∂l2=(Σ−1−Σ−1N1n=1∑N(xn−μ)(xn−μ)TΣ−1)T令其等于 0 0 0，解得极大似然估计为： Σ ^ = 1 N ∑ n = 1 N ( x n − μ ) ( x n − μ ) T \hat \Sigma = \frac{1}{N} \sum_{n = 1}^N (x^n - \mu)(x^n - \mu)^T Σ^=N1n=1∑N(xn−μ)(xn−μ)T将 μ ^ = x ‾ \hat \mu = \overline{x} μ^=x代入上式得： Σ ^ = 1 N ∑ n = 1 N ( x n − x ‾ ) ( x n − x ‾ ) T \hat \Sigma = \frac{1}{N} \sum_{n = 1}^N (x^n - \overline{x})(x^n - \overline{x})^T Σ^=N1∑n=1N(xn−x)(xn−x)T。

定义样本离差矩阵（又称为交叉乘积阵） V = ∑ n = 1 N ( x n − x ‾ ) ( x n − x ‾ ) T V = \sum_{n = 1}^N (x^n - \overline{x})(x^n - \overline{x})^T V=∑n=1N(xn−x)(xn−x)T，它是一个 K × K K \times K K×K的矩阵。

定义样本协差阵 S = 1 N V = 1 N ∑ n = 1 N ( x n − x ‾ ) ( x n − x ‾ ) T S = \frac{1}{N} V = \frac{1}{N} \sum_{n = 1}^N (x^n - \overline{x})(x^n - \overline{x})^T S=N1V=N1∑n=1N(xn−x)(xn−x)T，它也是一个 K × K K \times K K×K的矩阵。

计算结果

综上，多维正态分布的极大似然估计为： μ ^ = x ‾ ， Σ ^ = S \hat \mu = \overline{x}，\hat \Sigma = S μ^=x，Σ^=S其中， S = 1 N ∑ n = 1 N ( x n − x ‾ ) ( x n − x ‾ ) T S = \frac{1}{N} \sum_{n = 1}^N (x^n - \overline{x})(x^n - \overline{x})^T S=N1∑n=1N(xn−x)(xn−x)T， N N N为样本个数。

参考资料：
矩阵代数_pdf
矩阵求导术（上）
第一章矩阵代数_pdf
第二章多元正态分布的参数估计_ppt
第二章多元正态分布及参数的估计_pdf
多元正态分布参数的估计和数据的清洁与变换_pdf