Introduction to Pólya-Gamma distribution
持橋大地
統計数理研究所 数理・推論研究系
daichi@ism.ac.jp

SVM 2017
2017-9-1 (Fri)NAISTIntroduction to Pólya-Gamma distribution – p.1/35

Overview• 離散分布とロジスティック分布• ロジスティック回帰の完全ベイズ推定• Pólya-Gamma 分布による Gibbs サンプラー• ラプラス変換による証明• 多値ロジスティック回帰, Stick-breaking 過程

Introduction to Pólya-Gamma distribution – p.4/35

Background• 自然言語処理の離散分布

p = (p1 , p2 , · · · , pK )(1)

で, 確率 p1 , · · · , pK を直接モデル化するのは限界◦ 共分散が常に一定 (center ↔ centre, nurse ↔ hospital)!
◦k pk = 1 の制約 (分布を modulate するのが難しい)• LDA: トピック分布 {p(w|k)} を通じて, 「まとめて」共分散を表現• word2vec: 単語を個々にモデル化できる, 統計的背景が弱い

Introduction to Pólya-Gamma distribution – p.5/35

ロジスティック分布
exk,p(k) = x1xx2
Ke + e + ··· + e(x1 , · · · , xK ) ∼ N(0, Σ)(2)• アイテム (単語やトピック) 間の共分散を行列 Σ で表現できる• アイテム k に特徴ベクトル xk = (xk1 , · · · , xkD ) があり,重み β = (β1 , · · · , βD ) で回帰される場合は
Te β xk
p(k|x) = ! βT xk
ke· · · ロジスティック回帰モデル.

(3)

Introduction to Pólya-Gamma distribution – p.6/35

回帰係数 β の学習"!
D#exp
j=1 βj xkj"!
#p(k|x) = !
Dexp
kj=1 βj xkj(4)• Gradient を計算して最適化すればよいのでは?

↓教師あり学習なら (大体)Yes, 教師なし学習では No• 学習途中で, アイテム k の「特徴ベクトル」xk 自体が学習すべき

潜在変数◦ SGD 等で最適化してしまうと, あっという間に局所解にはまる◦ 正しい β の事後分布を求める必要.

Introduction to Pólya-Gamma distribution – p.7/35

例: Continuous-space topic models (NL213)
p(w) exp(θdT w)⃗p(w|d) = !
T⃗p(v)
exp(θvd v)(5)• 潜在層が連続値 (ガウス分布) の場合のニューラルネット型

トピックモデル• 文書 d の潜在ベクトル θd も, 単語ベクトル w⃗ も学習前には未知• θ, w⃗ の要素ごとに単純なランダムウォーク MCMC で学習(NL213).
• θ の「事後分布」を求めたい

Introduction to Pólya-Gamma distribution – p.8/35

ロジスティック回帰のベイズ推定• 最も単純な場合:$x ∼ N(0, 1)
y ∼ Bernoulli(σ(x))(σ(x) =1−x )1+e(6)• いま, 観測値 y = (0, 0, 1, 0) だったとする. x の推定値は?
• Gibbs サンプリング: Albert&Chib (1993) の有名な Probit の場合と同様に, Groenewald&Mokgatlhe (2005) などで提案

Introduction to Pólya-Gamma distribution – p.9/35

ロジスティック回帰のベイズ推定 (2)y = [0 0 1 0] の場合

y = [1 1 1 1] の場合

1500

1500

1000

1000500

5000-4-3-2-1012340-4-3-2-101234• x を点推定するのは, 大きく問題あり!
• 特に, 少量の観測値からの推測 (よくある), 社会科学・医学等の応用の場面

Introduction to Pólya-Gamma distribution – p.10/35

ロジスティック回帰のベイズ推定 (3)• x の事後分布を, サンプルではなく分布で推定できないか?
◦ p(x|y) はガウス分布に似ているが, 正確には正規分布ではない• 補助変数 ω を導入することにより, x をガウス分布として推定できる$Draw ω ∼ Pólya-Gamma (n, x)
Draw x ∼ N(µ(ω), V (ω))(7)

• Pólya-Gamma 分布とは?

Introduction to Pólya-Gamma distribution – p.11/35

論文• Polson&Scott, “Bayesian inference for logistic models using

Pólya-Gamma latent variables”,
arXiv (2011), JASA (2013)
Bayesian inference for logistic models usingPólya-Gamma latent variables
Nicholas G. Polson∗

05.0310v3 [stat.ME] 22 Jul 2013

University of Chicago

James G. Scott†Jesse Windle‡
University of Texas at Austin

First Draft: August 2011
This Draft: July 2013

