多變量常態分布的最大似然估計

本文的閱讀等級:中級

\mathbf{x} 為一 p 維連續型隨機向量。若 \mathbf{x} 服從 (非退化) 多變量常態分布,則機率 (概率) 密度函數完全由 p 維平均數向量 \boldsymbol{\mu}=E[\mathbf{x}]p\times p 階共變異數矩陣 \Sigma=\hbox{cov}[\mathbf{x}] 決定,如下:

\displaystyle  \mathcal{N}(\mathbf{x}\vert\boldsymbol{\mu},\Sigma)=\frac{1}{(2\pi)^{p/2}\vert\Sigma\vert^{1/2}}\exp\left\{-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^T\Sigma^{-1}(\mathbf{x}-\boldsymbol{\mu})\right\}

其中 \vert\Sigma\vert=\det \Sigma>0 (見“共變異數矩陣與常態分布”)。英國統計學家費雪 (Ronald Fisher) 認為機率分布只是一個抽象的數學模型,而我們所蒐集的數據僅能用來估計機率分布的參數。給定一筆取自常態分布的隨機樣本 \mathcal{X}=\{\mathbf{x}_1,\ldots,\mathbf{x}_n\},如何估計模型參數,即平均數向量 \boldsymbol{\mu} 和共變異數矩陣 \Sigma?本文介紹費雪提出的參數估計法,稱為最大似然估計 (maximum likelihood estimation)。根據共變異數矩陣的最大似然估計,我們引進皮爾生 (Pearson) 相關係數,並討論平均數向量的最大似然估計的分布。

 
最大似然估計

假設 p 維樣本 \mathcal{X}=\{\mathbf{x}_1,\ldots,\mathbf{x}_n\} 來自常態分布 \mathcal{N}(\boldsymbol{\mu},\Sigma),且 n>p,其中 \boldsymbol{\mu}\Sigma 是未知的參數。所謂似然函數就是在給定 \boldsymbol{\mu}\Sigma 的情況下,樣本 \mathcal{X} 出現的條件機率密度函數 p(\mathcal{X}|\boldsymbol{\mu},\Sigma)。因為樣本 \mathcal{X} 是隨機選取的,也就是說,\mathbf{x}_1,\ldots,\mathbf{x}_n 是獨立的觀察值,似然函數可表示為

\displaystyle\begin{aligned}  \mathcal{L}(\boldsymbol{\mu},\Sigma|\mathcal{X})&=p(\mathcal{X}|\boldsymbol{\mu},\Sigma)=\prod_{k=1}^n\mathcal{N}(\mathbf{x}_k\vert\boldsymbol{\mu},\Sigma)\\  &=\prod_{k=1}^n\frac{1}{(2\pi)^{p/2}\vert\Sigma\vert^{1/2}}\exp\left\{-\frac{1}{2}(\mathbf{x}_k-\boldsymbol{\mu})^T\Sigma^{-1}(\mathbf{x}_k-\boldsymbol{\mu})\right\}.  \end{aligned}

我們的目標要找出可使 \mathcal{L} 最大化的 \boldsymbol{\mu}\Sigma,稱為最大似然估計。對數是單調遞增函數,因此最大化 \mathcal{L} 等價於最大化 \log\mathcal{L},如下:

\displaystyle\begin{aligned}  \log \mathcal{L}(\boldsymbol{\mu},\Sigma|\mathcal{X})&=\sum_{k=1}^n\log \mathcal{N}(\mathbf{x}_k\vert\boldsymbol{\mu},\Sigma)\\  &=\sum_{k=1}^n\left(-\frac{p}{2}\log(2\pi)-\frac{1}{2}\log\vert\Sigma\vert-\frac{1}{2}(\mathbf{x}_k-\boldsymbol{\mu})^T\Sigma^{-1}(\mathbf{x}_k-\boldsymbol{\mu})\right).  \end{aligned}

運用微分學,求導函數 \frac{\partial \log \mathcal{L}}{\partial\boldsymbol{\mu}}\frac{\partial\log \mathcal{L}}{\partial \Sigma} 並設為零,再解出最佳條件式即可。最大似然估計法看似簡單,但過程中涉及複雜艱澀的代數推導,尤其是引用許多矩陣、跡數 (trace) 和行列式的導數公式。

 
使用純量的向量導數公式 (見“矩陣導數”,SV-8,SV-10),

