Platinum Blogger

2016年5月24日＠統計モデリング
統計モデリング
第六回 配布資料
文献:
a) C. P. Robert and G. Casella: Monte Carlo Statistical Methods. 2nd ed.,
Springer, 2004.
b) C. P. Robert and G. Casella: Introducing Monte Carlo Methods with R. Springer,
2010.
c) J. Albert: Bayesian Computation with R. 2nd ed. Springer, 2009
b,c は日本語訳 (ともに石田基広, 石田和枝訳, 丸善出版) あり
配布資料の一部は以下からもDLできます.
短縮URL http://tinyurl.com/lxb7kb8
担当：田中冬彦
今後の予定
第六回 数値技法（MCMC）
θ (1) , θ ( 2 ) ,  , θ ( M ) のヒストグラム
π (θ | x) の関数形
posterior distribution
0.4
0.3
distribution
0.2
10
M →∞
0
0.0
0.1
5
Frequency
15
sample histogram
-3
-2
-1
0
1
2
-2
-1
0
1
2
theta
第七回 グループ発表１
x
第八回
線形モデルのベイズ解析
第九回
一般化線形モデルのベイズ解析
出典：東京電力HP
今日の内容
１．モンテカルロ法
２．マルコフ核
３．マルコフ連鎖と定常分布
４．マルコフ連鎖モンテカルロ法
～Metropolis-Hastings Method
５．マルコフ連鎖モンテカルロ法
～Gibbs Sampler
６．WinBUGS
１. モンテカルロ法
Introduction
続・冬のカナタ （以下は架空のものです）
・村上冬樹さんは「泣けるシナリオ大賞」を目指して「続・冬
のカナタ」というシナリオを書き上げました.
・目標としては国民の８割以上が泣くシナリオを目指して
います.
・ 周囲の友人に「続・冬のカナタ」を読んでもらい, 泣いたか
どうかを聞きました.
考えてみよう
続・冬のカナタを読んで泣く人の割合（確率）を q とし設定は事前分布以
外, 前回と同じ. q の事後分布は以下のようになったという
π (q | x = 3) = Cq 3 tanh(1 − q ), 0 ≤ q ≤ 1
ここでCは積分定数.
冬樹さんは q に関する以下の確率を評価したいが, どうすればいいか途方に
暮れている・・・・
1
Pr( q ≥ 0.8 | x = 3) = ∫ π (q | x = 3)dq
0.8
統計の素養を身につけたあなたから冬樹さんへのアドバイス
モンテカルロ法(1/2)
ベイズ統計の難点 （~1980年）
θ 高次元のパラメータ
π (θ | x) の関数形が簡単な形でない場合、以下の積分は困難
Pr(θ ∈ A | x) = ∫ π (θ | x)dθ , E[θ | x] = ∫ θπ (θ | x)dθ
A
→ 今回紹介するモンテカルロ法、特にマルコフ連鎖モンテカルロ法の導入で、
高次元の積分計算も個人PCなどでできるようになった
モンテカルロ法(2/2)
モンテカルロサンプリング
(1)
( 2)
(M )
プログラムを実行して, θ , θ ,  , θ
のようなM個の値が出てきたとする.
うまいプログラムをつくると, あたかも事後分布（確率分布）から値が発生したかのよう
に見える. このように値を抽出する方法をモンテカルロサンプリングという
θ (1) , θ ( 2 ) ,  , θ ( M ) のヒストグラム
π (θ | x) の関数形
posterior distribution
distribution
0.2
0
0.0
0.1
10
M →∞
5
Frequency
0.3
15
0.4
sample histogram
-3
-2
-1
0
1
-2
2
-1
0
1
theta
x
ベイズ統計におけるモンテカルロ法
(1)
( 2)
,,θ ( M )
モンテカルロサンプル θ , θ
を用いて
事後分布（確率分布）に関する様々な量を近似的に計算する方法.
2
モンテカルロサンプルを用いた計算例
θ (1) , θ ( 2 ) ,  , θ ( M )
(a)
1
M
M
∑θ
( j)
j =1
π (θ | x) のモンテカルロサンプル
≈ ∫ θπ (θ | x)dθ
1
M
M
( j)
θ
(
) ≈ ∫ h(θ )π (θ | x)dθ
h
∑
j =1
(b) max{π (θ (1) | x ), π (θ ( 2 ) | x ),  , π (θ ( M ) | x )}
≈ π (θ | x)
のモード（最頻値）
(c)
0.8 から1の間に入る θ ( j ) の数
M
1
≈ ∫ π (θ | x)dθ
0.8
モンテカルロ法は確率分布に関する量の近似計算法
→ 事後分布でなく一般の確率分布として以下では説明
独立同一分布でのモンテカルロ法
基本的な確率分布の場合
モンテカルロサンプルは組み込み関数などから独立に発生
θ
(1)
,θ
1
M
i .i . d .
,,θ
(M )
M
) → ∫ h(θ ) f (θ )dθ
( 2)
∑ h(θ
( j)
j =1
~ f (θ )
M →∞
↑ 数値積分法と対比して, モンテカルロ積分法などという
注意点
・１，２次元の場合は棄却法で任意の確率分布から発生できる(数値積分でもよい）
・サンプリング数
M →∞
で収束するが効率の面から工夫が必要
・変数変換やImportance Sampling
Importance Sampling
工夫が必要な例
X ~ N (0,1)
Pr( X ≥ 20) = E[1{ X ≥ 20} ]
Q1. モンテカルロ積分で計算する方法を説明せよ.
また, その問題点を述べよ
次の恒等式を用いた変数変換が有効
Pr( X ≥ t ) = ∫
t
+∞
1/ t
2
1
e −v / 2 dv = ∫
0
2π
1 1 −1/ 2u 2
e
du
2
2π u
ただし, うまい変数変換を見つける必要がある
Importance Sampling (重点サンプリング)
h(θ ) f (θ )
1
=
≈
h
(
θ
)
f
(
θ
)
d
θ
g
(
θ
)
d
θ
∫
∫ g (θ )
M
h(θ ( j ) ) f (θ ( j ) )
∑
g (θ ( j ) )
j =1
M
Q2. Q1 の例で正規分布の代りにコーシー分布をつかってモンテカルロ
積分をすることができる. h, f, g から被積分関数を具体的に書き下せ.
２. マルコフ核
マルコフ核
X: 有限集合
以下を満たす条件付き分布の族をマルコフ核 （Markov kernel）という.
K ( x | y ) ≥ 0, ∀x ∈ X , ∀y ∈ X
∑ K ( x | y) = 1 ∀y ∈ X
x∈X
マルコフ核はX上の分布をX上の分布へ写す
~p ( x) = K ( x | y ) p( y )
∑
y∈X
~p = Kp
~p ( x) = p( x) = const.
∑
∑
x∈X
← 行列表記 [ K xy ] =
x∈X
K ( x | y)
例題：黄金のブロガー
無料ブログサイトJC2には１０万人程度のブロガー（ほぼ初心者）が登録
している. JC2ではブロガーをゴールド(G)、シルバー(S)、ブロンズ(B)に
クラス分けすることにした. すなわち, 下記の称号を付与し, 称号に応じ
て様々な特典を与えることにした.
ゴールドブロガー, シルバーブロガー, ブロンズブロガー
ただし, 毎週、前の週の閲覧実績や更新頻度などから以下のルールで
クラス間の遷移が発生する.
ゴールド：
７０％が残留；３０％がシルバーに降格
シルバー：
８０％が残留；１０％がブロンズに降格、１０％がゴールドに昇格
ブロンズ：
８０％が残留：２０％がシルバーに昇格
*実在のWebサイトとは一切関係ありません。
例題：黄金のブロガー
注： 登録ブロガーは一定数とし, 第 n 週目の各比率を
p ( n ) (G ), p ( n ) ( S ), p ( n ) ( B)
ゴールド
目指すズラ！
とあらわす.(称号はG,S,B).
第一週は, 全員, S からスタート.
K
K
p ( 2) ( S )
p (1) ( S ) = 1
p ( 2) (G )
G
S B
G
p ( 2) ( B)
S B
X = {G, S , B} 上の分布が毎週、変化していく
G
S B
計算してみよう
１．
( 2)
( 3)
p ( S ), p ( S )
を求めなさい.
p ( S ) = 0.69
p ( S ) = 0.8
( 3)
( 2)
p ( n ) ( x), ( x = G, S , B)
p ( n+1) ( x), ( x = G, S , B)
2．
p
p
から,
を計算する漸化式を書きなさい
( n +1)
(G ) = 0.7 p (G ) + 0.1 p ( S )
( n +1)
( S ) = 0.3 p (G ) + 0.8 p ( S ) + 0.2 p ( B)
p
(n)
( n +1)
(n)
(n)
(n)
( B) = 0.1 p ( S ) + 0.8 p ( B)
(n)
(n)
(n)
挑戦してみよう
３． (a) ２を行列Kとベクトルで書き, Kがマルコフ核になっていることを確認せよ.
p (n) = q
p ( n +1) = q
(b) 今,
の時は,
となり全体の比率は不変である
という. このようなq を求めなさい. (線形代数の問題！)
 q(G )   2 / 11


 
q =  q( S )  =  6 / 11
 q( B)   3 / 11

 

３. マルコフ連鎖と定常分布
マルコフ連鎖
{K ( x | y )}x , y∈X
初期値を
x0 ∈ X
x1 , x2 ,, xt ,
X上のマルコフ核
としてマルコフ核から次のように乱数（確率変数）の列が得られる
xt +1 | xt ~ K ( xt +1 | xt )
このように生成する確率変数の列がマルコフ連鎖 (Markov Chain)
（その意味で {K ( x | y )}x, y∈X をマルコフ連鎖とよぶことも）
マルコフ連鎖が与えられた時, 以下を満たす確率分布を定常分布という
∑ K ( x | y) p( y) = p( x), ∀x ∈ X
y∈X
例題：黄金のブロガー
JC2サイトでは先の行列Kを用いてランダムに
ブロガーの昇格・降格を決めていた。
私の苦労は
何だったの？
事情はさておき・・・
Xn
あるブロガーの n 週目の称号（クラス）
Pr( X n = y ); y = G, S , B
から, n+1 週目に各クラスにいる確率
の計算式は以下のようになる.
Pr( X n +1 = x) =
=
∑ Pr( X
y =G , S , B
n +1
= x | X n = y ) Pr( X n = y )
∑ K ( x | y) Pr( X
y =G , S , B
n
= y)
*実在のWebサイトとは一切関係ありません。
計算してみよう
１． 第一週目にシルバーからスタートしたブロガーが
第二、三週目にブロンズクラスにいる確率を求めよ.
Pr( X 2 = B) = 0.10
Pr( X 3 = B) = 0.69
2．
x1 , x2 ,, xt , という毎週のクラスはマルコフ連鎖となっている. この
マルコフ連鎖に関する定常分布を求めよ.（実はすでに求めている.）
 q(G )   2 / 11

 

q =  q( S )  =  6 / 11
 q( B)   3 / 11

 

基本的な事実
以下では定常分布が唯一存在するケースのみを考える
基本的な事実 １(厳密なStatementは文献参照)
{K ( x | y )}x , y∈X マルコフ連鎖
p∞ (x)
定常分布
p (t ) ( x) = Pr( X t = x) とおくと
p → p∞
( t)
したがって, 初期値
t →∞
x0 に依存しなくなる
考えてみよう
1*．ここでは定常分布が唯一存在することを仮定するが, そのための条件
に関して以下の問に答えよ
Q1. ２つ以上定常分布が存在する例をつくれ
Q2. |X| = ∞ の場合は定常分布が存在しないケースがありえる.
具体的に構成せよ
2． マルコフ核K から, 定常分布を数値的にシミュレーションする方法とし
てどのようなものが考えられるだろうか。
基本的な事実その２
基本的な事実 ２(厳密なStatementは文献参照)
{K ( x | y )}x , y∈X マルコフ連鎖
実現値
p∞ (x)
p∞ (x)
定常分布
x1 , x2 ,, xt , の経験分布は t → ∞ で定常分布
に収束. 特に初期値
x0 に依存しなくなる
p∞ (x) の関数形
x1 , x2 ,, xt , のヒストグラム
0.4
0.3
distribution
15
10
0
0.0
0.1
5
Frequency
posterior distribution
0.2
t →∞
sample histogram
-2
-1
0
x
1
2
-3
-2
-1
0
1
2
theta
MCMC法の核心
確率分布 p (x) を与えたとき、 その確率分布を定常分布にもつような
マルコフ連鎖 {K ( x | y )}x , y∈X を具体的に構成
（定常分布が唯一存在するから, もともとの確率分布に収束してくれる）
４. マルコフ連鎖モンテカルロ法
~ Metropolis-Hastings Method
マルコフ連鎖モンテカルロ法（ＭＣＭＣ法）
発生させたい分布 (Target distribution)
π ( x), ∑ π ( x) = 1, π ( x) ≥ 0
x
理論
∑ K ( x | y)π ( y) = π ( x), ∀x ∈ X
y∈X
を満たすマルコフ核 {K ( x | y )}x , y∈X を与え
x1 , x2 ,, xt ,を以下、順次発生
xt +1 | xt ~ K ( xt +1 | xt )
実際には（ほとんどのテキスト）
モンテカルロサンプル x1 , x2 ,  , xt ,  を発生させるアルゴリズム
（マルコフ核は実は煩雑）
資料ではアルゴリズムとマルコフ核の両方を紹介（時間の都合でとばすかも）
紹介するアルゴリズム
代表的な２つ
- MH法 (Metropolis-Hastings Method)
- Gibbs Sampler (Gibbs Sampling)
MH法
- 詳細つり合い条件 (detailed balance condition) を満たすマルコフ核
独立MH法, 一般MH法
Gibbs Sampler
- 定常条件を直接満たすマルコフ核（詳細つり合い条件は満たさないかも）
二段階GS, 多段階GS
詳細つり合い条件
詳細つりあい条件 (Detailed Balance Condition)
K ( y | x)π ( x) = K ( x | y )π ( y )
詳細つりあい条件を満たす分布
π
∀x∀y
は, 実は定常分布
Ｑ．詳細つりあい条件を満たす分布が定常分布になっていることを示せ
MCMC法の基本原理
π
について詳細つりあい条件を満たすマルコフ核を K ( x | y ) とする.
x0 から生成するマルコフ連鎖 x1 , x2 ,, xt ,
の経験分布 (ヒストグラム)は, t → ∞ で, π に収束
任意の初期値
Ｑ．「基本的な事実その２」を用いて, 上を示しなさい.
参考：詳細つり合い条件
状態数が２（Kが２×２行列） の場合、DBCと定常条件は同値
状態数が３以上の場合, DBCを満たさない場合がある。
例：
 1/ 3 1/ 3 1/ 4


K ( x | y) =  2 / 3 0 1 / 4 
 0 2 / 3 1/ 2


１．定常分布は唯一であり以下で与えられる（DBCをみたさない）
 3 / 10 


π (x) =  3 / 10 
 2/5 


２．DBCを満たす分布は存在しない
２点間のMH法 (1/2)
２点間のMH法
π ( y ) > π ( x)
π (y )
π (x)
今,
zt (=x, y) が所与の下
次の行き先 zt +1 を確率的に定める.
(a) 行き先候補
y
x
{x, y} の中からランダムに選ぶ
= x) = P( z prop = y ) = 1 / 2
z prop を
P( z prop
(b) バランスを崩さないように本当に移動するか, 踏みとどまるか決める
π ( z prop )
≥1
π ( zt )
π ( z prop )
<1
π ( zt )
zt +1 = z prop
←覚え方： 転職先の給与額が高ければ迷わず移動
π ( z prop )
確率 π ( z ) で zt +1 = z prop
t
π ( z prop )
1
−
z = zt
確率
π ( zt ) で t +1
Rでの実装例(1/2)
２点間のMH法の実装例
π (1) = 1 / 3, π (2) = 2 / 3
z=1,2 をMH法に基づいてランダムに動き回る
シミュレーション
プログラム例
f_prior <- function(x){
if(x==1){ return(1/3) }; if(x==2){ return(2/3) };
return(0);
}
Nsim <- 10^3;z <- array(0, Nsim);
# z= 1, 2 をうろうろする.
z[1] <- 1; # initial value
for(t in 1:(Nsim-1)){
z_prop <- sample(1:2, 1);
if ( f_prior(z_prop)/f_prior(z[t]) >= runif(1) ){
z[t+1] <- z_prop ; }else{z[t+1] <- z[t];}
}
π (y )
π (x)
y
x
Rでの実装例(2/2)
２点間MH法の実行結果をプロット
π (1) = 1 / 3
Two-point MH algorithm
2
1.0 1.4 1.8
z[1:50]
π (2) = 2 / 3
1
0
10
20
30
40
50
Two-point MH algorithm
1.0 1.4 1.8
z[950:1000]
time
950
プログラム例
960
970
980
990
1000
time
par(mfrow=c(2,1));
plot(1:50, z[1:50], type="b",xlab="time", main="Two-point MH algorithm");
plot(950:1000, z[950:1000], type="b",xlab="time", main="Two-point MH algorithm");
２点間のMH法 (2/2)
マルコフ核
先に紹介したアルゴリズムは
v, w ∈ x, y として以下のマルコフ核に相当
{
}

 π (v ) 
1
min 1,


2
 π ( w) 

K (v | w) =  1 1

π (u ) 
1 −

 +
∑
π ( w) 
 2 2 u: π ( u ) <1 
π ( w)

Q1.
∑ K (v | w) = 1
v≠w
v=w
を示せ
v
Q2.
K (v | w)π ( w) = K ( w | v )π (v ) （詳細つりあい条件）を示せ
Q2’. π (1) = 1 / 3 π (2) = 2 / 3 のとき,マルコフ核を行列であらわせ.
独立MH法 (1/2)
簡単のため有限集合
今,
Θ
上の分布を考える
zt ∈ Θ が所与の下 次の行き先 zt +1 を確率的に定める.
(a) 行き先候補
z prop ~ q ( z )
∑ q( z ) = 1, q( z ) ≥ 0
z∈Θ
q(z): 提案分布 (自由に設定)
(b) 本当に移動するか踏みとどまるか決める
ρ=
q( zt )π ( z prop )
q( z prop )π ( zt )
← 一般のMH法は詳細つりあい条件を満たすようにここが複雑になる
ρ ≥1
zt +1 = z prop
ρ <1
確率
確率
ρ
1− ρ
で
zt +1 = z prop
で
zt +1 = zt
独立MH法 (2/2)
マルコフ核
v, w ∈ Θ
（簡単のため離散分布で表現）

 q ( w)π (v) 
q (v) min 1,


 q (v)π ( w) 

K (v | w) = 

q ( w)π (u ) 

q (u )1 −
q ( v ) +
∑
q (u )π ( w) 
q ( w )π ( u )


<1
u:
q ( u )π ( w )

Q3.
∑ K (v | w) = 1
v≠w
v=w
を示せ
v∈Θ
Q4.
K (v | w)π ( w) = K ( w | v )π (v ) （詳細つりあい条件）を示せ
提案分布は現在の場所
zt ∈ Θ に依存してもよい（次に紹介）
一般MH法
今,
zt ∈ Θ が所与の下 次の行き先 zt +1 を条件付き分布で定める.
(a) 行き先候補
z prop ~ q ( z | zt )
(b) 本当に移動するか踏みとどまるか決める
ρ=
q( zt | z prop )π ( z prop )
q( z prop | zt )π ( zt )
ρ ≥1
zt +1 = z prop
ρ <1
確率
ρ
確率
1− ρ
zt +1 = z prop
で
で
zt +1 = zt
Q5. 独立MH法と同様にマルコフ核を書き下せ, また詳細つりあい条件を
満たすことを示せ
独立MH法の実装例
ベータ分布 Be(2.7, 6.3)
MCMC利用（MH法）とRの組み込み関数rbeta() の比較
200
100
50
0
0.0
0.2
0.4
0.6
0.8
X
100
200
300
histogram of rbeta
0
50
set.seed(1234);
a <- 2.7; b <- 6.3; c <- 2.669;
Nsim <- 5000;
X <- rep(runif(1), Nsim);
for (i in 2:Nsim){
Y <- runif(1);
rho <- dbeta(Y, a, b)/dbeta(X[i-1], a, b) ;
X[i] <- X[i-1] + (Y - X[i-1] )* (runif(1) < rho);
}
S <- Nsim/40;
MAX_Y <- S*max(dbeta((1:10000)/10000,a,b));
hist(X, col=2, breaks=40, ylim = c(0, MAX_Y),
main="histogram of MCMC");
curve(S*dbeta(x, a,b), add=T);
Frequency
プログラム例 (MCMC)
Frequency
300
histogram of MCMC
# Rの組み込み関数を用いる場合
hikaku <- rbeta(Nsim, a,b);
0.0
0.2
0.4
hikaku
0.6
0.8
理解度チェック：MCMC法
MCMC法を提案する場合、発生させたい確率分布（事後分布）につい
て、詳細つり合い条件を満たすマルコフ核を与える必要がある
モンテカルロサンプルの定常分布への収束は一般の場合で証明
できるが、早く収束するような工夫は個別に必要である
MH法のアルゴリズムをマルコフ核で書き下すと、発生させたい確
率分布は、常に詳細つりあい条件を満たしている。
独立ＭＨ法では、次の行き先を決める分布（提案分布）は現在の
状態と独立である
MCMC法でサンプル数を増やしていくと, ほとんど一定の値しか
出てこなくなり, 収束していく。
５. マルコフ連鎖モンテカルロ法
~ Gibbs Sampler
*Gibbs Sampling ともいう
Gibbs Sampler
Geman and Geman (1984)
応用でGibbs Samplerを使った
Gelfand and Smith (1990)
Gibbs Sampling を紹介, 統計学者の間で認知され広まる
基本的な使用の方針
周辺分布や同時分布が簡単に書けない（サンプリングが難しい）が、
条件付き分布はサンプリングしやすい（正規分布など）場合に使う
Two Stage Gibbs Sampling
Two Stage Gibbs Sampling
サンプリングしたい確率分布 （x, yの同時分布！）
f ( x, y )
二つの条件付き分布 （こちらは比較的簡単な形）
fY | X ( y | x), f X |Y ( x | y )
Algorithm
1. 初期値を選ぶ
X 0 = x0
2. 条件付き分布から逐次的に
xt , yt
Yt | X t −1 = xt −1 ~ fY | X (⋅ | xt −1 )
X t | Yt = yt ~ f X |Y (⋅ | yt )
を生成
Gibbs Samplingの実装例 (1/2)
確率分布 (GSで発生させたい分布)
exp(−θ / 2)
f (θ ) = C
(1 + θ 2 )ν
(*この例は棄却法などでも可能)
2
やりかた
1. 補助変数 (
η
ν =3
) を加えた積分表示を求めて, 同時分布
f (θ ,η )
を考える
 1 + θ 2  ν −1
f (θ ) = D ∫ exp(−θ / 2) exp −
η η dη
0
 2 
∞
2
2. 次の二つの条件付き分布を求める
 1 +η 2 
f (θ | η ) =
θ 
exp −
−1

 2
2π (1 + η )
1
3. Two Stage GS
ν
1  1 + θ 2  ν −1  1 + θ 2 

 η exp −
f (η | θ ) =
η 
Γ(ν )  2 
 2 
Gibbs Samplingの実装例 (2/2)
1000
Gibbs sim
プログラム例
Nsim <- 5000;
theta <- eta <- array(0, dim=c(Nsim, 1));
nu <- 3;
theta[1] <- rnorm(1);
eta[1] <- rgamma(1, shape=nu, rate= (1+theta[1]^2)/2);
for (i in 2:Nsim){
theta[i] <- rnorm(1, mean=0, sd=1/sqrt( 1 + eta[i-1]));
eta[i] <- rgamma(1, shape=nu, rate= (1+theta[i]^2)/2);
}
600
0
200
400
ν =3
Frequency
exp(−θ 2 / 2)
f (θ ) = C
(1 + θ 2 )ν
800
確率分布 (GSで発生させたい分布)
-3
-2
-1
0
1
2
3
theta
点線は理論曲線をヒストグラムにあわせ
て描いたもの
2次元同時事後分布を作る場合
同時分布を求めたい場合
確率分布
X i ~ N (θ , v)
θ ~ N (θ 0 , w)
v ~ IG (a, b)
i = 1,2,, n
事後分布 (θ , v) ~ π (θ , v | x) の直接のサンプリングはめんどう （MH法）
しかし, 次の二つの条件付き分布を使うのは（組み込み関数が使えるので）容易
θ | v, x ~ N (ξ ,V ), v | θ , x ~ IG( A, B)
(nw) xn + vθ0
vw
,ξ =
V=
v + nw
v + wn
やりかたは同じ
S + n( x − θ ) 2
A = n / 2 + a, B = b +
2
n
S = ∑ ( xn − x j ) 2
j =1
計算部分
補足（計算部分）
確率分布
X i ~ N (θ , v)
θ ~ N (θ 0 , w)
v ~ IG (a, b)
(θ , v) ~ π (θ , v | x)
i = 1,2,, n
 S + n( x − θ ) 2   (θ 0 − θ ) 2  1 −b / v
 a +1 e
 exp −
π (θ , v | x) ∝ n / 2 exp −
2v
2w  v
v

 
1

 1
  1
(θ − xn ) 2  exp − (θ − θ 0 ) 2 
π (θ | v, x) ∝ exp −

 2v / n
  2w

 1
∝ exp − (θ − ξ ) 2 

 2V
π (v | θ , x )
(n / v) xn + w−1θ 0
1 1 1
=
+ ,ξ =
n / v + w−1
V v/n w
  S + n( x − θ ) 2 b  
∝ n / 2+ a +1 exp− 
+ 
v
2v
v 
 
1
v | θ , x ~ IG ( A, B)
S + n( x − θ ) 2
A = n / 2 + a, B = b +
2
Multi Stage Gibbs Sampling
Multi Stage Gibbs Sampling
サンプリングしたい確率分布
f ( x1 , x2 ,, x p )
条件付き分布
f1 ( y | x2 ,, x p ), f 2 ( y | x1 , x3 ,, x p ) ,, f p ( y | x1 , x2 ,, x p −1 )
Algorithm
1. 初期値を選ぶ
(0)
x1 ,, x p
2. 条件付き分布から逐次的に
(0)
(t )
(t )
x1 , x2 ,, x p
x1(t ) ~ f1 (⋅ | x2(t −1) ,, x (pt −1) )
x2(t ) ~ f 2 (⋅ | x1(t ) , x3(t −1) ,, x (pt −1) )

x (pt ) ~ f p (⋅ | x1(t ) , x2(t ) ,, x (pt−)1 )
(t )
を生成
Gibbs Samplingの収束性
証明の概略
マルコフ連鎖
x ( 0) ,..., x (t ) , x (t +1) ,
のマルコフ核Kについて
fが定常分布であることを示せば十分. (詳細は文献参照)
p=3 で示す（一般の場合も同様）
Gibbs sampler で生成される
x (t )
は以下のマルコフ核をもつマルコフ連鎖
K ( z | y ) = f 3 ( z3 | z1 , z2 ) f 2 ( z2 | z1 , y3 ) f1 ( z1 | y2 , y3 )
また,
y ~ f ( y ) の時,
z の分布は
z ~ g ( z ) = ∫ K ( z | y ) f ( y )dy
ところがこれは計算すると
f (z )
に一致する.
したがって, f はマルコフ連鎖の定常分布
※GSのマルコフ核は一般に詳細つりあい条件は満たしてない
計算部分
p=3 で g ( z ) = ∫ K ( z | y ) f ( y )dy = f ( z ) を示す
やや一般的に以下を計算 (A は任意のz可測集合)
P( z ∈ A) = ∫ 1A ( z ) K ( z | y ) f ( y )dydz
= ∫ 1A ( z ) f 3 ( z3 | z1 , z2 ) f 2 ( z2 | z1 , y3 ) f1 ( z1 | y2 , y3 ) f ( y1 , y2 , y3 )dzdy
ここで,
∫∫ f ( z | y , y ) f ( y , y , y )dy dy = ∫ f ( z | y , y ) f ( y , y )dy
= ∫ f ( z , y , y )dy
1
1
2
3
1
2
3
1
2
1
1
1
2
2
3
3
2
3
2
2
= f ( z1 , y3 )
といった計算を繰り返すと
= ∫ 1A ( z ) f 3 ( z3 | z1 , z2 ) f 2 ( z2 | z1 , y3 ) f ( z1 , y3 )dz1dz2 dz3dy3
= ∫ 1A ( z ) f 3 ( z3 | z1 , z2 ) f ( z1 , z2 , y3 )dz1dz2 dz3dy3
= ∫ 1A ( z ) f ( z1 , z2 , z3 )dz1dz2 dz3
したがって
z ~ f ( z)
理解度チェック：MCMC法２
高次元のパラメータに関する積分であっても、組み込み関数で簡単
に独立乱数を発生できるなら、MCMC法は必要ない
f(x,y)でのTwo Stage Gibbs Sampling では、十分、xのサンプルを独立
に発生させた後で、条件付分布からy のサンプルを発生させる.
Two Stage Gibbs Sampling では、同時分布f(x,y)からのサンプリング
を利用して、条件付分布のモンテカルロサンプルを求める
Two Stage, Multi Stage Gibbs Sampling は, 条件付分布を利用して、
逐次的にサンプリングするという点で同じである。
条件付き確率分布が簡単な形で書けない場合は、Gibbs Sampling
を利用するしかない
６.WinBUGS
BUGS
BUGS=Bayesian inference Using Gibbs Sampling
WinBUGS = Windows上で動くバイナリー (closed source)
・統計モデル
・事前分布
・データ
・初期値
・固定パラメータ
MCMCサンプル
0.05155622 0.07840824
0.06736837 0.08869456
0.03744929 0.02129103
0.08448234 0.05873367
0.02704619 ....
基本的な使用の方針
メジャーな確率分布を複数組み合わせて解析する場合に使う
（たとえば正規分布のみでも複数組み合わせると条件付き分布の計算が面倒になってくる）
*Open SourceのBUGSもある(OpenBUGS)
Rとの連携
WinBUGSの欠点
・付属のGUIが使いにくい
・MCMC出力を使った処理が自由にできない
Rのインターフェース利用
Rプログラムの中で呼び出して使うことができる
(管理者として実行しないとエラーになることも)
・統計モデル
・事前分布
・データ
・初期値
・固定パラメータ
MCMCサンプル
0.05155622 0.07840824
0.06736837 0.08869456
0.03744929 0.02129103
0.08448234 0.05873367
0.02704619 ....
ExcelやMATLAB, Stata, etc. からも使える
Example: Coal Mining (1/2)
データ
イギリスの炭鉱での災害数
１８５１年～１９６２年 （Carlin et al., 1992）
4
3
2
1
0
> N <- 112;
> D <- c(4,5,4,1,0, 4,3,4,0,6,
3,3,4,0,2, 6,3,3,5,4, 5,3,1,4,4, 1,5,5,3,4, 2,5,2,2,3, 4,2,1,3,2,
1,1,1,1,1, 3,0,0,1,0, 1,1,0,0,3, 1,0,3,2,2,
0,1,1,1,0, 1,0,1,0,0, 0,2,1,0,0, 0,1,1,0,2,2,3,1,1,2, 1,1,1,1,2,
4,2,0,0,0, 1,4,0,0,0,1,0,0,0,0, 0,1,0,0,1, 0,0);
> plot(1851:1962, D, main="Example of Poisson Change
Point Model", xlab="Year", ylab="Num. of Accident")
Num. of Accident
5
6
Example of Poisson Change Point
1860
1880
1900
1920
1940
Year
ここでは変化点問題として、２つの独立なポアソンモデルからのデータ発生と
みなしてモデリングする （今回はWinBUGSの使用方法の説明として採用）
1960
Example: Coal Mining (2/2)
統計モデル
Dt ~ Po( µt ) log µt = b1 + H (t − τ )b2
b1 ~ N (0,106 ), b2 ~ N (0,106 ),
τ ~ U (1,112)
← BUGS コードでは, 分散の逆数を指定する
−6
ので 10
が出てくる
← 変化点は本来、離散値だが、連続値として一様分布
をふっている
model のBUGS言語での記述 (coalmining.txt)
model
{
for(year in 1:N){
D[year] ~ dpois(mu[year])
log(mu[year]) <- b[1] + step(year - changeyear) * b[2]
}
for(j in 1:2){
b[j] ~ dnorm(0.0, 1.0E-6)
}
changeyear ~ dunif(1, N)
}
R+WinBUGSの実装例
実装例:
・ bugs 関数 (codaPkg=TRUE)でWinBUGSを実行、さらに read.bugs
・ plot でトレースプロットと事後分布（を推定したもの）が表示される.
0
1
2
3
0.9 1.1 1.3
700
800
900
1000
0.8
1.0
1.2
1.4
Iterations
N = 500 Bandw idth = 0.02272
Trace of b[2]
Density of b[2]
0.0
-1.8
1.0
-1.2
2.0
-0.6
600
600
700
800
900
1000
-1.8
-1.6
-1.4
-1.2
-1.0
-0.8
-0.6
Iterations
N = 500 Bandw idth = 0.03829
Trace of changeyear
Density of changeyear
35
40
45
0.00 0.10 0.20
500
700
800
900
1000
35
40
45
Iterations
N = 500 Bandw idth = 0.5065
Trace of deviance
Density of deviance
0.10
0.20
600
355
500
345
0.00
coalmining.sim <- bugs(data, inits,
parameters, "coalmining.txt", n.chains=3,
n.iter=1000, bugs.directory="c:/Program
Files/WinBUGS14/", codaPkg=TRUE);
coalmining.coda <read.bugs(coalmining.sim);
plot(coalmining.coda);
summary(coalmining.coda);
500
335
data <- list("N", "D");
parameters <- c("changeyear", "b");
inits <- function(){ list(b=c(0,0),
changeyear=50) };
Density of b[1]
4
Trace of b[1]
500
600
700
800
Iterations
900
1000
335
340
345
350
N = 500 Bandw idth = 0.5497
355
MCMCサンプルの取り出し (1/2)
MCMCサンプルの取り出し方法 １
（収束が確認できた後で） bugs 関数 (codaPkg=FALSE)でWinBUGSを実行
→ 出力オブジェクトから取り出す
> coalmining.sim.dat <- bugs(data, inits, parameters, "coalmining.txt", n.chains=3, n.iter=1000,
bugs.directory="c:/Program Files/WinBUGS14/", codaPkg=FALSE);
> coalmining.sim.dat$sims.list$changeyear[1:10]
[1] 40.60 40.69 40.81 41.77 37.52 37.67 40.40 38.97 38.31 36.49
> coalmining.sim.dat$sims.list$b[1:10,]
[,1] [,2]
[1,] 1.362 -1.4190
[2,] 1.193 -1.1050
[3,] 1.095 -1.1670
[4,] 1.049 -0.9935
[5,] 1.230 -1.3540
[6,] 1.120 -1.2310
[7,] 1.127 -1.4230
[8,] 1.184 -1.3800
[9,] 1.275 -1.3890
[10,] 1.405 -1.3830
MCMCサンプルの取り出し (2/2)
MCサンプルの取り出し方法 ２
codaPkg=TRUE で実行, read.bugs()でリストを生成.
→ リストの中身を直接、参照する.
> coalmining.coda[1:3,]
[[1]]
b[1] b[2] changeyear deviance
[1,] 1.175 -1.134 37.20 336.6
[2,] 1.243 -1.538 41.94 339.3
[3,] 1.242 -1.481 38.11 340.2
[[2]]
b[1] b[2] changeyear deviance
[1,] 0.9676 -1.153 43.10 341.5
[2,] 1.0320 -1.225 40.31 335.6
[3,] 1.1270 -1.016 36.00 338.6
[[3]]
b[1] b[2] changeyear deviance
[1,] 1.256 -1.482 37.49 337.8
[2,] 1.300 -1.561 37.73 339.7
[3,] 1.376 -1.668 39.84 343.4
MCMCの収束診断 (1/2)
トレースプロット
時系列データとして各パラメータ値をプロットしたもの
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
x
a= 2.7 :b= 6.3 :c= 2.669
4500
4550
4600
4650
4700
4750
4800
step
トレースプロットの例
トレースプロットでのチェック
１．トレンドがはっきり見えると、確率分布から独立標本とほど遠いのでダメ
２．自己相関関数(時系列データとして)を計算してすそが０に近づかない場合はダメ
３．初期値を変えて数回走らせる
→ 時系列でプロットしたときに振る舞いが明らかに違うとダメ
MCMCの収束診断 (2/2)
注意点
・明らかに収束してないケースをはじくことは容易
・パラメータの次元が大きい場合すべてのパラメータの視覚的な収束診断は不可能
・収束してないけど、さまざまな判断基準をパスする可能性もある
（この条件を満たせば収束OK といえる、万能な判断基準はない！！）
理論とのギャップ
(適切なMCMCの設定の下では)MCMCサンプルの数 → ∞ で
サンプルの分布は、必ず求めたい事後分布に収束
→ しかし、MCMCを使う必要がある場合、収束の速さなどもうまく評価できない
まとめ
MCMCで分析したいとき
・WinBUGSでの解析例 (統計モデル, BUGS のコード)が次第に増えてきて
いる → Exampleで、解析したいケースと類似のケースを探してみる
・類似例が見つからない、もしくはうまくいかない場合
→ 自分で統計モデルをたてるしかない
→ WinBUGSを経由せず（Rなどで）いちからＧＳ, ＭＨ法などでプログラムを作る
のがベストかも【今回、紹介した以外にも方法はある】
ベイズモデリングの指針
・（共役事前分布を用いるなどして）陽な形で事後分布がかけるモデルで
やってみる
・階層ベイズにすることで、かなり柔軟な統計モデルがつくれる