ページ

2016年8月20日土曜日

【R】 ディリクレ分布

ディリクレ分布 Wikipedia

今晩は。確率分布マニアックスの時間がやってきました。今宵はディリクレ分布をご紹介します。この分布は、ベイジアン・ネットワークでよく用いられるとのことで、わたくしはこれまで一度も活用したことがない。そこでWikipediaで勉強してみました。



定義:密度関数

密度関数は次のように書ける。(ただし、xi>0、Σxi=1)

\[
f(x_{1},\cdots,x_{K} | \alpha_{1},\cdots,\alpha_{K}) = \frac{1}{B(\alpha)} \prod_{i=1}^{K} x_{i}^{\alpha_{i}-1}
\]

ここでB(α)は多変量ベータ関数。K=2のときを考えると普通のベータ関数になっている。

\[
B(\alpha) = \frac{\displaystyle \prod_{i=1}^{K} \Gamma(\alpha_{i})}{\Gamma \left(\displaystyle \sum_{i=1}^{K} \alpha_{i} \right)}
\]

ディリクレ分布についてもK=2の場合を考えるとベータ分布になっていることが確認できる。ベータ分布はベイズ統計で一変量統計の事前共役分布のデフォルトとして活躍しているので、その多変量版があるというのも自然だ。

ここで多項分布の確率関数との類似性に注目しておこう。多項分布の確率関数は、次のようにかける。(ただし、Σθi=1)

\[
P(x_{1},\cdots,x_{K} | \theta_{1},\cdots,\theta_{K}) = \frac{\left(\displaystyle \sum_{i=1}^{K} x_{i} \right)!}{\displaystyle \prod_{i=1}^{K} x_{i} !} \prod_{i=1}^{K} \theta_{i}^{x_{i}}
\]

見かけがほぼ同じで、違いと言えば一方が実数をとるガンマ関数を使っているのに対して、他方が自然数をとる階乗を使っているということに過ぎない。・・・確率変数xの位置は全然違うけど。

特徴

まずは期待値と共分散について。

\[
E(x_{i}) = \frac{\alpha_{i}}{A}
\]

ただし、
\[
A = \sum_{i=1}^{K} \alpha_{i}
\]