\displaystyle\begin{aligned}  \frac{\partial \log \mathcal{L}}{\partial\boldsymbol{\mu}}&=\frac{\partial}{\partial\boldsymbol{\mu}}\sum_{k=1}^n\left(-\frac{p}{2}\log(2\pi)-\frac{1}{2}\log\vert\Sigma\vert-\frac{1}{2}(\mathbf{x}_k-\boldsymbol{\mu})^T\Sigma^{-1}(\mathbf{x}_k-\boldsymbol{\mu})\right) \\  &=\sum_{k=1}^n\frac{\partial}{\partial\boldsymbol{\mu}}\left(-\frac{p}{2}\log(2\pi)-\frac{1}{2}\log\vert\Sigma\vert-\frac{1}{2}(\mathbf{x}_k-\boldsymbol{\mu})^T\Sigma^{-1}(\mathbf{x}_k-\boldsymbol{\mu})\right) \\  &=\sum_{k=1}^n\frac{\partial}{\partial\boldsymbol{\mu}}\left(-\frac{1}{2}(\mathbf{x}_k-\boldsymbol{\mu})^T\Sigma^{-1}(\mathbf{x}_k-\boldsymbol{\mu})\right) \\  &=\sum_{k=1}^n\frac{\partial}{\partial\boldsymbol{\mu}}\left(-\frac{1}{2}\boldsymbol{\mu}^T\Sigma^{-1}\boldsymbol{\mu}+\boldsymbol{\mu}^T\Sigma^{-1}\mathbf{x}_k  -\frac{1}{2}\mathbf{x}_k^T\Sigma^{-1}\mathbf{x}_k\right) \\  &=\sum_{k=1}^n\left(-\Sigma^{-1}\boldsymbol{\mu}+\Sigma^{-1}\mathbf{x}_k\right)\\  &=\Sigma^{-1}\sum_{k=1}^n(\mathbf{x}_k-\boldsymbol{\mu}).  \end{aligned}

\frac{\partial \log \mathcal{L}}{\partial\boldsymbol{\mu}}=\mathbf{0},可解出 \boldsymbol{\mu} 的最大似然估計即為樣本平均數向量,記作

\displaystyle  \hat{\boldsymbol{\mu}}=\overline{\mathbf{x}}=\frac{1}{n}\sum_{k=1}^n\mathbf{x}_k

 
接著推導共變異數矩陣 \Sigma 的最大似然估計。利用跡數運算化簡似然函數裡的二次型,如下:

\displaystyle\begin{aligned}  \sum_{k=1}^n(\mathbf{x}_i-\boldsymbol{\mu})^T\Sigma^{-1}(\mathbf{x}_k-\boldsymbol{\mu})  &=\sum_{k=1}^n\text{trace}\left((\mathbf{x}_k-\boldsymbol{\mu})^T\Sigma^{-1}(\mathbf{x}_k-\boldsymbol{\mu})\right)\\  &=\sum_{k=1}^n\text{trace}\left(\Sigma^{-1}(\mathbf{x}_k-\boldsymbol{\mu})(\mathbf{x}_k-\boldsymbol{\mu})^T\right)\\  &=\text{trace}\left(\Sigma^{-1}\sum_{k=1}^n(\mathbf{x}_k-\boldsymbol{\mu})(\mathbf{x}_k-\boldsymbol{\mu})^T\right),  \end{aligned}

上面使用了跡數的循環不變性和線性函數性質 (見“矩陣跡數的性質與應用”)。定義 p\times p 階樣本共變異數矩陣 (見“主成分分析”)

\displaystyle  S=\frac{1}{n-1}\sum_{k=1}^n(\mathbf{x}_k-\overline{\mathbf{x}})(\mathbf{x}_k-\overline{\mathbf{x}})^T

乘開上式即可確認 S(i,j)s_{ij} 為第 i 個變數和第 j 個變數的樣本共變異數 (見“樣本平均數、變異數和共變異數”):

\displaystyle    s_{ij}=\frac{1}{n-1}\sum_{k=1}^n(x_{ki}-\overline{x}_i)(x_{kj}-\overline{x}_j)

其中 x_{ki} 是數據點 \mathbf{x}_k 的第 i 元 (第 i 的變數的觀察值),\overline{x}_i=\frac{1}{n}\sum_{k=1}^nx_{ki} 是第 i 個變數的樣本平均數 (即 \overline{\mathbf{x}} 的第 i 元)。我們已經知道 \overline{\mathbf{x}}\boldsymbol{\mu} 是最大似然估計,因此尋找 \Sigma 使最大化 \mathcal{L}(\overline{\mathbf{x}},\Sigma|\mathcal{X}) 等價於最大化

