下で書いたような離散確率の時系列の他に, 特に社会科学などで, 観測値が ロジスティック正規分布に従っている場合が多くあると思います。 つまり具体的には, 観測値 y ∈ {1,0} (1:生起, 0:非生起)だとして,

y 〜 Bernoulli(σ(x)) = Bernoulli (1 / (1 + exp(-x)))
x 〜 N(0,σ²)

になっているようなモデル。回帰モデルでは x がさらに w^Txと回帰になっている場合を考えますが, 議論は基本的に同じです。

これは多項分布の場合はいわゆる対数線形モデルで, 自然言語処理では通常 gradientを計算してL-BFGSやSGDなどの最適化で解くことが多いと思います。 ただし, 最適化の前提となる共変量xが既知ではなく, 学習途中に決まる 潜在変数だったりすると, 最適化してしまうと最初に変な局所解にトラップされて しまい, 学習がうまく動かなくなります。

このため, 例えばBleiの Dynamic Topic Modelではロジスティック正規分布を 下から変分近似して解く方法が示されています。 僕もなぜか長い間, ロジスティック正規分布はこうやって近似するしか ベイズ推定する方法はないのかと思っていたのですが, 最近研究で必要があって 調べたところ, 補助変数を使ってまったく問題なく正しくベイズ推定できることを 知りました。以下, すでに知っているという方は飛ばして下さい。
論文は Groenewald&Mokgatlhe (2005) "Bayesian computation for logistic regression"です。 実は Paul Damien, Jon Wakefield, Stephen Walker(1999), "Gibbs sampling for Bayesian non-conjugate and hierarchical models by using auxiliary variables", JRSS B. の方が著者は有名で, 本質的に同じ内容が 一部書かれていますが, この論文は最初見たときさっぱり分からず, 最初の論文の方が回帰の場合, 多項分布の場合や順序変量の場合など一般的な場合を丁寧に扱っていてわかりやすい ので, そちらの方がお薦めです。

さて, 基本的な戦略はプロビット回帰の場合の有名な Albert&Chib (1993)と同じで, 「観測値を得るために本来あるはずの変数」を潜在変数にして, 補助変数として サンプリングすること。
以下は x を直接推定する場合を考えていますが, xが w^Txのように回帰になっている場合も同じで, 論文にはむしろこちらが書かれています。

具体的には, ロジスティック回帰で p=σ(x)=1/(1+exp(-x)) のとき, 観測値 y が1だということは, [0,1] の一様乱数 u があって u < p だったということなので,

p(y|x) = (σ(x))^y (1-σ(x))^1-y

の補助変数 u との同時確率は,

p(y,u|x)
= I[u≦σ(x)]^y I[u>σ(x)]^1-y …(1)
= I[x≧log (u/(1-u))]^y I[x<log (u/(1-u))]^1-y …(2)

になります。
よって, (1)式から u に依存する項を取り出せば, p(u|y,x)∝p(u)p(y,u|x) なので

• y = 1 のとき p(u|y,x)∝p(u) I[u≦σ(x)] = Unif(0,σ(x))
• y = 0 のとき p(u|y,x)∝p(u) I[u>σ(x)] = Unif(σ(x),1)

になります。 同様に x の場合は, 上の(2)式から

• y = 1のとき p(x|y,u) ∝ p(x) I[x≧log (u/(1-u))]
• y = 0のとき p(x|y,u) ∝ p(x) I[x< log (u/(1-u))]

になります。 これから, u の分布は σ(x) を端点の片方とした一様分布, x の事後分布は事前分布p(x)が正規分布N(0,1)だったとしたとき, 正規分布をyに応じて右側あるいは左側を切り落とした truncated Gaussian になる, ということがわかります。
一般には y は1つではなく, Y = [1 1 0 1] のように複数の観測値があるのが普通 なので, この場合は上式の掛け算を考えると, truncated Gaussian は

p(x|Y,U) ∝ p(x)I(a<x<b),
a = max_i log u_i/(1-u_i)  for y_i=1
b = min_i log u_i/(1-u_i)  for y_i=0

になります。なお, 論文ではこれが掛け算(Π)ではなく, 突然Σで書かれていて, 恐らく間違っているので注意が必要です。Σにしてしまうと, 上式のmaxやminが 導けないはずです。

