Mar 19, 2018

ディリクレ分布を直観的に理解する

サイコロを振って1,2,3,4,5,6の目が出た回数がそれぞれn1,n2,n3,n4,n5,n6だった時、それぞれの目が出る確率p1,p2,p3,p4,p5,p6は(n1,n2,n3,n4,n5,n6)をパラメーターとするディレクトリ分布に従う、という話を聞いたので、ちょっと勉強しておくことにした。

コインの表が出る確率がpの時、N回振って表が出る回数nはNとpをパラメーターとする二項分布に従う。
コインをN回振って表が出た回数がn回だと観測された時、点推定では表が出る確率はn/Nとなるが、ベイズ推定では表が出る確率(事後確率)pは(n+1, N-n+1)をパラメーターとするベータ分布に従う。
勝手な表記だが、回数がnになる確率P(n)が
P(n;N,p)=Binomial\_PMF(n;N,p)=\pmatrix{N \cr n}p^n(1-p)^{N-n}
(PMF=Probability mass function)という感じであり、事後確率がpになる確率の密度関数f(p)が
f(p;N,n)=Beta\_PDF(p;n+1,N-n+1)=\frac{\Gamma(N+1)}{\Gamma(n+1)\Gamma(N-n+1)}p^n(1-p)^{N-n}=\pmatrix{N \cr n}p^n(1-p)^{N-n}
(PDF: probability density function)ということになる。
この関係のことを、ベイズ推定の用語で、尤度関数が二項分布の場合、ベータ分布は共役事前分布(conjugate prior)であると言うらしい。

多項分布とディリクレ分布もその関係にある。
サイコロをN回振って1〜6の目が出る確率がp1〜p6の時、1〜6の目が出る回数n1〜n6はp1〜p6をパラメーターとする多項分布に従う。(n1+n2+n3+n4+n5+n6=N, p1+p2+p3+p4+p5+p6=1)
サイコロをN回振って1〜6の目が出た回数がそれぞれn1〜n6回だと観測された時、推定される1〜6の目が出る事後確率p1〜p6は(n1+1, n2+1, n3+1, n4+1, n5+1, n6+1)をパラメーターとするディリクレ分布に従う。
nをn1〜n6のベクトル、pをp1〜p6のベクトルとすると、回数がnになる確率P(n)が
P(n;p)=Multinomial\_PMF(n;p)=\frac{\left(\sum_{k=1}^{6}n_k\right)!}{\prod_{k=1}^{6}n_k!}\prod_{k=1}^{6}p_k^{n_k}
であり、事後確率がpになる確率の密度関数f(p)が
f(p;n)=Dirichlet\_PDF(p;n+1)=\frac{\Gamma\left((\sum_{k=1}^6n_k)+1\right)}{\prod_{k=1}^6\Gamma(n_k+1)}\prod_{k=1}^6p_k^{n_k}=\frac{\left(\sum_{k=1}^6n_k\right)!}{\prod_{k=1}^6n_k!}\prod_{k=1}^6p_k^{n_k}
ということである。

ディリクレ分布はベータ分布の確率変数の次元を拡張したものである。
上記のコインの例も多項分布とディリクレ分布で表現でき、コインの表が出る確率がp1、裏が出る確率がp2の時、N回振って表が出る回数n1と裏が出る確率n2は(p1, p2)をパラメーターとする多項分布に従う。
コインをN回振って表が出た回数n1回、裏が出た回数がn2だと観測された時、表が出る事後確率p1、裏が出る事後確率p2は(n1+1, n2+1)をパラメーターとするディリクレ分布に従う。

図1はパラメーターを(3n+1, 2n+1)、n=1〜10としたベータ分布のグラフである。 コインを投げて5回中3回、10回中6回、15回中9回、...、50回中30回、表が出た時の表が出る確率の分布に対応する。
図1
5回中3回でも表が出る確率が0.6である確率が最も高いが、他の確率である確率もそれなりに高いのに対し、回数を増すほど確率が0.6近辺に限定されていく様子がわかる。
最尤推定のような点推定では5回中3回でも50回中30回でも単に0.6であり、その尤もらしさが区別されない。

図2は同様にパラメーターを(3n+1, n+1)、n=1〜10としたベータ分布のグラフである。 コインを投げて4回中3回、8回中6回、...、40回中30回、表が出た時の表が出る確率の分布に対応する。 図2