Abstract
We propose a new data-augmentation strategy for fully Bayesian inference in models
with binomial likelihoods. The approach appeals to a new class of Pólya-Gamma distributions, which are constructed in detail. A variety of examples are presented to show
the versatility of the method, including logistic regression, negative binomial regression,• 激しく難解!!
Introduction to Pólya-Gamma distribution – p.12/35

Pólya-Gamma 分布∞1 %gk
PG(b, c) =,21c222π(k− 2 ) + 4π2

gk ∼ Ga(b, 1) .

(8)

k=1

• 無限個のガンマ乱数 Ga(b, 1) の重みつき和• Pólya の名前は, Pólya frequency function という解析学の概念に由来する (Karlin 1968)◦ 指数乱数の無限和で表現できる分布のクラス
454332211000510PG(1, 0)152005101520PG(4, 0)
Introduction to Pólya-Gamma distribution – p.13/35

ガンマ分布
0.25

10.20.8

0.15

0.6

0.1

0.4

0.05

0.2
000510015Ga(4, 1)

51015Ga(1, 1) = Exp(1)ba a−1 −bx
Ga(x|a, b) =x eΓ(a)(9)& ∞(10)• 注:Γ(a) =

xa−1 e−x dx .

0Introduction to Pólya-Gamma distribution – p.14/35

ガンマ分布 (2)• これから, ガンマ分布のラプラス変換は

E[e−tx ]x∼Ga(a,b) ==& ∞e−tx xa−1 e−bx dx(11)xa−1 e−(b+t)x dx

(12)

0& ∞
0Γ(a)
.
=a(b + t)(13)

1Ga(b, 1) に従うとき,
a& ∞−tx−t· a1 x· xb−1 e−x dx
E[e ] =e&0 ∞b−1 −( at +1)x=dxx e• 特に, 確率変数 X が X ∼

(14)
(15)

0= 'Γ(b)(t b1+ a. !

(16)
Introduction to Pólya-Gamma distribution – p.15/35

二値ロジスティック回帰• 簡単のため, まずは二値 {0, 1} のロジスティック回帰モデルについて考える.
• いま, N 個の観測値があるとする.

y = (y1 , y2 , · · · , yN ) ∈ {0, 1}N(17)• このとき, 潜在変数 x ∈ R からのロジスティック回帰

1exp(y = 1|x) = σ(x) ==−x
1+e
1 + ex(18)で y が生成されたとすると,

p(y|x) =N)
i=1(ex )np(yi |x) =(1+ex )N

ただし n は y の中で 1 が出た回数で, n :=

N%

(19)

I(yi = 1) .

i=1

Introduction to Pólya-Gamma distribution – p.16/35

二値ロジスティック回帰 (2)p(y|x) =N)
i=1(ex )np(yi |x) =(1+ex )N(20)• このとき x ∼ N(0, σ 2 ) とすれば, ベイズの定理からp(x|y) ∝ p(y|x)p(x) = e− 2σ12 x2(ex )n(1+ex )N(21)• x の事後分布が解析的な形にならない!!

↓どうするか?

Introduction to Pólya-Gamma distribution – p.17/35

Hyperbolic cosine 関数
ex +e−xcosh x =
2(22)

80
1.2
160

0.8
40

0.6
0.4

20
0.2
0-5-4-3-2-1012y = cosh(x)3450-5-4-3-2-1012

345y = cosh−1 (x)Introduction to Pólya-Gamma distribution – p.18/35

準備 (1)ex +e−xcosh x =2*+2x111e+1x=e + x =2e2 ex• よって,∴x
1 ex +1cosh=x
22 e2xe2−1 x
=2 x
cosh
2e +1(23)
(24)

(25)
(26)

• これから,x"#aax ) a22xe(eacosh−1= 2a x=
22(e +1)a(ex +1)a(27)が成り立つ.
Introduction to Pólya-Gamma distribution – p.19/35

準備 (2)• よって,xa#"ax
a2e2a (e )−1 x
a=2=2coshxa2(e +1)(ex +1)axae(ex )ax
−b−b=
cosh
··2x(ex +1)b2 e2bx
−b (a− 2b )x=2 e
cosh−b2(28)(29)(30)• これを使えば,(ex )nx−N (n− N−N x)
2coshp(y|x) ==2 e.
(1+ex )N2(31)Introduction to Pólya-Gamma distribution – p.20/35

二値ロジスティック回帰 (3)−Np(y|x) = 2n− Nx−N x()
2ecosh2(32)• 上のままでは Gibbs サンプラーの使える共役形式にならないが,ここで次の魔法のような定理を導入する

Theorem 1 ω が Pólya-Gamma 分布 PG(b, 0) に従うとき,
E[e,−tω]ω∼PG(b,0) = cosh−b,t2(33)txx2=すなわち t =とおけば,2222E[e− x2 ω]ω∼PG(b,0) = cosh−b x2.