というわけで, MATLABでプログラムを書いて確かめてみました。 スクリプトはこちらです。 blogis.m
なお, truncated Gaussianからの生成は, 単純に正規分布から生成してから, 範囲に入るまでrejectすればokです。原理的には, 幅が小さいときは効率が悪くなる ので, Chopin (2012)ではテーブルを使って高速化する方法が示されているようです。

y=[1]	y=[0 0 1 0]	y=[1 1 1 1 1 1 1 1 1 1 1 1 1]

下のyが観測値です。 これからわかるように, 綺麗にガウス分布状となる事後分布が推定できています。 y=111111111111のような分布は横が歪んでいるので, 単純な正規分布ではないことが わかります。
上で書いたように, 回帰係数 w を推定する場合も方法は同じなので(論文参照), これでロジスティック回帰も問題なくベイズ推定することができるので安心しました。

関東の T-PRIMAL に対抗して?先日立ち上がった関西の統計的機械学習勉強会 サークルk の帰り道で, *1 T-PRIMALと共催で学生向けの機械学習サマースクールができれば いいなというアイデアを思いついたので, 案をスライドにしてちょっと まとめてみました。

• mlssj.pdf

海外の MLSS (Machine Learning Summer School) の日本版という感じです。
今年はもう夏が近いので, できても2009年度以降だと思いますが, 関東の人が見ても 面白いかも知れないので, 公開してみました。

---

たまにはいわゆる日記を書いてみるテスト。(連休だし。)
普段はこういう日記は個人用の Changelog 日誌に書いているのですが, たまにはいいかなと。

有休を使って, 淡路島の 淡路夢舞台 に行ってきました。淡路花博の時には, 行こうと思いながら天邪鬼なため結局 行かなかったので。
関西からだと, JR舞子駅からJRバスの高速舞子バス停→高速バスで30分で 明石海峡大橋を渡って夢舞台に着きます。帰りはそうしましたが, 行きはそれでは 面白くないので, けいはんなプラザ(9:50発)から関西国際空港(11:30)→ 洲本パールライン で淡路島・洲本まで海路(11:50発)。

淡路島に上陸すると, まず花でいっぱいなのにびっくり。花博が終わった後も, 花の美しいまち, 淡路島という感じでキャンペーンを張ってさまざまな花を町中に植え ているようだ。
洲本港のあたりは昔鐘淵紡績(カネボウ)の工場があったらしく, 工場跡の赤煉瓦の 建物が新開発され, ごちそう館 御食国(みけつくに) や, 洲本アルチザンスクエア として新しいシンボル的場所になっているよう。 アルチザンスクエアはかなり綺麗な建物で, カフェ等が入っている他, 工房やギャラリーなどとしても利用されているようで, なかなか素晴らしい取り組みだと感じました。 洲本市酪農協同組合の出しているアイスクリームの店で飲んだ有機栽培コーヒーは超ウマ。 (洲本市酪農協同組合の取り組みについて) 僕はそういう教育を受けたので, 単に観光客というユーザーとしてというより, どういう風に地方の刷新と活性化を図るといいか, という行政側の視点からも とても興味深い。

洲本バスセンターから淡路島縦貫線で夢舞台までは約60分で, Westin 淡路の前に到着。 部屋はテラス付きで, すごい贅沢でした。ウィスラーの Westin よりも贅沢かも。 安藤建築のこの建物も, V字型のどちらからも海(と僕の部屋からは山)が見える設計 になっていて, 素晴らしい建築でした。
翌日は夢舞台を遊覧。 全編安藤忠雄の建築で素晴らしいのですが(これが大きな目的), 気になったのは, コンクリート打ち放しの耐用年数。花博から5年程度で, すでにある程度は雨で劣化が 始まっており, このまま20年など経った後に, よくあるコンクリート打ち放し建築 のように黒くなっていないか, かなり心配です。少なくとも今の美しさは保てない のではないかと。 (東京都庁はそういう懸念から, 壁面は花崗岩打ち込みになっている。) まあコンクリート打ち放しをやめてしまうと, 安藤建築でなくなってしまうかも 知れないので, 難しいところなのかとも思います。
夢舞台は今も花で一杯で, 素晴らしいところです。
最後に海の方へ橋を渡って 淡路交流の翼港 に行って, 海をのぞんでマターリ。本当に素晴らしい場所でした。