\displaystyle  \log \mathcal{L}(\overline{\mathbf{x}},\Sigma|\mathcal{X})=-\frac{np}{2}\log(2\pi)-\frac{n}{2}\log\vert\Sigma\vert-\frac{n-1}{2}\text{trace}\left(\Sigma^{-1}S\right)

使用跡數與行列式的導數公式 (見“跡數與行列式的導數”,det-2,tr-18),

\displaystyle  \frac{\partial \log \mathcal{L}(\overline{\mathbf{x}},\Sigma|\mathcal{X})}{\partial \Sigma}=-\frac{n}{2}\Sigma^{-1}+\frac{n-1}{2}\Sigma^{-1}S\Sigma^{-1}

\frac{\partial \log \mathcal{L}(\overline{\mathbf{x}},\Sigma|\mathcal{X})}{\partial \Sigma}=0,可得共變異數矩陣 \Sigma 的最大似然估計

\displaystyle  \hat{\Sigma}=\frac{n-1}{n}S=\frac{1}{n}\sum_{k=1}^n(\mathbf{x}_k-\overline{\mathbf{x}})(\mathbf{x}_k-\overline{\mathbf{x}})^T

 
樣本共變異數矩陣 S 和最大似然估計 \hat{\Sigma} 的計算可簡化如下:

\displaystyle\begin{aligned}  \sum_{k=1}^n(\mathbf{x}_k-\overline{\mathbf{x}})(\mathbf{x}_k-\overline{\mathbf{x}})^T  &=\sum_{k=1}^n\left(\mathbf{x}_k\mathbf{x}_i^T-\mathbf{x}_k\overline{\mathbf{x}}^T-\overline{\mathbf{x}}\mathbf{x}_k^T+\overline{\mathbf{x}}\overline{\mathbf{x}}^T\right)\\  &=\sum_{k=1}^n\mathbf{x}_k\mathbf{x}_k^T-\left(\sum_{k=1}^n\mathbf{x}_k\right)\overline{\mathbf{x}}^T-\overline{\mathbf{x}}\left(\sum_{k=1}^n\mathbf{x}_k^T\right)+\sum_{k=1}^n\overline{\mathbf{x}}\overline{\mathbf{x}}^T\\  &=\sum_{k=1}^n\mathbf{x}_k\mathbf{x}_k^T-n\overline{\mathbf{x}}\overline{\mathbf{x}}^T-n\overline{\mathbf{x}}\overline{\mathbf{x}}^T+n\overline{\mathbf{x}}\overline{\mathbf{x}}^T\\  &=\sum_{k=1}^n\mathbf{x}_k\mathbf{x}_k^T-n\overline{\mathbf{x}}\overline{\mathbf{x}}^T.  \end{aligned}

另一方面,寫出

\displaystyle  \sum_{k=1}^n(\mathbf{x}_k-\overline{\mathbf{x}})(\mathbf{x}_k-\overline{\mathbf{x}})^T=\begin{bmatrix}  \mathbf{x}_1-\overline{\mathbf{x}}&\cdots&\mathbf{x}_n-\overline{\mathbf{x}}  \end{bmatrix}\begin{bmatrix}  (\mathbf{x}_1-\overline{\mathbf{x}})^T\\  \vdots\\  (\mathbf{x}_n-\overline{\mathbf{x}})^T  \end{bmatrix}=AA^T

上面令 A=\begin{bmatrix}  \mathbf{x}_1-\overline{\mathbf{x}}&\cdots&\mathbf{x}_n-\overline{\mathbf{x}}  \end{bmatrix} 為一 p\times n 階矩陣。根據定義,

\displaystyle  (\mathbf{x}_1-\overline{\mathbf{x}})+\cdots+(\mathbf{x}_n-\overline{\mathbf{x}})=\sum_{k=1}^n\mathbf{x}_k-n\overline{\mathbf{x}}=\mathbf{0}

可知 \{\mathbf{x}_1-\mathbf{x},\ldots,\mathbf{x}_n-\overline{\mathbf{x}}\} 是一線性相關集,並推得 (見“利用 Gramian 矩陣證明行秩等於列秩”)

\displaystyle  \text{rank}\hat{\Sigma}=\text{rank}(AA^T)=\text{rank}A<n

因為這個緣故,我們要求 n>p,如此由連續隨機樣本 \mathcal{X} 得到的 p\times p\hat{\Sigma} 為可逆矩陣的機率等於 1

 
皮爾生相關係數

給定隨機向量 \mathbf{x}=(x_1,\ldots,x_p)^T 的共變異數矩陣 \Sigma=[\sigma_{ij}],隨機變數 x_ix_j 的相關係數定義為 (見“共變異數矩陣的性質”)