(34)

Introduction to Pólya-Gamma distribution – p.21/35

注意2E[e− x2 ω]ω∼PG(b,0) = cosh−b x2

は, ガウス分布の分散 (精度) に関する scale mixture を表している
2− x2 ω2x− 2/ω=eなので, 精度 ω = σ −1 を PG 分布に従って様々に変えたときの, ガウス分布 x ∼ N(0, σ 2 ) の確率密度の期待値• e• 参考: t 分布& ∞N(x|µ, (ηλ)−1 )Ga(η|ν/2, ν/2)dη(35)0= St(x|µ, λ, ν)(36)Introduction to Pólya-Gamma distribution – p.22/35

二値ロジスティック回帰 (4)(34) 式をもとの式に代入すれば,p(y|x) =すなわち,(ex )n−N(1+ex )N=2n− Nx(2 )e& ∞2e− x2 ωPG(ω|b, 0)dω(37)02n− Nx − x2 ω()
2p(y, ω|x) ∝ ee
PG(ω|N, 0)./0Gaussian(38)• ここで, 青い部分で次の定理が成り立つことを利用する.

Theorem 22e− x2 ωPG(ω|b, 0) ∼ PG(ω|b, x)◦ PG 事前分布 PG(ω|b, 0) は, 「尤度」e

事後分布 PG(ω|b, x) になる

2− x2 ω(39)によって,Introduction to Pólya-Gamma distribution – p.23/35

二値ロジスティック回帰 (5)• ゆえに,p(x, ω|y) ∝ p(x)p(y, ω|x)

(40)

Gaussian/0.
2x − x2 ω− 2σ12 x2 (n− N)2∝ee
PG(ω|N, 0)
e./0PG(ω|N, x)(41)• これから, 次の 2 ステップを繰り返す Gibbs サンプラーで, x は正しい分布に収束する. !
(i) x の下で, ω を PG(N, x) から draw(ii) ω の下で, x を (41) 式の Gaussian から draw

Introduction to Pólya-Gamma distribution – p.24/35

定理 1 (PG 分布のラプラス変換) の証明∞1 %gkx ∼ PG(b, 0) =1 22π 2(k−2)k=1& ∞E[e−tx ] =e−tx PG(x|b, 0)dxのとき,(42)

(43)

0

=
==& ∞

e−t!

0∞ & ∞)

xk=1 0∞)*
k=11+= cosh−b1
2π 2,b−1x
k=1 (k− 1 )22"∞e!
− 1+

#

t2
2
2π (k− 1 )
2

Γ(b)
t2
2π 2 (k− 21 )x.
2

xb−1 e−x dx+b(準備より)#x(44)

dx(45)(準備より)(46)(47)
Introduction to Pólya-Gamma distribution – p.25/35

Weierstrass の分解定理+∞ *2)
4
x1+ 2cosh x =π (2k − 1)2(48)k=1x :=,xを代入して,2,+∞ *x)x42
cosh
1+ 2=2π (2k − 1)2k=112∞)
x=1+2 (k − 1 )22π2
k=1

よって,cosh−b,x=2∞)

k=1

1

x
1+2π 2 (k − 21 )22−b(49)

(50)

.

(51)

Introduction to Pólya-Gamma distribution – p.26/35

定理 2 (PG 分布の事後分布) の証明• ガウス分布の分散 (精度) の, PG(b, 0) による scale mixture を考える.
22− c2 x− c2 xee
PG(x|b, 0)
PG(x|b, 0)'c(f = 3∞=c2−b− 2 xcosh
PG(x|b, 0)dx
e20(52)• ラプラス変換すると,

E[e−tx ]x∼f =& ∞0e$%c2− t+ 2 x

PG(x|b, 0)dx' (−b ccosh' (b ccosh 2*4+=
c2t+ 2bcosh
2(53)

2(54)

Introduction to Pólya-Gamma distribution – p.27/35

定理 2 (PG 分布の事後分布) の証明 (2)• Weierstrass の分解定理より,

E[e−tx ]x∼f =⎛k=1=∞)

k=1=⎞bc2 /21+ 2∞2)
2π (k− 12 ) ⎟⎜⎝1∞ *)k=1• ここでX=X +t*(55)⎠c2 /2+t1+ 222π (k− 12 )1c2 /2 2b∞) y+ c21+ y2=2y1 + c /2+t
k=1y222π z + c /2
2π 2 z + c2 /2 + t

X +tX+−1+b·y2y+ c2 +t+*"#1 2z = k−2*+−1t= 1+Xなので,2b(56)(57)

Introduction to Pólya-Gamma distribution – p.28/35

定理 2 (PG 分布の事後分布) の証明 (3)E[e−tx]x∼f ==∞)

