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

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

$\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）
よくよく考えてみたら周辺確率密度関数を考えた方が確率論的には筋が良さそうな気がしましたが、そこを議論しようとすると大学教養レベルを超えてしまうので困りました。

主成分分析の正規分布を使った説明について

古典的なデータ分析、特に次元圧縮の手法として主成分分析というものがあります。
さまざまな流儀がありますが今回はサンプル点から計算した共分散行列の固有値を大きい順にいくつか選び、それらの固有値に対応する固有ベクトルを次元圧縮した空間の基底として採用するものを主成分分析と呼ぶことにします．
主成分分析においてなぜ共分散行列の固有値を考えるのかについて，正規分布を使ってざっくりと説明をしたいと思います．
（本来は別の考え方をするのですが、目を瞑っておきます。今回の議論だと少し問題が起こってしまうので本当はよろしくないのですが……）

与えられた標本の集合を $\mathbf{X} = \left\{ x_1 , x_2 , ... , x_n \right\}$ とします．ただし，各 $x_i$ は$D$次元ユークリッド空間の元とします．
このとき標本が正規分布に従って生成されたと考えることにします．
標本集合から平均 $m$ と共分散行列 $S$ を推定します．

$\displaystyle m = \frac{1}{N} \sum^N_{i=1} x_i$.

$\displaystyle S = \frac{1}{N} \sum^N_{i=1} (x_i - m)(x_i - m)^\top$.

この平均と共分散行列を用いて正規分布の確率密度関数は次のように書くことができます．

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

この正規分布の等確率集合を考えます．条件 $ 0 < c < \frac{1}{(\sqrt{2\pi})^D \sqrt{\det S}}$ を満たす実数 $c$ を一つ固定します．確率密度 $c$ の等位集合を定める方程式は $p(x) = c$， すなわち

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

です．この式を整理すると，

$\displaystyle (x-m)^\top S^{-1}(x-m) = -2\log\left( c (\sqrt{2\pi})^D \sqrt{\det S} \right)$

となります．ここでこの式の右辺を $d$ という文字で表すことにする．定数 $c$ の定義から $d>0$ が成り立ちます．
さらに，$y := x-m$ とすると上の式は次のように書き換えられます．

                      $\displaystyle y^\top S^{-1} y = d$.                             (1) 

式(1)がどのような集合を定めるかを調べるために行列 $S$ の性質を少し調べておきましょう．

命題

共分散行列$S$は実対称行列である．

定義どおりに $S$ の成分を書き下せばすぐに分かる．■

線型代数の結果から行列 $S$ はある直交行列 $P$ によって対角化可能，すなわち $P^\top S P$ は対角行列となり，しかも $S$ の固有値はすべて実数になります．

命題

共分散行列$S$は半正定値行列である．

任意のベクトル $x$ に対して $ x^\top S x \ge 0$ が成り立つことを示す．

定義より

$\displaystyle x^\top S x = \frac{1}{N} \sum^N_{i=1} x^\top (x_i -m) (x_i-m)^\top x$

である．ここで，

$\displaystyle a_i = x^\top (x_i-m)$

と置く．この $a$ は実数になる．

この式の両辺の転置を取って $a_i = (x_i-m)^\top x$ となるから $ x^\top (x_i -m) (x_i-m)^\top x = a_i^2 \ge 0$ となるので，

$\displaystyle x^\top S x = \frac{1}{N} \sum^N_{i=1} a_i^2 \ge 0$

が成り立つ．■

対称性と半正定値性から共分散行列$S$の任意の固有値は非負であることが分かります．

ここで困ったことが起こりました．式(1)には $S$ の逆行列が出てきていますが，半正定値では逆行列の存在が保証できません．仕方がないので $S$ は正定値だと仮定して話を進めることにしましょう．
（追記：$S$のランクが落ちていてもその分だけ低い次元の空間にデータが住んでいることがわかる（ http://lembretesombras.blogspot.jp/2017/05/blog-post_27.html ）ので考える空間を制限してやれば問題ありません．）

$S$の逆行列の対角化について考察しましょう．
まず，$S$を対角化すると

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

となるのでした．ここで $\lambda_i$ を $S$ の固有値として $\Lambda = \left( \lambda_i \right)_{i = 1,...,D}$ と定義します．
上の式の逆行列を計算すると $ \Lambda^{-1} = P^\top S^{-1} P$ となります．したがって，$S^{-1}$ の固有値は $S$ の固有値を使って $1/\lambda_i$ と表すことができます．（ $S$ は正定値だと仮定したので固有値はすべて正なので逆数を取っても大丈夫です．）

式(1)に戻ります．新しい変数 $z$ を導入して $z$ と $y$ の関係を $y = Pz$ と定めます．
このとき，式(1)を書き換えると

$\displaystyle y^\top S^{-1} y = (Pz)^\top S^{-1} (Pz) = z^\top (P^\top S^{-1} P) z = z^\top \Lambda^{-1} z$

となります．
いま，$z$ の第 $i$ 成分を $z_i$ と書くことにすると，式(1)は

$ \displaystyle \frac{1}{\lambda_1} z_1^2 + \frac{1}{\lambda_2}z_2^2 + \cdots + \frac{1}{\lambda_D} z_D^2 = d$

と書きことができます．$d>0$ なのでこの式は楕円体を定めます．
各 $z_i$ 軸の長さは $\sqrt{ \lambda_i }$ に比例します．
したがって，$ \lambda_i $ が大きいほどデータの散らばりが大きいということが分かります．
（この $z_i$ 軸は $\lambda_i$ に対応する $S$ の固有ベクトルで定まります．$y = Pz$ としたので，$z$ の空間を $P$ で回転させてやれば $y$ の空間になります．）

以上のように正規分布を使って説明すると主成分分析が分かりやすいような気もするのですが，共分散行列に正定値性を仮定しないといけないあたり議論として筋が悪いような気もします．