\displaystyle  \rho_{ij}=\frac{\sigma_{ij}}{\sqrt{\sigma_{ii}}\sqrt{\sigma_{jj}}}

以最大似然估計 \hat{\Sigma}=[\hat{\sigma}_{ij}] 替換 \Sigma,其中

\displaystyle  \hat{\sigma}_{ij}=\frac{n-1}{n}s_{ij}=\frac{1}{n}\sum_{k=1}^n(x_{ki}-\overline{x}_i)(x_{kj}-\overline{x}_j)

可得相關係數 \rho_{ij} 的最大似然估計

\displaystyle\begin{aligned}  \hat{\rho}_{ij}&=\frac{\hat{\sigma}_{ij}}{\sqrt{\hat{\sigma}_{ii}}\sqrt{\hat{\sigma}_{jj}}}\\  &=\frac{\sum_{k=1}^n(x_{ki}-\overline{x}_i)(x_{kj}-\overline{x}_j)}{\sqrt{\sum_{k=1}^n(x_{ki}-\overline{x}_i)^2}\sqrt{\sum_{k=1}^n(x_{kj}-\overline{x}_j)^2}}\\  &=\frac{\sum_{k=1}^nx_{ki}x_{kj}-n\overline{x}_i\overline{x}_j}{\sqrt{\sum_{k=1}^nx_{ki}^2-n\overline{x}_i^2}\sqrt{\sum_{k=1}^nx_{kj}^2-n\overline{x}_j^2}},  \end{aligned}

稱為皮爾生相關係數 (見“相關係數”),經常表示為 r_{ij}。皮爾生相關係數的幾何意義請見“數據矩陣的列與行”。

 
樣本平均數向量的分布

按照費雪的觀點,樣本是從所有可能出現的量測隨機選取,由此得到的估計量也具有隨機性。換句話說,估計量本身就是一個隨機變數。下面證明樣本平均數向量 \overline{\mathbf{x}} 服從常態分佈 \mathcal{N}(\boldsymbol{\mu},(1/n)\Sigma)。為了簡化推導過程,我們計算隨機向量 \overline{\mathbf{x}} 的動差生成函數 (moment generating function)。由於數據點 \mathbf{x}_1,\ldots,\mathbf{x}_n 具有獨立性,且每一 \mathbf{x}_k 服從 \mathcal{N}(\boldsymbol{\mu},\Sigma),其動差生成函數為 (見“多變量常態分布”)

\displaystyle  m_{\mathbf{x}_k}(\mathbf{t})=\exp\left(\mathbf{t}^T\boldsymbol{\mu}+\frac{1}{2}\mathbf{t}^T\Sigma\mathbf{t}\right)

可得

\displaystyle\begin{aligned}  m_{\overline{\mathbf{x}}}(\mathbf{t})&=E\left[\exp\left(\mathbf{t}^T\overline{\mathbf{x}}\right)\right]\\  &=E\left[\exp\left(\mathbf{t}^T\left(\frac{1}{n}\sum_{k=1}^n\mathbf{x}_k\right)\right)\right]\\  &=E\left[\exp\left(\sum_{k=1}^n\frac{1}{n}\mathbf{t}^T\mathbf{x}_k\right)\right]\\  &=E\left[\prod_{k=1}^n\exp\left(\frac{1}{n}\mathbf{t}^T\mathbf{x}_k\right)\right]\\  &=\prod_{k=1}^nE\left[\exp\left(\frac{1}{n}\mathbf{t}^T\mathbf{x}_k\right)\right]\\  &=\prod_{k=1}^nm_{\mathbf{x}_k}\left(\frac{1}{n}\mathbf{t}\right)\\  &=\prod_{k=1}^n\exp\left(\frac{1}{n}\mathbf{t}^T\boldsymbol{\mu}+\frac{1}{2n^2}\mathbf{t}^T\Sigma\mathbf{t}\right)\\  &=\exp\left(\mathbf{t}^T\boldsymbol{\mu}+\frac{1}{2}\mathbf{t}^T\left(\frac{1}{n}\Sigma\right)\mathbf{t}\right),  \end{aligned}

證明 \overline{\mathbf{x}} 亦為常態分布,平均數向量為 E[\overline{\mathbf{x}}]=\boldsymbol{\mu},共變異數矩陣為 \hbox{cov}[\overline{\mathbf{x}}]=(1/n)\Sigma

 
因為選取的樣本具有隨機性,因此單問一個估計量的準確性沒有任何意義,我們應該以估計量的機率分布來評判。設樣本 \mathcal{X} 服從由參數 (向量或矩陣) \boldsymbol{\theta} 決定的機率分布,且 \mathbf{w}=h(\mathcal{X})\boldsymbol{\theta} 的估計量。若 E[\mathbf{w}]=\boldsymbol{\theta},則稱之為無偏 (unbiased) 估計,否則稱為有偏 (biased) 估計。因為 E[\overline{\mathbf{x}}]=\boldsymbol{\mu},樣本平均數向量 \overline{\mathbf{x}}\boldsymbol{\mu} 的一個無偏估計。然而,

