共分散行列が正則でない場合の正規分布の確率密度関数の計算

共分散行列が正則でないと正規分布の確率密度関数の値が計算できないという問題があります．これは正規分布の確率密度関数が次のように定義されるからです：

$\displaystyle p(x) = \frac{1}{(\sqrt{2\pi})^n \sqrt{\det S}}  \exp\left( -\frac{1}{2}(x-m)^\top S^{-1}(x-m) \right)$．

ただし，$m$ は平均ベクトル，$S$ は共分散行列です．

この式からわかるように $S$ の逆行列が必要になります．

一般に共分散行列は半正定値なので値が $0$ となる固有値が存在してランクが落ちる可能性があります．

このような場合にも確率密度関数を計算したいので，$S$ のランクが落ちる場合にどのようなことが起こっているかを次の命題で説明します．

命題

観測したサンプルを $x_1, x_2, ... , x_N \in \mathbb{R}^n$ とする．

サンプルから計算した平均は$0$であると仮定する．

共分散行列を $S = \frac{1}{N}\sum^N_{i=1} x_i x_i^\top$ で定義する．

$V$ を $\mathbb{R}^n$ の部分線型空間として $\dim V = k \le n$ とする．

このとき，「任意の $i$ について $x_i \in V$」と「$\mathrm{rank}\,S = n-k$」は同値である．

（証明）

$S$ の固有値を $\lambda_1, \lambda_2, ... , \lambda_n$ として 

$ \displaystyle \Lambda = \left( \begin{array}{c} \lambda_1 & & \\ & \lambda_2 & \\ & & \ddots \\ & & & \lambda_n \end{array} \right)$

とする．このとき，ある直交行列 $P$ が存在して，

$\displaystyle \Lambda = P^\top S P$

と $S$ を分解することができる．

ここで，$y_i = P^\top x_i$ として，

$\displaystyle \Lambda = P^\top S P = P^\top \left( \frac{1}{N} \sum^N_{i=1} x_i x_i^\top \right) P = \frac{1}{N} \sum^N_{i=1} (P^\top x_i)(P^\top x_i)^\top = \frac{1}{N} \sum^N_{i=1} y_i y_i^\top$

が成り立つ．従って，

$\displaystyle \Lambda = \frac{1}{N} \sum^N_{i=1} y_i y_i^\top$

の対角成分に着目すると，$1 \le p \le n$ について

$\displaystyle \lambda_p = \frac{1}{N} \sum^N_{i=1} y_i^{(p)}y_i^{(p)}$

が成り立つ．ここで，$y_i = ( y_i^{(1)}, \ldots, y_i^{(n)})$ です．

従って，

$\displaystyle \lambda_p = 0 \Leftrightarrow \forall i \, , y_i^{(p)} = 0$

が成り立ち，

$\displaystyle \lambda_p \ne 0 \Leftrightarrow \exists i \, , y_i^{(p)} \ne 0$

が成り立つので，命題の主張が導かれる．■

以下では上の命題の状況に従うことにします．
上の命題によれば $S$ のランクが落ちている場合にはサンプルは $\mathbb{R}^n$ の部分空間 $V$ の上にあることがわかります．少なくともデータを観測する限りでは $V$ の外には正規分布は広がっていないので，この部分の確率密度は $0$ としてしまいましょう．（観測できていないだけで本当は広がっていると考えることもできるので，こう考えてしまって良いのかという疑問は残りますが……）

$S$ の固有値を $\lambda_1, \ldots, \lambda_n$ として，$1 \le p \le k$ について $\lambda_p \ne 0$ として，$k+1 \le q \le n$ について $\lambda_q = 0$ と仮定しましょう．
このとき，固有値 $\lambda_{k+1}, \ldots, \lambda_{n}$ に対応する座標を無視して正規分布の確率密度関数を $y$ について書き直すと，

$\displaystyle p(y) = \frac{1}{(\sqrt{2\pi})^n \sqrt{\lambda_1 \cdots \lambda_k}}\exp\left( -\frac{1}{2}\left( \frac{1}{\lambda_1}y^{(1)}y^{(1)} + \cdots + \frac{1}{\lambda_k}y^{(k)}y^{(k)} \right) \right)$

となります．
以上，多少強引でしたが共分散行列のランクが落ちていても固有空間を使えば確率密度が計算できることがわかりました．

追記（2017/5/28）
よくよく考えてみたら周辺確率密度関数を考えた方が確率論的には筋が良さそうな気がしましたが、そこを議論しようとすると大学教養レベルを超えてしまうので困りました。