\[
\mbox{Cov}(x_{i},x_{j}) = \left\{ \begin{array}{l l} \displaystyle -\frac{\alpha_{i}\alpha_{j}}{A^{2}(A+1)} & i \neq j \\ \displaystyle \frac{\alpha_{i}(A - \alpha_{i})}{A^{2}(A+1)} & i = j \end{array} \right.
\]

パラメータの解釈について

Bela A. Frigyik; Amol Kapila; Maya R. Gupta (2010)によると、ディリクレ分布は1mの紙テープをK個にちぎる様子に見立てることが出来る。K=3のときに、パラメータαの違いによって右図のようなデータが生成されている。(1,1,1)は特にどの部分にも肩入れせずにいい加減にちぎっている。仕事が雑である。(0.1,0.1,0.1)はそれに輪をかけていい加減である。時に全然切れていないときがある。

それに対して(10,10,10)は職人である(前二者に比べればだが)。ほぼ3つの部分を切り出している。(2,5,15)はムラがある。最後の1片を切り出すことには熱心だ。

なので、パラメータαは「Concentration」と呼ばれている。多変量正規分布でのConcentration行列とは・・・関係があるかどうか不明だが・・・とりえず別物と認識しておこう。

ベイジアン・ネットワークへの応用。特にLDA(Latent Diriclet Allocation)と呼ばれる技術については稿を改めて紹介することにしよう。


追記:さらにマニアの世界に

https://www.amazon.co.jp/Continuous-Multivariate-Distributions-Applications-Probability/dp/0471183873/ref=sr_1_9?ie=UTF8&qid=1471693226&sr=8-9&keywords=Multivariate+Distribution


右のような本があります。著者のKotzさんたちはたぶん変態です。この種の本(1冊が小さい辞典ほどの厚さがある)を仲間たちとともに、5,6冊作り上げています。この中からディリクレ分布に関するマニアな記述を紹介していきます。

といっても、まずは基本から。ディリクレ分布はごく当たり前のカイ自乗分布から構成される。

\[
Z_{0}, Z_{1},\cdots,Z_{j},\cdots,Z_{m} \sim \chi^{2}(v_{j})
\]
ただし(自由度が整数でなくても良い)、

\[
v_{j} (\in R) \geq 0
\]


このとき次の量、Yjの同時確率分布を考える。

\[
Y_{j} = \frac{Z_{j}}{\displaystyle \sum_{i=0}^{m} Z_{i}} \; (j=1,2,\cdots,m)
\]

(j=0スタートではないのは気にかかると思うが、比率を計算したいということを思い出す。Y0は他の変数から求められるから、わざわざ書かなくていい。)むしろY0は元の情報を覚えておくために次のように定義しておく。

\[
Y_{0} = \sum_{i=0}^{m} Z_{i}
\]


その密度関数は(ヤコビアンを使った合成式を応用する)、
\[
P (y_{0},y_{1},\cdots,y_{m}) = \left[ 2^{(1/2) \sum_{j=0}^{m} v_{j}} \prod_{j=0}^{m} \Gamma \left( \frac{1}{2} v_{j} \right) \right]^{-1} \left[ \prod_{j=1}^{m} y_{j}^{(v_{j}/2)-1} \right]
\]
\[
\times y_{0}^{(1/2) \sum_{j=0}^{m} (v_{j}-1)} \left( 1-\sum_{j=1}^{m}y_{j} \right)^{(v_{0}/2)-1} e^{-(y_{0}/2)}
\]

これをy0について積分(integrating out)して、

\[
P(y_{1},\cdots,y_{m}) = \frac{\displaystyle \Gamma \left( \frac{1}{2} \sum_{j=0}^{m} v_{j} \right)}{\displaystyle \prod_{j=0}^{m} \Gamma \left( \frac{1}{2} v_{j} \right)} \left[ \prod_{j=1}^{m} y_{j}^{(v_{j}/2)-1} \right] \left( 1-\sum_{j=1}^{m}y_{j} \right)^{(v_{0}/2)-1}
\]

この式はディリクレさんが1839年に計算した積分式の特殊例なので、それを記念してこの分布の名前が「ディリクレ分布」と(1946年にクラメールさんによって)命名されたとのことです。

ここで、αj=vj/2と置き換えて、標準型とする。

\[
P(y_{1},\cdots,y_{m}) = \frac{\displaystyle \Gamma \left( \sum_{j=0}^{m} \alpha_{j} \right)}{\displaystyle \prod_{j=0}^{m} \Gamma \left( \alpha_{j} \right)} \left[ \prod_{j=1}^{m} y_{j}^{\alpha_{j}-1} \right] \left( 1-\sum_{j=1}^{m}y_{j} \right)^{\alpha_{0}-1}
\]
ただし(Yjは比率であることを思い出す)、
\[
y_{1},\cdots,y_{m} \geq 0, \sum_{j=1}^{m} y_{j} \leq 1
\]

(冒頭の密度関数と形が違うんでねえの、と思う人もあろうと思う。むしろこちらの方が正しくて、冒頭の式は1-Σyを便宜的にxKと置いて表現を整えたものに過ぎない。)

(平均と共分散の導出)

一般的なモーメントを計算してみる。

\[
\mu_{t_{1},\cdots,t_{m}} = E \left[ \prod_{j=1}^{m} Y_{j}^{t_{j}} \right]
\]
\[
= \frac{\displaystyle \Gamma \left( \sum_{j=0}^{m} \alpha_{j} \right)}{\displaystyle \prod_{j=0}^{m} \Gamma \left( \alpha_{j} \right)} \int \int \cdots \int \left( 1-\sum_{j=1}^{m}y_{j} \right)^{\alpha_{0}-1} \left[ \prod_{j=1}^{m} y_{j}^{\alpha_{j}+t_{j}-1} \right] d \boldsymbol{y}
\]

ところで、密度関数を積分すると1になることを利用して、

\[
\frac{\displaystyle \Gamma \left( \sum_{j=0}^{m} \alpha_{j} \right)}{\displaystyle \prod_{j=0}^{m} \Gamma \left( \alpha_{j} \right)} \int \int \cdots \int \left[ \prod_{j=1}^{m} y_{j}^{\alpha_{j}-1} \right] \left( 1-\sum_{j=1}^{m}y_{j} \right)^{\alpha_{0}-1} d \boldsymbol{y} = 1
\]

だから、

\[
\int \int \cdots \int \left[ \prod_{j=1}^{m} y_{j}^{\alpha_{j}-1} \right] \left( 1-\sum_{j=1}^{m}y_{j} \right)^{\alpha_{0}-1} d \boldsymbol{y} = \frac{\displaystyle \prod_{j=0}^{m} \Gamma \left( \alpha_{j} \right)}{\displaystyle \Gamma \left( \sum_{j=0}^{m} \alpha_{j} \right)}
\]

αj → αj+tjと置き換えて(t0=0とみなす)、

\[
\int \int \cdots \int \left[ \prod_{j=1}^{m} y_{j}^{\alpha_{j}+t_{j}-1} \right] \left( 1-\sum_{j=1}^{m}y_{j} \right)^{\alpha_{0}-1} d \boldsymbol{y} = \frac{\displaystyle \prod_{j=0}^{m} \Gamma \left( \alpha_{j}+t_{j} \right)}{\displaystyle \Gamma \left( \sum_{j=0}^{m} \alpha_{j}+t_{j} \right)}
\]

これをμの式に代入して、

\[
\mu_{t_{1},\cdots,t_{m}} = \frac{\displaystyle \Gamma \left( \sum_{j=0}^{m} \alpha_{j} \right)}{\displaystyle \prod_{j=0}^{m} \Gamma \left( \alpha_{j} \right)} \frac{\displaystyle \prod_{j=0}^{m} \Gamma \left( \alpha_{j}+t_{j} \right)}{\displaystyle \Gamma \left( \sum_{j=0}^{m} \alpha_{j}+t_{j} \right)}
\]

平均は、ある一つの変数についてのみti=1として、他はtj=0とすることによって、
\[
E \left[ Y_{i} \right] = \frac{\alpha_{i}}{\displaystyle \sum_{j=0}^{m} \alpha_{j}} = \frac{\alpha_{i}}{A}
\]


二次モーメントは、変数iについてti=1、変数jについてtj=1として、他はtj=0とすることによって、
\[
\frac{\alpha_{i} \alpha_{j}}{(A+1)A}
\]

これから共分散は、
\[
\mbox{Cov}_{ij} = \frac{\alpha_{i} \alpha_{j}}{(A+1)A} - \frac{\alpha_{i}}{A} \frac{\alpha_{j}}{A} = - \frac{\alpha_{i}\alpha_{j}}{A^{2}(A+1)}
\]

普通の分散は面白いことに、共分散とは異なる式になる。ガンマ関数を二回外すかどうかで違いが出る。
\[
\mbox{Var}_{i} = \frac{\alpha_{i} (A-\alpha_{i})}{A^{2}(A+1)}
\]

相関係数は、

\[
r_{ij} = \frac{\mbox{Cov}_{ij}}{\sqrt{\mbox{Var}_{i}}\sqrt{\mbox{Var}_{j}}} = - \sqrt{\displaystyle \frac{\alpha_{i}\alpha_{j}}{(A-\alpha_{i})(A-\alpha_{j})} }
\]



Rにおけるディリクレ分布
Dirichlet {MCMCpack}
Dirichlet {lca}
Package ‘DirichletReg’
dist.Dirichlet {LaplacesDemon}

ちょっと調べただけで、Rパッケージでディリクレ分布を取り扱うものは上のリストのように見つかる。ここでは最も有名なMCMCpack(Markov Chain Monte Carlo Package)を検討してみよう。


サンプル・コードを動かすとこうなる。

> density <- ddirichlet(c(.1,.2,.7), c(1,1,1))
> density
[1] 2
> 
> draws <- rdirichlet(20, c(1,1,1) )
> draws
            [,1]       [,2]       [,3]
 [1,] 0.31737823 0.20280556 0.47981621
 [2,] 0.05150236 0.47300579 0.47549185
 [3,] 0.56493603 0.26624087 0.16882309
 [4,] 0.20959677 0.14230472 0.64809851
 [5,] 0.08069126 0.74290071 0.17640803
 [6,] 0.31567819 0.26100751 0.42331430
 [7,] 0.26357681 0.39343370 0.34298949
 [8,] 0.39649968 0.31422641 0.28927391
 [9,] 0.25399596 0.72440238 0.02160166
[10,] 0.74225166 0.16009573 0.09765261
[11,] 0.44249988 0.32435725 0.23314287
[12,] 0.16561792 0.39196033 0.44242174
[13,] 0.60624552 0.02430257 0.36945191
[14,] 0.72258069 0.02265145 0.25476786
[15,] 0.04803853 0.20016514 0.75179633
[16,] 0.23066619 0.04585696 0.72347686
[17,] 0.07885911 0.10094257 0.82019831
[18,] 0.57113408 0.03846995 0.39039597
[19,] 0.18952093 0.02232665 0.78815241
[20,] 0.06723889 0.48816993 0.44459118


これだけだとうまく動いているかどうかわからないので、まず密度関数の計算ddirichletにいろいろな数字を入れてみる。

x0 <- c()
for(x00 in seq(.1,.9,length=30)){
 for(y00 in seq(.1,.9,length=30)){
  x0 <- rbind(x0,c(x00,y00,1-x00-y00))
 }
}
density <- ddirichlet(x0, c(6,2,2))
dg <- cbind(x0[,1:2],density)
colnames(dg) <- c("x00","y00","dens.")
plot(dg[,2],dg[,3],type="l",xlab="y",ylab="dens.")
これはα=(6,2,2)について、x,yをそれぞれ動かしてみたときのy側から見たプロット。平均が0.2のはずだが、だいたいそのあたりに山が来ている。 次に生成した乱数から密度推定してみる。
library(KernSmooth)
library(MASS)
draws <- rdirichlet(600, c(6,2,2) )
f1 <- bkde2D(draws[,1:2],bandwidth=c(width.SJ(draws[,1]),width.SJ(draws[,2])))
persp(f1$fhat,phi=30,theta=60,d=5,xlab="x",ylab="y",zlab="dens.")
やはりα=(6,2,2)について、カーネル密度推定した(KernSmooth上のbkde2D関数を使う)。xについて0.6、yについて0.2あたりにピークが来ているので、きっと正しく動いているのだろう。

今回はここまで。