カーネル密度推定法 - 機械学習に詳しくなりたいブログ

パターン認識と機械学習（下）のカーネル法を勉強しています。カーネル法とは直接関係ないようですが、カーネル密度推定法（パターン認識と機械学習（上）で紹介されています）を勉強したので記事にします。

カーネル密度推定法とは、サンプルデータから元の確率分布を推定するノンパラメトリック手法です。パラメトリックな手法とは、何らかの確率分布を推定する際に、先にモデルを決めてしまうようなやり方です。例えば正規分布であることを仮定すると、サンプルの標本平均と標本分散は正規分布のパラメータの推定量となります。このような方法では、仮定が間違っていれば正確な推定をすることができません。これに対してノンパラメトリックな手法では、分布の形状を仮定せずに推定します。

ここからが本題。ある未知の確率分布 $p(\mathbf{x})$ からサンプルが得られて、このデータから $p(\mathbf{x})$ を推定したいとします。 $\mathbf{x}$ を含む小さな領域 $R$ を考えると、この領域内でデータが生起する確率 $P$ は、

$\displaystyle P = \int_{R} p(\mathbf{x})d\mathbf{x} \tag{1}$

です。 $p(\mathbf{x})$ から $N$ 個のサンプルが得られたとき、各サンプルが領域 $R$ 内である確率は $P$ ですから、 $R$ 内のデータ個数 $K$ は二項分布

${}_N \mathrm{C} _K P^{K}(1-P)^{N-K} \tag{2}$

に従います。（ $N$ 回の試行で確率 $P$ の事象が $K$ 回発生する二項分布）　二項分布の期待値、分散の公式より、 $E[K] = NP$ 、 $V[K] = NP(1-P)$ です。また、期待値と分散は $E[aX] =aE[X]$ 、 $V[aX] =a^{2}V[X]$ の関係がありますから、

$\begin{eqnarray*} \displaystyle E\left[ \frac{K}{N} \right] &=& P \tag{3} \\ \displaystyle V\left[ \frac{K}{N} \right] &=& \frac{P(1-P)}{N} \tag{4} \end{eqnarray*}$

です。式(4)より、 $N$ が大きくなれば、 $\displaystyle \frac{K}{N}$ の分散が小さくなることがわかります。 $\displaystyle \frac{K}{N}$ の期待値は $P$ ですから、分散が小さいとき、

$K \simeq NP \tag{5}$

と書けます。また、 $R$ の範囲内で $p(\mathbf{x})$ が一定とみなせるほど小さいと仮定すると、 $R$ の体積を $V$ として

$P \simeq p(\mathbf{x})V \tag{6}$

が成り立ちます。この式(6)はちょっと戸惑いましたが、1変数の確率密度関数をイメージして、面積で考えたら理解できました。

さて、式(5)、(6)より、ある領域 $R$ 内での $p(\mathbf{x})$ の推定値は

$\displaystyle p(\mathbf{x}) = \frac{K}{NV} \tag{7}$

となります。ここで、 $V$ は領域 $R$ の体積ですから領域の場所に応じて変化する値ですが、カーネル密度推定法では $V$ を固定して考えます。すると $NV$ は定数になります。 $K$ は $R$ 内のデータの個数でしたので、僕の解釈ですが、式(7)は範囲を十分に小さくしたヒストグラムのようなものと考えると、なんかすっきりしました。単純にヒストグラムの値が確率密度の推定値って言ってるだけですよね。多分。

f:id:opabinia2:20190224002058p:plain 適当な正規分布からサンプリングしたデータをヒストグラムにすると普通は左図のようなものになりますが、範囲を十分に小さくすると右図のようになります。で、これだと範囲を狭くしすぎてデータが存在しない部分は確率密度が不明なので、カーネル関数を使って推定するっていうのがカーネル密度推定法です。じゃあ左図のヒストグラムのままでこれを推定値とすればいいじゃんって話なんですけど、1変数程度の場合はこれでもいいそうですが、確率変数の次元が大きくなると通用しなくなるそうです。*1

さて、肝心な推定方法です。以下のような関数を定義します。

$k(\mathbf{u})= \begin{cases} \displaystyle 1, | u_{i} | \le \frac{1}{2} \\ 0, それ以外 \end{cases} \tag{8}$

確率変数の次元を $D$ としたとき、 $i =1, \cdots,D$ です。式(8)は、 $\mathbf{u}$ が原点を中心とする単位立方体内にあるときのみ1となる関数です。次に $\displaystyle k\left( \frac{\mathbf{x}-\mathbf{x}_n}{h} \right)$ を考えます。ここで $\mathbf{x}_n$ は推定したい確率分布から得られたサンプルデータです。これは $\mathbf{x}$ を中心とする一辺が $h$ の立方体内部に、 $\mathbf{x}_n$ があるときに1となりますので、得られたサンプルデータが、この領域内にある個数 $K$ は

$\displaystyle K=\sum_{n=1}^{N} k \left( \frac{\mathbf{x}-\mathbf{x}_n}{h} \right) \tag{9}$

と書けます。式(9)を式(7)に代入すれば、この領域の体積は $V=h^{D}$ ですから、

$\displaystyle p(\mathbf{x}) =\frac{1}{N} \sum_{n=1}^{N} \frac{1}{h^{D}}k \left( \frac{\mathbf{x}-\mathbf{x}_n}{h} \right) \tag{10}$

となります。この式(10)が推定した確率密度関数になります。実は $k(\mathbf{u})$ がカーネル関数と呼ばれるもので、以下を満たせば任意の関数を使うことができます。

$\begin{eqnarray*} k(\mathbf{u}) &\ge&0 \tag{11} \\ \int k(\mathbf{u}) d\mathbf{u} &=& 1 \tag{12} \end{eqnarray*}$

この条件は式(10)が確率密度関数となるために必要な条件です。負の値があったり、全区間で積分して1にならなければ確率密度関数にはなりません。

さて、一般的にはガウスカーネルが用いられるそうです。ガウスカーネルを用いると、式(10)の確率密度関数は次のようになります。

$\displaystyle p(\mathbf{x}) =\frac{1}{N} \sum_{n=1}^{N} \frac{1}{(2 \pi h^{2})^{\frac{D}{2}}} \exp \left( -\frac{ \| \mathbf{x}-\mathbf{x}_n \|^{2}}{2h^{2}} \right) \tag{13}$