k=1∞)

1t1+22π 2 z + c2−b(1 + d−1t)kk=1• ここで2−b(58)#2"
2#c1+dk = 2π 2 k−22(59)"−b(1 + d−1)="k1

#b1dk t + 1(60)1はガンマ分布Ga(b, 1) のラプラス変換だから, 変換を元にdk
戻せば
∞%

gk
f=, gk ∼ Ga(b, 1)
1 2c2
22π (k− 2 ) + 2
k=1(61)∞gk
1 %=, gk ∼ Ga(b, 1) = PG(b, c) . !
21c2
22π(k− 2 ) + 4π2
k=1

Introduction to Pólya-Gamma distribution – p.29/35(62)Pólya-Gamma Gibbs sampling: ω の事後分布• 元に戻って, ロジスティック回帰の学習アルゴリズムを導出する.
• 特徴ベクトル xi = (xi1 , · · · , xid ) をもつアイテムについて, ni 回の試行についての観測値 yi がyi ∼ Binom(ni , σ(xTi β))(63)で得られたとする. (i = 1, · · · , N )• このとき, 前の議論から

p(yi |xi , β) ∝ enx2(y− 2i )x − 2i ωie

PG(ωi |ni , 0)(64)なので, ωi の事後分布は

ωi | β ∼ PG(ni , xTi β) .

(65)

Introduction to Pólya-Gamma distribution – p.30/35

Pólya-Gamma Gibbs sampling: β の事後分布• 同様に β については, データ Y = (y1 , · · · , yN ) 全体に関わる

のでN)

p(yi |xi , β) ∝ ei=1ni
Ti=1 (yi − 2 )xi β"Ne(xTβ)2iωii=1 −2"N(66)• よって, κi = yi −ni /2, zi = κi /ωi とおいて
N%

<N ;"
#T
2%ni(xi β)
Tlog p(yi |xi , β) ∝yi −xi β −
ωi
22i=1i=1*+<N ;%
κi T
ωi
T2(xi β) − 2 xi β−=2
ωi
i=1=>
*+N
2%
ωi
κi
∝−xTi β −2ωi

(67)

(68)

(69)

i=11∝ − (z − Xβ)T Ω(z − Xβ)2(70)Introduction to Pólya-Gamma distribution – p.31/35

Pólya-Gamma Gibbs sampling: β の事後分布 (2)log p(β|X, Ω) ∝ (z −Xβ)T Ω(z −Xβ) + β T Σβ(71)TTTT)β+2β(XΩz)+zΩz= β T (XΩX+
Σ./
0-(72)(β − m)T V (β − m) = β T V β − 2β T (V m) + mT V m

(73)

V

• 一方,と比較して,

V m = X T Ωz(74)∴ m = V −1 X T Ωz(75)β |X, Ω ∼ N (V −1 X T Ωz, V −1 ) .

(76)• よって,Introduction to Pólya-Gamma distribution – p.32/35

多値ロジスティック回帰• xi から出力カテゴリ j が選ばれる確率はT

exp(xTi βj − cij )p(j|xi ) = !K=TTβ − c )xβk1+exp(xieiji jk=1%
cij = logexp(xTi βk )e xi β j(77)k̸=j• 上の確率を θij とおくと,

p(βj |β−j , y) ∝KN ))=j)p(j|xi )I(yi =j) p(ℓ ̸= j|xi )I(yi ̸(78)i=1 k=1=N)

θij iN)θij ii=1=i=1⎛!

I(y =j) ⎜I(y =j)⎝xTi βℓxTi βk⎞ee·"x T βkk̸
=jie
⎟k̸
=j
⎠
! K xT βkik=1 e=j)((1 − θij )wi ) I(yi̸

(79)

(80)

Introduction to Pólya-Gamma distribution – p.33/35

多値ロジスティック回帰 (2)p(βj |β−j , y) ∝N)

I(y =j)θij i=j)(1 − θij ) I(yi̸

(81)

i=1

• したがって尤度項はN)eκij θij e2
θij− 2ωij

PG(ωij |ni , 0)(82)

ωij | βj ∼ PG(ni , θij )βj | Ωj ∼ N(mj (Ω), Vj (Ω)) .

(83)

i=1

• これから, 更新式は$Introduction to Pólya-Gamma distribution – p.34/35

まとめ• ロジスティック回帰モデルの完全なベイズ推定◦ Pólya-Gamma 分布に従う補助変数 ω の下で, パラメータはガウス分布に従う• ロジスティック回帰を, 統計モデルのコンポーネントとして使うことができる◦ 教師なし学習でも OK◦ 豊富な共分散構造・連続値モデル• 多値および Stick-breaking 過程への展開, 様々な modulation

Introduction to Pólya-Gamma distribution – p.35/35