ノンパラメトリックベイズ(0.5)Dirichlet Process Mixture ModelとInfinite GMM -

Dirichlet process mixture model(DPMM)とinfinite Gaussian mixture model。まだ実装していないので(0.5)。infinite GMMを実装したら、次はinfinite HMMを理解したい

普通、サンプルデータを混合ガウス分布でfittingしようとする時、混合数は、人間が手動で与えるけど、混合数はいくらでも大きくなり得る(つまり∞)ようにしておいて、うまく推定する方法(母集団の分布として、無限混合モデルを仮定するようなものなので、パラメータが原理的には無限個あるけど、パラメトリックモデルではあると思う。混合数を決めないので、ノンパラメトリックらしいけど、しっくりこないネーミングではある)。Dirichlet process mixture model(DPMM)を最初に考えたのが誰かはよく分からない(多分、1974年のAntoniakの論文?)けど、Rasmussenという人が、2000年にハイパーパラメータも込みで推定する方法を考えて、Infinite Gaussian Mixture Modelという名前が付いた。これ以後、色んなInfinite XXXが生まれることになったっぽい。

The Inﬁnite Gaussian Mixture Model
http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.45.9111

多分、DPMMだと、正規分布の混合モデルである必要もないし、理屈上は、Infinite GMMでも、正規分布以外の混合モデルでも使うことはできるのだろうと思う。が、実用上は正規分布以外を考えることは、あんまりない気もする。Rasmussenの論文だと、MCMCでパラメータ・ハイパーパラメータの分布をサンプリングしてるだけのように見えるけど、多くのケースでは、"一番いい"パラメータ・ハイパーパラメータが欲しい気がする(つまり、ベイズ推定じゃなく、単にMAP推定したい)。これについては、MCMCで十分な数のサンプリングしてやって、その中で確率密度が最大になってるものを選べば、近似的に、真の最大確率密度を与える点を得られるというのが常套手法のよう

Rasmussenの論文では、ハイパーパラメータ(事前分布のパラメータをハイパーパラメータと呼ぶらしい)αの事前分布の導出が省略されすぎてて理解不能なのだけど

Hyperparameter estimation in Dirichlet process mixture models
http://citeseerx.ist.psu.edu/viewdoc/citations;jsessionid=A8264D1702C8FD7E1DDAEF8B1F6BD6B9?doi=10.1.1.51.6826

を読むといい気がする。この論文の一番最初に、Antoniakが与えた式として
$P(k|\alpha,n) = c_n(k) n! \alpha^k \frac{\Gamma(\alpha)}{\Gamma(\alpha+n)}$
が出てくる。Antoniakの論文は数学の論文なので全体的に読みたくないが、この式は、Ewensの抽出公式というのと、第一種Stiring数の組み合わせ論的解釈からでる。

第一種Stirling数の組み合わせ数学における意味
https://ja.wikipedia.org/wiki/%E3%82%B9%E3%82%BF%E3%83%BC%E3%83%AA%E3%83%B3%E3%82%B0%E6%95%B0#.E7.B5.84.E3.81.BF.E5.90.88.E3.82.8F.E3.81.9B.E6.95.B0.E5.AD.A6.E3.81.AB.E3.81.8A.E3.81.91.E3.82.8B.E6.84.8F.E5.91.B3

Ewens's sampling formula
https://en.wikipedia.org/wiki/Ewens%27s_sampling_formula

#Ewensの抽出公式は、1972年に、集団遺伝学の研究で導かれたものらしい。

これで、例えば、3つの要素{1,2,3}があったらクラスタリングの仕方は{{1,2,3}},{{1},{2,3}},{{2},{1,3}},{{3},{1,2}},{{1},{2},{3}}の4パターンあって、ハイパーパラメータαを決めると、それぞれの分割の出現確率が決まる。αは集団遺伝学の文脈では、突然変異"率"を表し、大きいほど、沢山のクラスターに分割されやすくなる(α→∞の極限で{{1},{2},{3}}への分割確率は1となる)。

正規分布のパラメータの事前分布としては、正規逆ガンマ分布というものを使う(これが、Dirichlet processの基底分布)。多変数(でfull共分散)の場合は、正規逆Wishart分布。この分布にも、パラメータ(ハイパーパラメータ)が入っており、これも決める必要がある。このへんの事前分布の決め方は、"In general the form of the priors are chosen to have (hopefully) reasonable modelling properties, with an eye to mathematical convenience (through the use of conjugate priors)"とか書いてあって、計算に都合のよさそうな形を選んだということ(?)

実装には、ガンマ分布に従う乱数の生成が必要になるけど、例えば、numpyを使う場合、
numpy.random.gamma
https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.gamma.html
の確率密度関数は、Rasmussenの論文の定義とは異なるものになっているに注意が必要。
$p(x|k,\theta) = x^{k-1} \frac{e^{-x/\theta}}{\theta^k \Gamma(k)}$
がnumpyの確率密度関数で、論文のは
$G(x|k,\theta) = \frac{x^{k/2-1} e^{-k x/2\theta}}{\Gamma(k/2) (2\theta/k)^{k/2}}$
となっている。両者は
$G(x|k,\theta) = p(x|k/2 , 2 \theta/k)$
という関係で結びついている。