1
Random Generation of
Bm×J
Contingency Tables
ver. 5.0
東京大学 松井知己
東海大学 松井泰子
東京理科大学 小野陽子
2
Random Generation of
｛１,2｝m×｛1,2,...,n｝
Contingency Tables
東京大学 松井知己
東海大学 松井泰子
東京理科大学 小野陽子
3
Path Coupling 法を用いた
多元分割表生成のためのマルコフ連鎖設計
東京大学 松井知己
東海大学 松井泰子
東京理科大学 小野陽子
4
本研究の目的
題目：Path Coupling 法を用いた
多元分割表生成のためのマルコフ連鎖設計
扱う分割表の形状： 2×2×･･･×2×K
2×2×K：[太田,青木，広津(昨日)]
マルコフ基底：自明なものが存在
3×3×K：[青木,竹村(昨日)]
極限分布： (1)一様分布 (2)超幾何分布
目的：1つサンプリングする迄に必要な
推移回数の算定 ＝mixing time の算定
5
MCMC法
多元分割表の行と列の独立性の検定
検定理論
正確検定
医療統計
p値の計算
帰無仮説の棄却を決める閾値
MCMC法（Markov Chain Monte Carlo Method)
Markov chain で sampling を行い
Monte Carlo 法を行う
mixing time の速いMarkov chain の必要性
6
contingency table
{1,2}2×{1,2,3,4,5,6}の例
列 (column)
7
contingency table (example)
２
６
８
５
６
１１
５ １
３ ６
８ ７
２ ７
４ ４
６ １１
４
３
７
１
３
４
５ ９
２ ０
７ ９
０ ３
８ ８
８ １１
２６
２０
１８
３３
７ ７ ８ ５ ５ １２
１２ ７ １０ ６ １０ ８
N：table sum
(N=97=26+20+18+33)
n: # of columns (n=6)
8
contingency tables
*
*
８
*
*
１１
* *
* *
８ ７
* *
* *
６ １１
*
*
７
*
*
４
* *
* *
７ ９
* *
* *
８ １１
２６
２０
１８
３３
７ ７ ８ ５ ５ １２
１２ ７ １０ ６ １０ ８
N：table sum
(N=97=26+20+18+33)
n: # of columns (n=6)
9
分布関数
周辺和を固定したとき、
T =（ nijk ）という分割表が出現する確率が、
（1）一様分布すると仮定 ←Diaconis-Efron検定
Pr[T]＝1/ |Ω|
Ω：与えられた表と周辺和が等しい分割表の集合.
（２）超幾何分布すると仮定 ←対数線形モデル
Pr[T]＝Pr[nijk]＝
Δ
Πi Πj Πk nijk !
Δ：総和が1となる係数
10
main result
目的
｛1,2｝m×｛1,2,...,n｝分割表を ランダム生成する
（定常分布が一様分布と超幾何分布の）
Markov chain の提案
結果
提案する Markov chain の mixing time は,
・分割表の列数,
・ln(分割表のセル中の数値の平均),
・ln (精度の逆数)
の多項式でおさえられる.
方法
path coupling 法 [Bubley and Dyer 1997]
11
2元分割表 (contingency table)
2元分割表
eye Blue
hair
Blond
3
Brown
2
Red
total
4
9
Brown Green total
4
4
3
2
10
8
3
11
5
10
12
30
行と列の独立性を検定したい.
total: 周辺和
12
行と列の独立性
行と列の独立性を検定したい. H0: pij = pi･ p･j
eye Blue
hair
Blond
Brown
Red
3
2
4
9
eye Blue
hair
Blond
Brown
Red
p･ j
9･10
/900
9･8
/900
9･12
/900
9/30
Brown
4
4
3
11
Brown
11･10
/900
１1･８
/900
11･12
/900
11/30
Green
3
2
5
10
10
8
12
30
Green
p i･
10･10
/900
10･8
/900
10･12
/900
10/30
10/30
8/30
12/30
1
pij
13
行と列の独立性の検定
行と列の独立性を検定したい.
eye Blue
hair
Blond
Brown
Red
n+j
3
2
4
9
eye Blue
hair
Blond
Brown
Red
9･10
/30
9･8
/30
9･12
/30
9
Brown
4
4
3
11
Brown
11･10
/30
１1･８
/30
11･12
/30
11
Green
ni+
3
2
5
10
10
8
12
30
nij
Green
10･10
/30
10･8
/30
10･12
/30
10
10
8
12
30
mij
14
正確検定
X2=Σi Σj（nij – mij)2/mij ： χ2値
サンプル数が多ければ, X2は 自由度が
(行数-1)(列数-1)のχ2分布に従う.
サンプル数が少ないとき（正確検定)：
Ω：与えられた表と周辺和が等しい2元分割表の集合.
T∈Ω, X2(T)：Tのχ2値
Σ{Pr[T]： X2<X2(T), T∈Ω} < α
→ 有意水準αで棄却
モンテカルロ法
Ω中の表をランダムに生成して,
Σ{Pr[T]： X2<X2(T), T∈Ω} の値を計算する.
15
main theorem
提案するMarkov chain の反復回数について,以下の
定理を得た.
定理
τ(ε)≦(1/2)n(n-1) ln (ndεｰ1）
ε：分布に対する精度
(正確な定義は後述)
τ(ε): mixing time
(Markov chain の反復回数)
(正確な定義は後述)
N：分割表中の数値の総和
n, m : 分割表のサイズは{1,2}m×{1,2,...,n}
d=N/(n2m ) ：セル内の数値の平均
16
２元分割表の最近の研究（主にアルゴリズム分野)
2元分割表(m×n)：
N: セル内の数値の総和
一様分布：
Diaconis, Saloff-Coste (1995)
poly(N, ln(1/ε))･････(n, m を定数と見なしたとき)
Chung, Graham, Yau (1996)
poly(m, n, N, ln(1/ε))･･･(周辺和が十分大きいとき)
Dyer, Greenhill (2000)
poly(n, lnN, ln(1/ε)) ･･････････(2×n分割表)
Open problem: ∃? poly(m, n, ln N, ln(1/ε))
多項超幾何分布：
Markov chainを用いずにO(N)で perfect sampling 可能
17
3元分割表の最近の研究（主にアルゴリズム分野)
3元分割表(n×m×l)：
Irving, Jerrum (1994)
周辺和を満たす表の存在判定問題は NP-完全
Markov 基底の生成：
Diaconis and Strumfels (1998)
グレブナー基底を用いた基底の生成法
Aoki, Takemura (2002)
3×3×J の基底の提案
多項超幾何分布：
Jerrum, Sinclair, Vigoda (2001)
poly(m,l, ln N, ln(1/ε))
････2部グラフの完全マッチングのサンプリング
（２×I×J ,周辺和は0or1の中の更に特殊ケース)
18
4. 一様分布を定常分布に持つ
Markov chain の提案
an extension of
Dyer and Greenhill’s chain [２０００]
19
Markov chain M0 (extension of Diaconis and Saloff-Coste’s chain)
M0: 列{j’, j’’} （j’≠j’’)をランダムに選ぶ
２
６
５
６
５
３
２
４
１
６
７
４
５
２
０
８
９
０
３
８
０ +１
０ -１
０ -１
０ +１
０ -１
０ +１
０ +１
０ -１
０
０
０
０
０
０
０
０
０ ‐１
０ +１
０ +１
０ -１
０ +１
０ -１
０ -１
０ +１
０
０
０
０
０
０
０
０
２
６
５
６
１
６
７
４
５
２
０
８
９
０
３
８
２
６
５
６
１
６
７
４
５
２
０
８
９
０
３
８
６
２
１
５
４
３
１
３
３
４
２
２
1/2
1/2
推移できないときは,
列選択から再試行
４
４
３
３
５
２
０
４
20
new Markov chain M1 (extension of Dyer and Greenhill’s chain)
M1: 列{j’, j’’} （j’≠j’’)をランダムに選ぶ
３
５
４
７
２
６
５
６
５
３
２
４
１
６
７
４
５
３
２
４
０
７
８
３
４
３
１
３
２
６
５
６
４
３
１
３
５
２
０
８
５
３
２
４
５
２
０
８
-θ
+θ
+θ
-θ
０ +θ ０
０ -θ ０
０ -θ ０
０ +θ ０
０
０
０
０
０
０
０
０
０
８
７
４
４
３
１
３
５
３
２
４
５
２
０
８
３
４
５
６
９
０
３
８
θ=
９
０
３
８
１
６
７
４
９
０
３
８
-1
0
４
３
１
３
５
２
０
８
９
０
３
８
1
2
１
７
６
５
５
３
２
４
２
５
６
５
４
３
１
３
５
２
０
８
９
０
３
８
21
new Markov chain M1
M1: 列{j’, j’’} （j’≠j’’)をランダムに選ぶ
３
５
４
７
２
６
５
６
５
３
２
４
１
６
７
４
５
３
２
４
０
７
８
３
４
３
１
３
２
６
５
６
４
３
１
３
５
２
０
８
５
３
２
４
５
２
０
８
９
０
３
８
０ +θ ０
０ -θ ０
０ -θ ０
０ +θ ０
０
０
０
０
０
０
０
０
1/4 ０
８
1/4
1/4 1/4
７
４
１ ４ ５ ９
１ ５ ２ ４
６ ３ ２ ０
７ ３ ５ ３
７ １ ０ ３
６ ２ ６ １
４ ３ ８ ８
５ ４ ５ ３
５
３
２
４
５
２
０
８
３
４
５
６
９
０
３
８
９
０
３
８
-θ
+θ
+θ
-θ
４
３
１
３
５
２
０
８
９
０
３
８
22
geometric interpretation
分割表の集合
＝24次元空間中の多面体中の整数格子点
…
…
…
…
x 1 x2
x 7 x8
x13
x19
x6
x18
x24
x24
x4
O
x3
x2
x1
23
geometric interpretation of M0
new Markov chain M0
j’
j’’
２ ５ １ ４ ５ ９
６ ３ ６ ３ ２ ０
５ ２ ７ １ ０ ３
６ ４ ４ ３ ８ ８
1/2
x24
1/2
x4
O
x3
x2
x1
24
geometric interpretation of M1
new Markov chain M1
j’
j’’
２ ５ １ ４ ５ ９
６ ３ ６ ３ ２ ０
1/4
５ ２ ７ １ ０ ３
６ ４ ４ ３ ８ ８
1/4
x24
1/4
1/4
x4
O
x3
x2
x1
θ
25
properties of new Markov chain
提案する Markov chain M1 の特徴
Markov chain M1 は, エルゴード的.
1. finite
2. a-periodic
3. irreducible （｛1,2｝m×｛1,2,...,n｝という性質)
定理
M1 の 定常分布は, 与えられた分割表全体の
集合上の一様分布である.
（証明略）
26
3. coupling method
[Aldous 1983]
Markov chain の mixing time を算定する手法
27
coupling sample path
2つのサンプルパスから， 新しいサンプルパス
を作る．
状態空間
x
y
0
時刻
28
coupling sample path
2つのサンプルパスから， 新しいサンプルパス
を作る．
状態空間
状態が最初に一致した時点で
他方のsample path に乗り換える
x
y
0
時刻
Markov chain の無記憶性
→ y から出発したsample path と同じに見える！
29
coupling sample path
仮想的にすべての状態からサンプルパスを発
生させる．
状態空間
x
0
時刻
パスが一つになった
→ 初期状態に依存していない． → 定常分布
30
coupling (definition)
coupling の定義：
与えられた Markov chain M の coupling
(Ω：Mの状態空間の集合)
1. Markov chain
2. 状態空間：Ω×Ω=Ω2
3. coupling の推移確率は次を満たす.
∀ (x , y ) ∈Ω2,
PM(x, x’)=∑｛P((x, y),(x’, y’)) | y’∈Ω｝
PM(y, y’)=∑｛P((x, y),(x’, y’)) | x’∈Ω｝
P((x, y),(x’, y’)) ：coupling での (x, y)から(x’, y’)への
推移確率
PM(x, x’)：Mにおける x から x’への推移確率
31
coupling (marginal distribution)
∀ (x , y ) ∈Ω2,
PM(x, x’)=∑｛P((x, y),(x’, y’)) | y’∈Ω｝
PM(y, y’)=∑｛P((x, y),(x’, y’)) | x’∈Ω｝
y∈Ω
x∈Ω
Ω2
32
coupling (as a simulation)
初期状態(x, y)からのcoupling の推移
=Markov chain M を 異なる初期状態 xとy から推移
したときのシミュレーション
y∈Ω
x∈Ω
(x,y)∈Ω2
33
coupling (初期状態への依存)
任意の初期状態 (x, y) から， 対角要素 (x’,x’) へ
coupling で推移するまでの時間を算定する．
初期状態に依存しない．
y∈Ω
x∈Ω
x’
(x,y)∈Ω2
34
coupling (non-trivial coupling)
non-trivial coupling:
case x＝y: Ω*上にMを貼り付ける
P((x, x),(x’, x’)) = PM(x, x’)
Ω*
P((x, x),(x’, y’)) = 0
(x’≠y’)
case x≠y :
P((x, y),(x’, y’)) = PM(x, x’) PM(y, y’)
とは限らない.
coupling method: 上記のような性質を満たす( M の)
coupling の, 吸収的な同値類 Ω*=｛(x, x) | x∈Ω｝に吸
収されるまでの反復回数を算定して, Markov chain M
の定常分布への収束に必要な反復回数（mixing
time） を算定する (詳細略).
35
coupling lemma (Aldous)
coupling lemma [Aldous 1983]
M: Ω 上の Markov chain
d(x, y)： x から yへの距離（非負整数）
D：Ωの直径（最遠状態対間の距離)
∃coupling : 0<∃β< 1, ∀（x, y)∈ Ω2,
（x, y) →（X’, Y’), E[d(X’, Y’)] ≦βd(x, y)
このとき以下が成り立つ
τ（ε) ≦ （1-β）-1 ln(Dε-1)
証明：略（後で path coupling 略証するので, その特殊ケース）
36
2. path coupling method
[Bubley and Dyer 1997]
Markov chain の mixing time を算定する手法
coupling method [Aldous 1983]の改訂版
●アルゴリズムの分野では, Markov chain の mixing time
を算定する基本的な方法として定着しつつある
37
length on a directed graph
y
x
Ω上に距離を定義する.
有向グラフG=（Ω, A )：
Ω：頂点集合, A：有向枝集合
∀（x, y)∈Ω2， d(x, y)： x から yへの距離
（G上の有向最短路の距離（距離は枝数))
∀ （x,y) ∈Ω2, d(x,y)=0 ⇔ x = y
0<∃β< 1, ∀（x, y)∈A,
（x, y） →（X’, Y’）, E[d(X’, Y’)] ≦βd(x, y)
（注：（X’, Y’)∈Aである必要は無い.)
y
x
y
Y’
x
X’
＝
β
X’ =Y’
38
convergence ratio
推移によって，距離が0に収束する.
0<∃β< 1, ∀（x, y)∈A,
（x, y) →（X’, Y’), E[d(X’, Y’)] ≦βd(x, y)
x
y
βd(x, y) X’
Y’
d(x, y)
z D
＝
β
Z’
β2d(x, y)
ε
Dβt ＝ε
t ＝(ln β-1)-1 ln (Dε-1 )
≦(1-β)-1 ln (Dε-1 )
（β ≒１）
39
path coupling theorem (Bubley and Dyer)
Theorem [Bubley and Dyer 1997]
M: Ω 上の Markov chain
G=(Ω, A )： 有向グラフ（完全とは限らない）
d(x, y)： G上での x から yへの距離(非負整数）
D：Gの直径（最遠頂点対間の距離)
∃coupling : 0<∃β< 1, ∀（x, y)∈A,
（x, y) →（X’, Y’), E[d(X’, Y’)] ≦βd(x, y)
このとき以下が成り立つ
＝
β
τ（ε) ≦ （1－β）－1 ln(Dε－1)
τ（ε)：mixing time (次ページで定義)
40
mixing time
π：Ω→[0,1]：Markov chain Mの定常分布
π’：Ω→[0,1] ：Ω上の(任意の)確率分布
DTV(π,π’)=(1/2)∑｛|π(x)-π’（x) |：x∈Ω｝
total variation
Px(t)(y)：Markov chain M が, 初期状態 x から
t 回目の推移で y に到達する確率.
τx（ε）= min｛t ：∀s≧t, DTV（π, P(s)）≦ε｝
x
τ（ε)= max ｛ τx（ε） ：x∈Ω｝：mixing time
41
distance
G=(Ω, A )： (x, y)∈A
coupling: 0<∃β< 1, ∀（x, y)∈A,
（x, y) →（X’, Y’), E[d(X’, Y’)] ≦βd(x, y)
（注：（X’, Y’)∈Aである必要は無い.)
y
x
y
Y’
x
X’ =Y’
X’
任意の状態対（x, y) から出発して, coupling が
推移すると, 状態対の距離は0に収束する.
すなわち, Ω*=｛(x, x) | x∈Ω｝ 中に収束する.
42
Proof (1/2)
補題： DTV(π,π’)≡(1/2)∑｛|π(x)－π’（x) | ：x∈Ω｝
＝max{π（A）－π’（A） ： A⊆Ω}
ただし
π（A）＝∑{π(x)：x∈A}, π’（A）＝∑{π’(x)：x∈A}.
X(t)：初期状態ｘからMarkov chain Mで推移した際の時点ｔ
での状態
Y(t)：初期状態をMの定常分布で選び, 推移した際の時点ｔ
での状態（実は定常分布に従う事は明らか）
時刻 t =（1－β）-1 ln(Dε-1)において,
∀ A⊆Ω, π(A) － Pr[X(t) ∈A] ≦εを示せばよい
（Pr[Y(t) ∈A]－ Pr[X(t) ∈A] ≦εでもよい）
43
Proof （2/2）
∀ A⊆Ω, π(A) － Pr[X(t) ∈A] ≦εを示す.
Pr[X(t) ∈A]≧Pr[Y(t) ∈A∧X(t)＝Y(t) ]
≧Pr[Y(t) ∈A]－Pr[X(t)≠ Y(t) ]
Y(t)の定義
＝π（A）－Pr[X(t)≠ Y(t) ]
距離の整数性
毎回
β倍
≧ π（A）－∑ (x’,y’) d(x’,y’)Pr[（X(t),Y(t))=(x’,y’)]
グラフの直径
＝π（A） －E[d(X(t), Y(t) )]
≧ π（A）－d(x,Y(0))βｔ ≧ π（A）－Dβｔ
≧ π（A）－De－ｔ (1－β)＝ π（A）－De－ln(D / ε)
ln β≦β－1
＝ π（A）－ε
－1 ln(Dε-1)
ln βt≦t (β－1)
時刻
t
=（1－
β
）
t
－t (1－β)
β ≦e
44
1. mixing time
実際に 有向グラフと coupling を作る.
45
directed graph
G=(Ω, A )： (x, y)∈A⇔2列での±1で推移可能
２
６
x５
６
５
３
２
４
１
６
７
４
４
３
１
３
５
２
０
８
９
０
３
８
x
y
０ ‐１
０ +１
０ +１
０ -１
０ +１
０ -１
０ -１
０ +１
２
y ６
５
６
０
０
０
０
４
４
３
３
０
０
０
０
１
６
７
４
５
２
０
４
５
２
０
８
９
０
３
８
46
coupling example (case 0)
M1のcoupling ：
x
x1
２ ５ １ ４ ５ ９
６ ３ ６ ３ ２ ０
３ ５ ０ ４ ５ ９
５ ３ ７ ３ ２ ０
２ ５ １ ４ ５ ９
６ ３ ６ ３ ２ ０
４ ２ ８ １ ０ ３
７ ４ ３ ３ ８ ８
５ ２ ７ １ ０ ３
６ ４ ４ ３ ８ ８
1/4
５ ２ ７ １ ０ ３
６ ４ ４ ３ ８ ８
０ ５ ３ ４ ５ ９
８ ３ ４ ３ ２ ０
１ ５ ２ ４ ５ ９
７ ３ ５ ３ ２ ０
x
case 0: d(x,y)=0 2
７ ２ ５ １ ０ ３
４ ４ ６ ３ ８ ８
６ ２ ６ １ ０ ３
５ ４ ５ ３ ８ ８
x x 1 x2 x3 x4
すなわち x=y のとき.
(x,x)
x3
¼
¼
(x1,x1) (x2,x2)
¼
¼
(x3,x3) (x4,x4)
x
x1
x2
x3
x4
¼
x4
47
coupling (case 0)
G=(Ω, A )： (x, y)∈A⇔M1で推移可能
case 0 : d(x, y)=0 (x＝y)
case 1 : d(x, y)=1 (x≠y) xとyは2列異なっている
2列をｘとｙが同じ列から選ぶ
case 2 : d(x, y)=1 (x≠y)
2列をｘとｙが異なる2列を選ぶ
case 3 : d(x, y)=1 (x≠y)
2列を、
1列はｘとｙが異なる列から、
1列はｘとｙが同じ列から選ぶ
２５１４５９
６３６３２０
５２７１０３ x
６４４３８８
１６１４５９
７２６３２０
y
６１７１０３
５５４３８８
48
coupling (case 1)
G=(Ω, A )： (x, y)∈A⇔M1で推移可能
case 0 : d(x, y)=0 (x＝y)
case 1 : d(x, y)=1 (x≠y) xとyは2列異なっている
2列をｘとｙが同じ列から選ぶ
case 2 : d(x, y)=1 (x≠y)
2列をｘとｙが異なる2列を選ぶ
case 3 : d(x, y)=1 (x≠y)
2列を、
1列はｘとｙが異なる列から、
1列はｘとｙが同じ列から選ぶ
２５１４５９
６３６３２０
５２７１０３ x
６４４３８８
j’ j’’
１６１４５９
７２６３２０
y
６１７１０３
５５４３８８
{ j’,j’’}∩{1,2}=φ
coupling example (case 1)
d(x,y)=1
２５１４５９
６３６３２０
５２７１０３ x
６４４３８８
j’ j’’
１６１４５９
７２６３２０
y
６１７１０３
５５４３８８
x1
½
49
½
２５１５４９
６３６２３０
５２７０１３
６４４４７８
y1
½ ½
１６１５４９
７２６２３０
６１７０１３
５５４４７８
x2
２５１４５９
６３６３２０
５２７１０３
６４４３８８
y2
１６１４５９
７２６３２０
６１７１０３
５５４３８８
(x2, y2)
(x1, y1) ½ ½
(x, y)
d(x, y) = d(x1, y1) = d(x2, y2) = 1 (距離は変わらない)
50
{ j’,j’’}∩{1,2}=φ
coupling example (case 1)
d(x,y)=1
½
½
(x
,
y
)
1 1
(x, y)
x x1 x2 y y1 y2
x
x1
x2
y
y1
y2
½
½
(x2, y2)
51
coupling (case 2)
G=(Ω, A )： (x, y)∈A⇔M1で推移可能
case 0 : d(x, y)=0 (x＝y)
case 1 : d(x, y)=1 (x≠y) xとyは2列異なっている
2列をｘとｙが同じ列から選ぶ
case 2 : d(x, y)=1 (x≠y)
2列をｘとｙが異なる2列を選ぶ
case 3 : d(x, y)=1 (x≠y)
2列を、
1列はｘとｙが異なる列から、
1列はｘとｙが同じ列から選ぶ
２５１４５９
６３６３２０
５２７１０３ x
６４４３８８
j’ j’’
１６１４５９
７２６３２０
y
６１７１０３
５５４３８８
coupling example (case 2)
d(x,y)=1
２５１４５９
６３６３２０
５２７１０３ x
６４４３８８
j’
j’’
１６１４５９
７２６３２０
y
６１７１０３
５５４３８８
{ j’,j’’}={1,2}
1/7
52
1/7
x1
x7
６１１４５９
０７１４５９
２７６３２０
８１６３２０
１ ６ ７ １ ０ ３ ‥‥ ７ ０ ７ １ ０ ３
10０ ４ ３ ８ ８
４６４３８８
y1
1/7 1/7
y7
６１１４５９
０７１４５９
２７６３２０
８１６３２０
１ ６ ７ １ ０ ３ ‥‥ ７ ０ ７ １ ０ ３
10０ ４ ３ ８ ８
４６４３８８
1/7
1/7
(x7, y7)
(x1, y1)
(x, y)
d(x1, y1) =‥= d(x7, y7) = 0 (短くなる!)
coupling example (case 2)
{ j’,j’’}={1,2
}
53
1/7
d(x,y)=1
(x1, y1)
(x, y)
x x1 x2 … y y1 y2
x
x1
x2
…
y
y1
y2
…
1/7
(x7, y7)
…
…
1/7
54
coupling (case 3)
G=(Ω, A )： (x, y)∈A⇔M0で推移可能
case 0 : d(x, y)=0 (x＝y)
case 1 : d(x, y)=1 (x≠y) xとyは2列異なっている
2列をｘとｙが同じ列から選ぶ
case 2 : d(x, y)=1 (x≠y)
2列をｘとｙが異なる2列を選ぶ
case 3 : d(x, y)=1 (x≠y)
2列を、
1列はｘとｙが異なる列から、
1列はｘとｙが同じ列から選ぶ
２５１４５９
６３６３２０
５２７１０３ x
６４４３８８
j’
j’’
１６１４５９
７２６３２０
y
６１７１０３
５５４３８８
55
|{
j’,j’’}∩{1,2}|=1
coupling example (case 3)
1/4
1/4
1/4
1/4
３５０４５９ ２５１４５９ １５２４５９ ０５３４５９
５３７３２０ ６３６３２０ ７３５３２０ ８３４３２０
４２８１０３ ５２７１０３ ６２６１０３ ７２５１０３
７４３３８８ ６４４３８８ ５４５３８８ ４４６３８８
j’
２６０４５９
６２７３２０
５１８１０３
６５３３８８
j’’
x
上下段の個数は最大1しか違わない
y
１６１４５９
７２６３２０
６１７１０３
５５４３８８
1/3
1/3
０６２４５９
８２５３２０
７１６１０３
４５５３８８
1/3
|{ j’,j’’}∩{1,2}|=1
1/4
1/4
1/4
1/4
２５１４５９ １５２４５９
６３６３２０ ７３５３２０
５２７１０３ ６２６１０３
６４４３８８ ５４５３８８
56
coupling example (case 3)
３５０４５９
５３７３２０
４２８１０３
７４３３８８
j’
２６０４５９
６２７３２０
５１８１０３
６５３３８８
j’’
１６１４５９
７２６３２０
６１７１０３
５５４３８８
1/3
1/3
０５３４５９
８３４３２０
７２５１０３
４４６３８８
距離=1
０６２４５９
８２５３２０
７１６１０３
４５５３８８
1/3
coupling example (case 3)
1/4
1/4
x1
３５０４５９
５３７３２０
４２８１０３
７４３３８８
x=x2
２５１４５９
６３６３２０
５２７１０３
６４４３８８
|{ j’,j’’}∩{1,2}|=1
1/4
x4
1/4
x3
１５２４５９ ０５３４５９
７３５３２０ ８３４３２０
６２６１０３ ７２５１０３
５４５３８８ ４４６３８８
(x1,y1)
(x3,y3 )
(x2,y2 )
(x4,y3 )
(x3,y2 )
(x2,y1 )
(x,y)
２６０４５９ １６１４５９ ０６２４５９
６２７３２０ ７２６３２０ ８２５３２０
５１８１０３ ６１７１０３ ７１６１０３
６５３３８８ ５５４３８８ ４５５３８８
y1
y=y2
y3
1/3
1/3
1/3
距離=1
57
coupling example (case 3)
1/4
1/4
x1
３５０４５９
５３７３２０
４２８１０３
７４３３８８
x=x2
２５１４５９
６３６３２０
５２７１０３
６４４３８８
|{ j’,j’’}∩{1,2}|=1
1/4
x4
1/4
x3
１５２４５９ ０５３４５９
７３５３２０ ８３４３２０
６２６１０３ ７２５１０３
５４５３８８ ４４６３８８
２６０４５９
６２７３２０
５１８１０３
６５３３８８
2/12
2/12 1/12 3/12
(x,y)
１６１４５９ ０６２４５９
７２６３２０ ８２５３２０
６１７１０３ ７１６１０３
５５４３８８ ４５５３８８
y1
y=y2
3/12
1/12
y3
1/3
1/3
1/3
距離=1
58
59
coupling example (case 3)
|{ j’,j’’}∩{1,2}|=1
1/3
1/3
1/3
y1 y=y2 y3
1/4
1/4
x1
x=x2
3/12
1/12
2/12
2/12
1/4
1/4
x3
1/12
3/12
x4
d(x1, y1) =‥= d(x4, y3) = 1 (距離は変わらない)
60
coupling Theorem
x と y が1,2列で異なるとき, ｛j’,j’’｝として｛1, 2｝
を選ぶならば, 距離は0になり, 他の場合は距
離は変わらない（1のままである).
(case 1)
２５１４５９
６３６３２０
x
５２７１０３
６４４３８８
j’ j’’
１６１４５９
y ７２６３２０
６１７１０３
５５４３８８
距離不変(=1)
(case 2)
２５１４５９
６３６３２０
５２７１０３
６４４３８８
j’ j’’
１６１４５９
７２６３２０
６１７１０３
５５４３８８
距離=０
(case 3)
２５１４５９
６３６３２０
５２７１０３
６４４３８８
j’
j’’
１６１４５９
７２６３２０
６１７１０３
５５４３８８
距離不変(=1)
61
coupling Theorem
x と y が1,2列で異なるとき, ｛j’,j’’｝として｛1, 2｝
を選ぶならば, 距離は0になり, 他の場合は距
離は変わらない（1のままである).
｛j’,j’’｝として｛1, 2｝を選ぶ確率は 2/(n(n-1)).
推移後の距離の期待値:
1×（1- 2/(n(n-1)))＋0× 2/(n(n-1))
=1- 2/(n(n-1))
62
directed graph （diameter)
G=(Ω, A )： (x, y)∈A⇔M0 で推移可能
(x, y)∈A ならば, x と y は 2m×2のセルにおいて,
±1違う.
ゆえに, N/2m回以内の推移で, 任意の状態間
を移ることができる.
D=Gの直径 ≦ N/2m = dn
d = N/(n2m ) : セル内の数値の平均
63
main theorem
Theorem
M1: new Markov chain (Ω : 状態空間）
G=(Ω, A )： 有向グラフ
d(x, y)： G上での x から yへの距離
D：Gの直径≦ N/2m = dn
d=N/(n2m ) ：セル内の数値の平均
∃coupling : 0<β< 1, ∀（x, y)∈A,
（x, y) →（X’, Y’), E[d(X’, Y’)] ≦βd(x, y)
＝
β=1- 2/(n(n-1))
β
このとき以下が成り立つ
τ（ε) ≦ (1/2) n(n-1) ln(dnε－1)
64
main theorem
Theorem
M1: new Markov chain (Ω : 状態空間）
G=(Ω, A )： 有向グラフ
d(x, y)： G上での x から yへの距離
D：Gの直径≦ N/2m = dn
d=N/(n2m ) ：セル内の数値の平均
0< ∃ β< 1, ∀（x, y)∈A,
（x, y) →（X’, Y’), E[d(X’, Y’)] ≦βd(x, y)
＝
β=1- 2/(n(n-1))
β
このとき以下が成り立つ
τ（ε) ≦ (1/2) n(n-1) ln(dnε－1)
65
number of iterations
実際の反復回数は?
発表での数値例 ε=10-10
τ（ε) ≦ (1/2) n(n-1) ln((N/2m) /ε)
=(1/2)6(6-1)ln((97/26 ) 1010 )=355.7
セルの数値の平均値を10程度とする
d = N/(n2m ) = 10 とする
τ（ε) ≦ (1/2) n(n-1) ln(dn/ε)
= (1/2) n(n-1) ln (10n/ε)
単位千回
n
5 10 20 40 50 100
ε10-10 0.3 1.2 5.3 22 35 148
10-12 0.3 1.5 6.3 26 42 171
66
0. 超幾何分布を定常分布に持つ
Markov chain の提案
超幾何分布のためのMarkov chain
67
超幾何分布
３元分割表の集合上の超幾何分布
Pr[T]＝Pr[nijk]＝
Δ
Πi Πj Πknijk !
Δ：総和が1となるように定めた係数
68
new Markov chain
M1: 列{j’, j’’} （j’≠j’’)をランダムに選ぶ
３
５
４
７
２
６
５
６
５
３
２
４
１
６
７
４
５
３
２
４
０
７
８
３
４
３
１
３
２
６
５
６
４
３
１
３
５
２
０
８
５
３
２
４
５
２
０
８
９
０
３
８
１
６
７
４
９
０
３
８
-θ
+θ
+θ
-θ
-1
0
４
３
１
３
５
２
０
８
９
０
３
８
０ +θ ０
０ -θ ０
０ -θ ０
０ +θ ０
０
０
０
０
０
０
０
０
０
８
７
４
４
３
１
３
５
３
２
４
５
２
０
８
３
４
５
６
９
０
３
８
1
2
１
７
６
５
５
３
２
４
２
５
６
５
４
３
１
３
５
２
０
８
９
０
３
８
69
new Markov chain
M1: 列{j’, j’’} （j’≠j’’)をランダムに選ぶ
Δ’/(0!3!8!4!7!5!4!6!)
Δ’/(3!0!5!7!4!8!7!3!)
Δ’：総和が1
となるよう
定めた係数
Δ’/(3!0!5!7!
4!8!7!3!)
３
５
４
７
５
３
２
４
０
７
８
３
４
３
１
３
２
６
５
６
５
２
０
８
５
３
２
４
９ Δ’/(2!1!6!6!
０
5!7!6!4!)
３
８
１
６
７
４
４
３
１
３
５
２
０
８
９
０
３
８
１
７
６
５
５
３
２
４
２
５
６
５
０
８
７
４
４
３
１
３
５
３
２
４
５
２
０
８
３
４
５
６
９
０
３
８
４
３
１
３
５
２
０
８
９
０
３
８
70
Mixing time
一様分布と同様にグラフを定め, coupling を作
ることができる．
Theorem
τ（ε) ≦ (1/2) n(n-1) ln(dnε－1)
証明のキー：
coupling を作る際, 一様分布と違って, 推移確
率を陽に書けない. 推移確率の存在性を証
明する必要がある.
coupling example (case 3)
Δ’/2!1!6!6!5!7!6!4!
Δ’/3!0!5!7!4!8!7!3!
３５０４５９
５３７３２０
４２８１０３
７４３３８８
?
?
２６０４５９
６２７３２０
５１８１０３
６５３３８８
Δ’’/2!0!6!7!
5!8!6!3!
２５１４５９
６３６３２０
５２７１０３
６４４３８８
?
(x,y)
１６１４５９
７２６３２０
６１７１０３
５５４３８８
Δ’/1!2!7!5!
6!6!5!5!
１５２４５９
７３５３２０
６２６１０３
５４５３８８
?
?
Δ’/0!3!8!4!
7!5!4!6!
71
０５３４５９
８３４３２０
７２５１０３
４４６３８８
?
値は唯一なので
０ ６ ２ ４ ５ ９ 非負性の証明
８ ２ ５ ３ ２ ０ が重要となる
７１６１０３
４５５３８８
Δ’’/0!2!8!5!
7!6!4!5!
Δ’’/1!1!7!6!6!7!5!4!
Δ’/2!1!6!6!5!7!6!4!
coupling example
Δ’/3!0!5!7!4!8!7!3!
３５０４５９
５３７３２０
４２８１０３
７４３３８８
?
２６０４５９
６２７３２０
５１８１０３
６５３３８８
Δ’’/2!0!6!7!
5!8!6!3!
２５１４５９
６３６３２０
５２７１０３
６４４３８８
(x,y)
１６１４５９
７２６３２０
６１７１０３
５５４３８８
Δ’/1!2!7!5!
6!6!5!5!
１５２４５９
７３５３２０
６２６１０３
５４５３８８
?
Δ’/0!3!8!4!
7!5!4!6!
72
０５３４５９
８３４３２０
７２５１０３
４４６３８８
?
大小関係が
逆転するのは
最大1回
０６２４５９
８２５３２０
７１６１０３
４５５３８８
Δ’’/0!2!8!5!
7!6!4!5!
Δ’’/1!1!7!6!6!7!5!4!
Δ’/2!1!6!6!5!7!6!4!
coupling example
Δ’/3!0!5!7!4!8!7!3!
３５０４５９
５３７３２０
４２８１０３
７４３３８８
?
２６０４５９
６２７３２０
５１８１０３
６５３３８８
Δ’’/2!0!6!7!
5!8!6!3!
２５１４５９
６３６３２０
５２７１０３
６４４３８８
Δ’/1!2!7!5!
6!6!5!5!
１５２４５９
７３５３２０
６２６１０３
５４５３８８
非負
(x,y)
１６１４５９
７２６３２０
６１７１０３
５５４３８８
Δ’/0!3!8!4!
7!5!4!6!
73
０５３４５９
８３４３２０
７２５１０３
４４６３８８
?
大小関係が
逆転するのは
最大1回
０６２４５９
８２５３２０
７１６１０３
４５５３８８
Δ’’/0!2!8!5!
7!6!4!5!
Δ’’/1!1!7!6!6!7!5!4!
74
coupling
非負性の証明
ここの非負性は自明
Ａ
Ｃ
Ｅ
逆転
Ｂ
Ｄ
Ｂ－Ａ≧0
Ｃ＋Ｅ－Ｄ
≧Ｃ－Ｄ≧0
75
Concluding
2×2×・・・×2×Kの分割表をサンプリングす
る Markov chain として,極限分布に一様分布
と超幾何分布を持つものを提案し, その
mixing time を算定した.
Theorem
τ（ε) ≦ （1/2） n(n-1) ln(dnε－1)
τ(ε): mixing time
d=N/(n2m ) ：セル内の数値の平均
n, m : 分割表のサイズは{1,2}m×{1,2,...,n}
N：分割表中の数値の総和
76
終わり
77
future work
Counting （D&G の拡張で可能と思われる）
計算機実験
78
main result
Purpose
random generation of ｛1,2｝m×｛1,2,...,n｝
contingency tabel by using Markov chain
Monte Calro method
Result
The mixing time of the Markov chain is
bounded by the polynomial of # of columns
and the logarithm of table sum.
Method
path coupling method [Bubley and Dyer 1997]
79
contingency table
{1,2}2×{1,2,3,4,5,6}の例
80
contingency table
２
６
８
５
６
１１
５ １
３ ６
８ ７
２ ７
４ ４
６ １１
４
３
７
１
３
４
５ ９
２ ０
７ ９
０ ３
８ ８
８ １１
２６
２０
１８
２７
７ ７ ８ ５ ５ １１
１２ ７ １０ ６ １０ ８
N：table sum
(N=91)
n: # of columns (n=6)
81
main theorem
Theorem
τ(ε)≦n(n-1) ln (([N/2m]+1)εｰ1）/2
ε：error bound
(precise definion → later)
τ(ε): mixing time
(# of iterations of MCMC method)
(precise definion → later)
N：table sum
n, m : the size of table is {1,2}m×{1,2,...,n}
82
coupling (trivial coupling)
trivial coupling:
P((x, y),(x’, y’)) = PM(x, x’) PM(y, y’)
∀ x ,∀ y∈Ω,
PM(x, x’)=∑｛P((x, y),(x’, y’)) | y’∈Ω｝
PM(y, y’)=∑｛P((x, y),(x’, y’)) | x’∈Ω｝
P((x, y),(x’, y’)) ：coupling における (x, y)から(x’, y’)
への推移確率
PM(x, x’)：Mにおける x から x’への推移確率
© Copyright 2024 ExpyDoc