図3はパラメーターを(2+1, 3+1, 5+1)とした3次元のディリクレ分布のグラフである。 3面しか無いサイコロを10回振って各面が2回、3回、5回出た時の各面が出る確率の分布に相当する。
図3
ちょっとややこしいが、三角形の各頂点がサイコロの各面に対応し、三角形内の各頂点への近さが各面の出る確率に対応し、上方向の高さがその確率の組み合わせになる確率であり、高いほど高温の色になっている。右下の頂点が1つ目の面、右上の頂点が2つ目の面、左の頂点が3つ目の面に対応する。XYZ空間で右下の頂点が(1,0,0)、右上の頂点が(0,1,0)、左の頂点が(0,0,1)にあるとすれば、三角形内の点のX座標、Y座標、Z座標が各面の出る確率である。
上から見ると、図4のようなヒートマップになる。
図4
(0.2, 0.3, 0.5)に対応しそうな所が頂点になっている。

図5はパラメーターを(20+1, 30+1, 50+1)とした3次元のディリクレ分布のグラフである。 3面しか無いサイコロを100回振って各面が20回、30回、50回出た時の各面が出る確率の分布に相当する。
図5
図6は上から見た図である。
図6
それぞれの面が出た割合が同じでも、回数を重ねた方が確率が高い範囲が限定される様子がわかる。

See more ...

Posted at 18:17 in 数学 | WriteBacks (0)
WriteBacks

Feb 01, 2018

(Emacs)Jediで数字に不正な補完候補が出る

筆者はEmacsでPythonのコードを書く時に補完機能としてJedi+CompanyModeを使っているが、筆者の環境では数字に対しておかしな補完が働く。

例えば、

このように、数字を打つと"and"などが補完候補として現れ、Enterキーを押すと

このように、打った数字が消えて"and"になってしまう。

そして慌てて"and"を消してもう一度数字を入力してEnterを押すと、また"and"になるのである。
これは不便である。

CompanyModeでなくAutoCompleteだと同じことにはならないが、やはりおかしな補完が起こることがある。

このように、数字を1文字打つと"if"などが補完候補として現れ、Enterキーを押すと

このように、"if"が足される。