\displaystyle\begin{aligned}  E\left[\hat{\Sigma}\right]&=E\left[\frac{1}{n}\left(\sum_{k=1}^n\mathbf{x}_k\mathbf{x}_k^T-n\overline{\mathbf{x}}\overline{\mathbf{x}}^T\right)\right]\\  &=\frac{1}{n}\sum_{k=1}^nE\left[\mathbf{x}_k\mathbf{x}_k^T\right]-E\left[\overline{\mathbf{x}}\overline{\mathbf{x}}^T\right]\\  &=\frac{1}{n}\sum_{k=1}^n\left(\Sigma+\boldsymbol{\mu}\boldsymbol{\mu}^T\right)-\left(\frac{1}{n}\Sigma+\boldsymbol{\mu}\boldsymbol{\mu}^T\right)\\  &=\Sigma+\boldsymbol{\mu}\boldsymbol{\mu}^T-\frac{1}{n}\Sigma-\boldsymbol{\mu}\boldsymbol{\mu}^T\\  &=\frac{n-1}{n}\Sigma,  \end{aligned}

上面使用了共變異數矩陣的等價表達式 (見“共變異數矩陣的性質”)

\displaystyle  \hbox{cov}[\mathbf{x}_k]=\Sigma=E\left[\mathbf{x}_k\mathbf{x}_k^T\right]-\boldsymbol{\mu}\boldsymbol{\mu}^T,~~~k=1,\ldots,n

以及

\displaystyle  \hbox{cov}[\overline{\mathbf{x}}]=\frac{1}{n}\Sigma=E\left[\overline{\mathbf{x}}\overline{\mathbf{x}}^T\right]-\boldsymbol{\mu}\boldsymbol{\mu}^T

最大似然估計 \hat{\Sigma}\Sigma 的一個有偏估計。在實際應用中,我們常以樣本共變異數矩陣 S=\frac{n}{n-1}\hat{\Sigma} 取代最大似然估計 \hat{\Sigma},原因在於 S\Sigma 的一個無偏估計。另外,向量空間觀點也提供了一個採用樣本共變異數矩陣 S 的理由,稱為自由度 (degrees of freedom),詳細討論請見“樣本平均數、變異數和共變異數”。

 
除了前述無偏性,統計學還使用其他準則來評判估計量的良窳。例如,一致性 (consistency) 是指樣本愈大,得到的估計量就愈接近參數的真實值;有效性 (efficiency) 是指二個無偏估計量有較小平方誤差者。最大似然估計具備一致性,且當樣本趨於無窮大時,沒有其他具一致性的估計量比最大似然估計更有效[1]

 
引用來源:
[1] 維基百科:Maximum likelihood

Advertisements
This entry was posted in 機率統計 and tagged , , , , . Bookmark the permalink.

3 則回應給 多變量常態分布的最大似然估計

  1. 徐建平 說:

    請問一下,中間有一段推導我看不太懂"利用跡數運算化簡似然函數裡的二次型",在第一行的推導中trace運算為何可以直接就套用?

  2. ccjou 說:

    二次型 \mathbf{x}^TA\mathbf{x} 是一個純量,可視之為 1\times 1 階矩陣,故可取跡數計算。多變量統計學經常使用這個運算技巧。例如,設 \mathbf{x}=\begin{bmatrix} a\\ b \end{bmatrix},則 a^2+b^2=\mathbf{x}^T\mathbf{x}=\hbox{trace}(\mathbf{x}^T\mathbf{x})=\hbox{trace}(\mathbf{x}\mathbf{x}^T)=\hbox{trace}\begin{bmatrix} a^2&ab\\ ab&b^2 \end{bmatrix}

發表迴響

在下方填入你的資料或按右方圖示以社群網站登入:

WordPress.com Logo

您的留言將使用 WordPress.com 帳號。 登出 / 變更 )

Twitter picture

您的留言將使用 Twitter 帳號。 登出 / 變更 )

Facebook照片

您的留言將使用 Facebook 帳號。 登出 / 變更 )

Google+ photo

您的留言將使用 Google+ 帳號。 登出 / 變更 )

連結到 %s