ノンパラメトリック密度推定(ヒストグラム密度推定、カーネル密度推定)
【概要】
【目次】
はじめに
確率モデルを使って推論計算をしていると、解析的な計算が難しい(めんどくさい)場面が出てくると思います。
そこで、MCMCや逐次モンテカルロのようなサンプル集合で確率密度を表現してやろうという手法がよく用いられると思います。
サンプル集合で確率分布を表現することで、その確率分布の特徴(平均や最頻値など)を捉えることはできますが、そのままでは再サンプルができません(逐次計算ができない)。サンプル集合から密度関数を推定する必要がでてきます。
そこで今回は、このような既知の確率分布を仮定しないシーンで用いられる密度関数法であるヒストグラム密度推定(HDE; Histogram Density Estimation)とカーネル密度推定(KDE; Kernel Density Estimation)を実装して確認しました。
まぁHDEは実用的じゃないしKDEはScipyなどで実装があるので実用的にはスクラッチで実装なんて必要ないんですが、あくまで練習ということで。
ノンパラメトリックな密度推定
確率密度を推定するためには、大きく分けてパラメトリックな密度推定とノンパラメトリックな密度推定の二つがあります(参考文献[1])。
パラメトリックな密度推定とは、手持ちのデータを既知の確率密関数に当てはめることを差します。例えば、正規分布を仮定して平均と分散を推定するというのは代表的なパラメトリック密度推定の一つです。
一方、ノンパラメトリックな密度推定とは、既知の確率密度関数を仮定せず、手持ちのデータだけから確率密度を推定します。ヒストグラムを書くイメージです。
今回はノンパラメトリックな密度推定について確認して実装してみます。
ヒストグラム密度推定(Histogram Density Estimation)
最も単純な手法で、データのヒストグラムに基づきます。
手元にN個のデータ()あるとします。そして、確率密度を推定する区間をM区間に分割し、分割幅をそれぞれとします。データ区間iに含まれるデータの個数をとします。 なお、の幅は一定であることは前提にしていません。密度の低い領域はを広くとるみたいなことは有効そうですよね。
このとき、i番目の区間の確率密度は次の通りです。
ヒストグラムそのものです。pythonだとmatplotlibなどでのhist関数におけるdensityオプションなどでおなじみです。
これはとても簡単なのですが、ヒストグラムのbinの幅によって得られる密度関数が大きく変わってきます。例えば[0,1)の範囲でN=100のデータの場合、bin幅を0.1にして密度推定した場合以下のような結果になりました。(実装については次章)
それなりに真の密度関数を推定できていそうですが、だいぶ粗いですね。ということで、bin幅を0.01にしてみます。
これだと、binが狭すぎて推定結果はノイジーですね。
ということで、HDEは単純でしたが、実用するには大量のデータとパラメータ調整が必要になります。また、密度が計算できる点が瓶の幅に依存してしまうので、そういった点でも実用的ではなさそうですね。
カーネル密度推定(Kernel Density Estimation)
まず、任意の位置xでの確率密度をとして、xを中心にした微小領域Vを考えます。このVの領域に含まれるデータの個数をKとします。すると、は次のようになります。
ただしこれは、Vが小さく、Kが十分に含まれるという仮定に基づくので、データが少なかったりVの指定で推定結果が大きく変わってしまうのはHDEと同じです。
HDEと大きく違うのは、予め区間を切るのではなく、xを連続的に動かし、その周辺の密度を逐次計算していくところかなと思います(イメージ図を載せましたが、全然伝わらなさそうな図になってしまった)。
このKを計算するために、次の関数を定義します。
これは、Parzen窓と呼ばれる関数です(参考文献[1])。
この関数を使うと密度関数は次のように書くことができます。
xを中心として幅hの領域にデータがあれば1を無ければ0を返すので、これで各がxを中心にした領域に含まれるかがわかります。これを総和することでKが計算できるわけです。 なお、データの次元は1次元とするので、領域Vは幅hとなります。
ここで、上記の密度推定の式のΣの中は、ターゲット位置xについて各の影響度を総和しているとみなせます。この場合、近いか遠いかを閾値hを使って不連続に切っています。そのため、後で実装結果を示しますが、推定結果はギザギザとなる傾向があります。
の影響力を連続的に伝える仕組みとして、次のようにガウスカーネルを利用するという方法が取れます。
先のHDEで使ったデータを利用して密度推定しみました。hは0.05としています。
このように、KDE(Parzen)は推定結果がバタついていますが、KDE(Gauss)は滑らかですね。hを調整していくとKDE(Parzen)も安定していくポイントはありますが、ガウスカーネルを使った方がだいたいどの場面でも滑らかな推定結果が得られるようです。
実装
具体的な実装を添付します。計算効率などあまり深く考えていないので、あくまで参考程度にしかならないですが。
HDEのクラスはこのように書きました。
KDEのクラス。
おわりに
ということで今回は、ノンパラメトリックな確率密度推定法であるヒストグラム密度推定(Histogram Density Estimation)とカーネル密度推定(Kernel Density Estimation)を実装してみました。
実装はしたんですが、あくまでも中身の理解ということで、実用的にはscikit-learnやscipyなどを利用するのが良いと思います。今回書いたのは1次元にしか対応していないですし。
参考資料
- 作者:C.M. ビショップ
- 発売日: 2012/04/05
- メディア: 単行本(ソフトカバー)