筆者の環境は以下である。
・macOS Sierra バージョン10.12.6
・Emacs 24.5.1 (https://emacsformacosx.comからダウンロード)
 ※同サイトからダウンロードしたEmacs 25.2でも起こる
・jedi 20160425.2156
・company-jedi 20151216.1921

.emacsは次のものだけにしても再現する。

(require 'package)
(package-initialize)

;; for Jedi + company
(require 'company)
(add-hook 'python-mode-hook 'company-mode)
(add-to-list 'company-backends 'company-jedi)

;; for Jedi + auto-complete / company-mode
(require 'jedi-core)
(setq jedi:complete-on-dot t)
(add-hook 'python-mode-hook 'jedi:setup)

Webで調べても、同じ症状の話はほとんど見つからない。
唯一見つけたのは、https://github.com/jorgenschaefer/elpy/issues/1115である。昨年に見つけて時々チェックしているが、今年に入って、筆者と同じ症状の報告が追加された。しかし、私には再現できていないという書き込みや、Jediの問題ではない、Python的に間違った補完ではない(??)、Elpyで回避可能といった書き込みがあり、すぐに修正されそうには見えない。

とりあえずの対策として .emacs に

(setq company-idle-delay 0.5)

を加えて、キー入力から0.5秒間は補完が働かないようにした。

See more ...

Posted at 23:10 in PC一般 | WriteBacks (0)
WriteBacks

Jan 29, 2018

ベイジアンネットワークの例題をやり直す

以前にやったFamily Out Problemの事後確率計算は、ベイジアンネットワークらしくない計算だったことに気付いた。ネットワーク全体をまとめて計算するのでは、ネットワークの構造を無視して確率モデル全体に対してベイズの定理を使ってるだけである。ベイジアンネットワークにしてるからには、新たな情報が追加される度に、ノード毎の事前確率が事後確率に更新されるべきであろう。

問題を再掲する。


Family Out Problem

これのP(fo | lo, hb)を計算する。

この初期状態の事前確率は、単純に上から全パターンの確率を足し合わせて、
(P(fo), P(bp), P(lo), P(do), P(hb)) = (0.15, 0.01, 0.1325, 0.3958, 0.2831)
と求まる。

次に、loが確定した(lo=True, ¬lo=False)として、それに合わないP(fo)を事後確率に更新する。
ベイズの定理より、条件付き確率P(fo | lo)は
P(fo|lo) = \frac{P(lo|fo)P(fo)}{P(lo)} = \frac{P(lo|fo)P(fo)}{P(lo|fo)P(fo) + P(lo|\lnot fo)P(\lnot fo)}
= (0.15 * 0.6) / (0.15 * 0.6 + 0.85 * 0.05) = 0.6792...
であり、これに事後確率P(lo)=1.0を掛けたものがfoの事後確率なので、これを新たなP(fo)とする。
P(do), P(hb)はP(fo)に依存するので、これらも計算し直すと、
(P(fo), P(bp), P(lo), P(do), P(hb)) = (0.6792, 0.01, 1.0, 0.7103, 0.5001)
となる。

次に、hbが確定したとして、それに合わせるべく、関係する事前確率を事後確率に更新する。

まず、hbが直接依存するdoの条件付き確率は、しつこいがベイズの定理より、
P(do|hb) = \frac{P(hb|do)P(do)}{P(hb)} = \frac{P(hb|do)P(do)}{P(hb|do)P(do) + P(hb|\lnot do)P(\lnot do)}
= (0.7 * 0.7103) / (0.7 * 0.7103 + 0.01 * 0.2897) = 0.9942...
であり、hbの事後確率が1.0なので、doの事後確率がこれになる。

次に、doが依存するfoの事後確率は、¬doの事後確率 ≠ 0.0なので、¬doを考慮して P(fo|do)P(do) + P(fo|¬do)P(¬do) と計算することになる。見た目にもかなりしつこいがベイズの定理より、条件付き確率は
P(fo|do) = \frac{P(do|fo)P(fo)}{P(do)} = \frac{P(do|fo,bp)P(fo)P(bp)+P(do|fo,\lnot bp)P(fo)P(\lnot bp)}{P(do|fo,bp)P(fo)P(bp)+P(do|fo,\lnot bp)P(fo)P(\lnot bp)+P(do|\lnot fo,bp)P(\lnot fo)P(bp)+P(do|\lnot fo,\lnot bp)P(\lnot fo)P(\lnot bp)}
= (0.99 * 0.6792 * 0.01 + 0.90 * 0.6792 * 0.99) / (0.99 * 0.6792 * 0.01 + 0.90 * 0.6792 * 0.99 + 0.97 * 0.3208 * 0.01 + 0.3 * 0.3208 * 0.99) = 0.86147...
P(fo|\lnot do) = \frac{P(\lnot do|fo)P(fo)}{P(\lnot do)} = \frac{P(\lnot do|fo,bp)P(fo)P(bp)+P(\lnot do|fo,\lnot bp)P(fo)P(\lnot bp)}{P(\lnot do|fo,bp)P(fo)P(bp)+P(\lnot do|fo,\lnot bp)P(fo)P(\lnot bp)+P(\lnot do|\lnot fo,bp)P(\lnot fo)P(bp)+P(\lnot do|\lnot fo,\lnot bp)P(\lnot fo)P(\lnot bp)}
= (0.01 * 0.6792 * 0.01 + 0.10 * 0.6792 * 0.99) / (0.01 * 0.6792 * 0.01 + 0.10 * 0.6792 * 0.99 + 0.03 * 0.3208 * 0.01 + 0.7 * 0.3208 * 0.99) = 0.23232...
なので、foの事後確率はdo及び¬doの事後確率を用いて、
P(fo|do)P(do) + P(fo|¬do)P(¬do) = 0.86147 * 0.9942 + 0.23232 * (1 - 0.9942) = 0.8578...
と計算できる。

bpの事後確率も同様に計算できるので、ついでに計算してまとめると、
(P(fo), P(bp), P(lo), P(do), P(hb)) = (0.8579, 0.0138, 1.0, 0.9942, 1.0)
となる。

See more ...

Posted at 18:00 in 数学 | WriteBacks (0)
WriteBacks