820 likes | 911 Views
Random Generation of B m × J Contingency Tables ver. 5.0. 東京大学 松井知己 東海大学 松井泰子 東京理科大学 小野陽子. Random Generation of {1 ,2 } m × { 1,2,..., n } Contingency Tables. 東京大学 松井知己 東海大学 松井泰子 東京理科大学 小野陽子. Path Coupling 法を用いた 多元分割表生成のためのマルコフ連鎖設計. 東京大学 松井知己 東海大学 松井泰子 東京理科大学 小野陽子. 本研究の目的.
E N D
Random Generation of Bm×J Contingency Tablesver. 5.0 東京大学 松井知己 東海大学 松井泰子 東京理科大学 小野陽子
Random Generation of {1,2}m×{1,2,...,n}Contingency Tables 東京大学 松井知己 東海大学 松井泰子 東京理科大学 小野陽子
Path Coupling 法を用いた多元分割表生成のためのマルコフ連鎖設計 東京大学 松井知己 東海大学 松井泰子 東京理科大学 小野陽子
本研究の目的 • 題目:Path Coupling 法を用いた多元分割表生成のためのマルコフ連鎖設計 • 扱う分割表の形状: 2×2×・・・×2×K • 2×2×K:[太田,青木,広津(昨日)] • マルコフ基底:自明なものが存在 • 3×3×K:[青木,竹村(昨日)] • 極限分布: (1)一様分布 (2)超幾何分布 • 目的:1つサンプリングする迄に必要な • 推移回数の算定 =mixing time の算定
MCMC法 • 多元分割表の行と列の独立性の検定 • 検定理論 • 正確検定 • 医療統計 • p値の計算 • 帰無仮説の棄却を決める閾値 • MCMC法(Markov Chain Monte Carlo Method) • Markov chain で sampling を行い • Monte Carlo 法を行う • mixing time の速いMarkov chain の必要性
列 (column) contingency table • {1,2}2×{1,2,3,4,5,6}の例
contingency table (example) 26 20 • 2 5 1 4 5 9 • 6 3 6 3 2 0 • 5 2 7 1 0 3 • 6 4 4 3 8 8 • N:table sum (N=97=26+20+18+33) • n: # of columns (n=6) 8 8 7 7 7 9 18 33 11 6 11 4 8 11 7 7 8 5 5 12 12 7 10 6 10 8
contingency tables 26 20 • * * * * * * • * * * * * * • * * * * * * • * * * * * * • N:table sum (N=97=26+20+18+33) • n: # of columns (n=6) 8 8 7 7 7 9 18 33 11 6 11 4 8 11 7 7 8 5 5 12 12 7 10 6 10 8
分布関数 • 周辺和を固定したとき、 • T =( nijk)という分割表が出現する確率が、 • (1)一様分布すると仮定 ←Diaconis-Efron検定 • Pr[T]=1/ |Ω| • Ω:与えられた表と周辺和が等しい分割表の集合. • (2)超幾何分布すると仮定 ←対数線形モデル • Δ • Pr[T]=Pr[nijk]= • Πi ΠjΠk nijk ! • Δ:総和が1となる係数
main result • 目的 • {1,2}m×{1,2,...,n}分割表を ランダム生成する (定常分布が一様分布と超幾何分布の) • Markov chain の提案 • 結果 • 提案する Markov chain の mixing time は, • ・分割表の列数, • ・ln(分割表のセル中の数値の平均), • ・ln (精度の逆数) • の多項式でおさえられる. • 方法 • path coupling 法 [Bubley and Dyer 1997]
2元分割表 (contingency table) • 2元分割表 • 行と列の独立性を検定したい. • total: 周辺和
行と列の独立性 • 行と列の独立性を検定したい. H0: pij =pi・p・j pij
eye hair Blue Brown Green Blond 3 4 3 10 Brown 2 4 2 8 Red 4 3 5 12 9 11 10 30 eye hair Blue Brown Green Blond 9・10 /30 11・10 /30 10・10 /30 10 Brown 9・8 /30 11・8 /30 10・8 /30 8 Red 9・12 /30 11・12 /30 10・12 /30 12 9 11 10 30 行と列の独立性の検定 • 行と列の独立性を検定したい. ni+ nij n+j mij
正確検定 • 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∈Ω} の値を計算する.
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 ):セル内の数値の平均
2元分割表の最近の研究(主にアルゴリズム分野) • 2元分割表(m×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 可能 N: セル内の数値の総和
3元分割表の最近の研究(主にアルゴリズム分野)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部グラフの完全マッチングのサンプリング • (2×I×J ,周辺和は0or1の中の更に特殊ケース)
4. 一様分布を定常分布に持つMarkov chain の提案 an extension of Dyer and Greenhill’s chain[2000]
2 5 1 4 5 9 6 3 6 3 2 0 5 2 7 1 0 3 6 4 4 3 8 8 0 ‐1 0 +1 0 0 0 +1 0 -1 0 0 0 +1 0 -1 0 0 0 -1 0 +1 0 0 0 +1 0 -1 0 0 0 -1 0 +1 0 0 0 -1 0 +1 0 0 0 +1 0 -1 0 0 2 4 1 5 5 9 6 4 6 2 2 0 2 6 1 3 5 9 6 2 6 4 2 0 5 3 7 0 0 3 6 3 4 4 8 8 5 1 7 2 0 3 6 5 4 2 8 8 Markov chain M0 (extension of Diaconis and Saloff-Coste’s chain) • M0: 列{j’, j’’} (j’≠j’’)をランダムに選ぶ 1/2 1/2 推移できないときは, 列選択から再試行
-θ 0 +θ0 0 0 +θ 0 -θ0 0 0 2 5 1 4 5 9 6 3 6 3 2 0 +θ0 -θ0 0 0 -θ 0 +θ0 0 0 5 2 7 1 0 3 6 4 4 3 8 8 3 5 0 4 5 9 5 3 7 3 2 0 2 5 1 4 5 9 6 3 6 3 2 0 0 5 3 4 5 9 8 3 4 3 2 0 1 5 2 4 5 9 7 3 5 3 2 0 6 2 6 1 0 3 5 4 5 3 8 8 4 2 8 1 0 3 7 4 3 3 8 8 7 2 5 1 0 3 4 4 6 3 8 8 5 2 7 1 0 3 6 4 4 3 8 8 new Markov chain M1 (extension of Dyer and Greenhill’s chain) • M1: 列{j’, j’’} (j’≠j’’)をランダムに選ぶ θ= 2 -1 1 0
-θ 0 +θ0 0 0 +θ 0 -θ0 0 0 2 5 1 4 5 9 6 3 6 3 2 0 +θ0 -θ0 0 0 -θ 0 +θ0 0 0 5 2 7 1 0 3 6 4 4 3 8 8 2 5 1 4 5 9 6 3 6 3 2 0 3 5 0 4 5 9 5 3 7 3 2 0 0 5 3 4 5 9 8 3 4 3 2 0 1 5 2 4 5 9 7 3 5 3 2 0 6 2 6 1 0 3 5 4 5 3 8 8 5 2 7 1 0 3 6 4 4 3 8 8 7 2 5 1 0 3 4 4 6 3 8 8 4 2 8 1 0 3 7 4 3 3 8 8 new Markov chain M1 • M1: 列{j’, j’’} (j’≠j’’)をランダムに選ぶ 1/4 1/4 1/4 1/4
x1 x2 … x6 x7 x8 … x13 … x18 x19 … x24 geometric interpretation • 分割表の集合 • =24次元空間中の多面体中の整数格子点 x24 … x4 x3 x2 x1 O
1/2 1/2 geometric interpretation of M0 • new Markov chain M0 j’ j’’ 2 5 1 4 5 9 6 3 6 3 2 0 5 2 7 1 0 3 6 4 4 3 8 8 x24 … x4 x3 x2 x1 O
θ 1/4 1/4 1/4 1/4 geometric interpretation of M1 • new Markov chain M1 j’ j’’ 2 5 1 4 5 9 6 3 6 3 2 0 5 2 7 1 0 3 6 4 4 3 8 8 x24 … x4 x3 x2 x1 O
properties of new Markov chain • 提案する Markov chain M1の特徴 • Markov chain M1は, エルゴード的. • finite • a-periodic • irreducible ({1,2}m×{1,2,...,n}という性質) • 定理 • M1 の 定常分布は, 与えられた分割表全体の集合上の一様分布である. • (証明略)
3. coupling method[Aldous 1983] Markov chain の mixing time を算定する手法
coupling sample path • 2つのサンプルパスから, 新しいサンプルパスを作る. 状態空間 x y 0 時刻
Markov chainの無記憶性 →yから出発したsample path と同じに見える! coupling sample path • 2つのサンプルパスから, 新しいサンプルパスを作る. 状態が最初に一致した時点で 他方のsample path に乗り換える 状態空間 x y 0 時刻
coupling sample path • 仮想的にすべての状態からサンプルパスを発生させる. 状態空間 x 0 時刻 パスが一つになった → 初期状態に依存していない. → 定常分布
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’への推移確率
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
coupling (as a simulation) • 初期状態(x, y)からのcoupling の推移 • =Markov chainM を 異なる初期状態 xとy から推移したときのシミュレーション y∈Ω x∈Ω (x,y)∈Ω2
coupling (初期状態への依存) • 任意の初期状態 (x, y) から, 対角要素 (x’,x’) へcoupling で推移するまでの時間を算定する. • 初期状態に依存しない. y∈Ω x∈Ω x’ (x,y)∈Ω2
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) を算定する (詳細略). Ω*
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 略証するので, その特殊ケース)
2. path couplingmethod[Bubley and Dyer 1997] Markov chain の mixing time を算定する手法 coupling method[Aldous 1983]の改訂版 ●アルゴリズムの分野では, Markov chain の mixing time を算定する基本的な方法として定着しつつある
= y y β Y’ x x X’ =Y’ X’ y length on a directed graph 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である必要は無い.)
= β convergence ratio • 推移によって,距離が0に収束する. • 0<∃β< 1, ∀(x, y)∈A, • (x, y) →(X’, Y’), E[d(X’, Y’)] ≦βd(x, y) y x z D d(x, y) Y’ X’ Z’ βd(x, y) Dβt =ε t=(ln β-1)-1 ln (Dε-1 ) ≦(1-β)-1ln (Dε-1 ) β2d(x, y) ε (β ≒1)
= β 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 (次ページで定義)
mixing time • π:Ω→[0,1]:Markov chain Mの定常分布 • π’:Ω→[0,1] :Ω上の(任意の)確率分布 • DTV(π,π’)=(1/2)∑{|π(x)-π’(x) |:x∈Ω} • total variation • P(t)(y):Markov chain Mが, 初期状態 x から • x • t 回目の推移で y に到達する確率. • τx(ε)= min{t:∀s≧t, DTV(π, P(s))≦ε} • x • τ(ε)= max { τx(ε) :x∈Ω}:mixing time
y y Y’ x x X’ =Y’ X’ 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である必要は無い.) • 任意の状態対(x, y) から出発して, coupling が推移すると, 状態対の距離は0に収束する. すなわち, Ω*={(x, x) | x∈Ω} 中に収束する.
Proof (1/2) • 補題: DTV(π,π’)≡(1/2)∑{|π(x)-π’(x) | :x∈Ω} • =max{π(A)-π’(A) : A⊆Ω} • ただし • π(A)=∑{π(x):x∈A}, π’(A)=∑{π’(x):x∈A}. • X(t):初期状態xからMarkov chain Mで推移した際の時点t での状態 • Y(t):初期状態をMの定常分布で選び, 推移した際の時点t での状態(実は定常分布に従う事は明らか) • 時刻t =(1-β)-1 ln(Dε-1)において, • ∀ A⊆Ω, π(A) - Pr[X(t) ∈A] ≦εを示せばよい • (Pr[Y(t) ∈A]- Pr[X(t) ∈A] ≦εでもよい)
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) ] • =π(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))βt≧ π(A)-Dβt • ≧ π(A)-De-t(1-β)= π(A)-De-ln(D / ε) • = π(A)-ε Y(t)の定義 距離の整数性 毎回 β倍 グラフの直径 ln β≦β-1 ln βt≦t (β-1) βt≦e -t (1-β) 時刻t =(1-β)-1 ln(Dε-1)
1. mixing time 実際に 有向グラフと coupling を作る.
2 5 1 4 5 9 6 3 6 3 2 0 5 2 7 1 0 3 6 4 4 3 8 8 0 ‐1 0 +1 0 0 0 +1 0 -1 0 0 0 +1 0 -1 0 0 0 -1 0 +1 0 0 2 4 1 5 5 9 6 4 6 2 2 0 5 3 7 0 0 3 6 3 4 4 8 8 directed graph • G=(Ω, A ): (x, y)∈A⇔2列での±1で推移可能 x x y y
2 5 1 4 5 9 6 3 6 3 2 0 5 2 7 1 0 3 6 4 4 3 8 8 3 5 0 4 5 9 5 3 7 3 2 0 0 5 3 4 5 9 8 3 4 3 2 0 ¼ ¼ 4 2 8 1 0 3 7 4 3 3 8 8 7 2 5 1 0 3 4 4 6 3 8 8 ¼ ¼ 2 5 1 4 5 9 6 3 6 3 2 0 1 5 2 4 5 9 7 3 5 3 2 0 5 2 7 1 0 3 6 4 4 3 8 8 6 2 6 1 0 3 5 4 5 3 8 8 coupling example (case 0) x1 x3 • M1のcoupling : • case 0: d(x,y)=0 • すなわち x=y のとき. • (x,x) (x1,x1) (x2,x2) • (x3,x3) (x4,x4) x 1/4 x2 x4 x x1 x2 x3 x4 ¼ x x1 x2 x3 x4
x 1 6 1 4 5 9 7 2 6 3 2 0 2 5 1 4 5 9 6 3 6 3 2 0 6 1 7 1 0 3 5 5 4 3 8 8 5 2 7 1 0 3 6 4 4 3 8 8 y 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列をxとyが同じ列から選ぶ • case 2 :d(x, y)=1 (x≠y) • 2列をxとyが異なる2列を選ぶ • case 3 :d(x, y)=1 (x≠y) • 2列を、 • 1列はxとyが異なる列から、 • 1列はxとyが同じ列から選ぶ
x 1 6 1 4 5 9 7 2 6 3 2 0 2 5 1 4 5 9 6 3 6 3 2 0 j’ j’’ 6 1 7 1 0 3 5 5 4 3 8 8 5 2 7 1 0 3 6 4 4 3 8 8 y 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列をxとyが同じ列から選ぶ • case 2 : d(x, y)=1 (x≠y) • 2列をxとyが異なる2列を選ぶ • case 3 : d(x, y)=1 (x≠y) • 2列を、 • 1列はxとyが異なる列から、 • 1列はxとyが同じ列から選ぶ
1 6 1 4 5 9 7 2 6 3 2 0 1 6 1 5 4 9 7 2 6 2 3 0 1 6 1 4 5 9 7 2 6 3 2 0 2 5 1 4 5 9 6 3 6 3 2 0 2 5 1 5 4 9 6 3 6 2 3 0 2 5 1 4 5 9 6 3 6 3 2 0 5 2 7 1 0 3 6 4 4 3 8 8 6 1 7 1 0 3 5 5 4 3 8 8 6 1 7 0 1 3 5 5 4 4 7 8 6 1 7 1 0 3 5 5 4 3 8 8 5 2 7 1 0 3 6 4 4 3 8 8 5 2 7 0 1 3 6 4 4 4 7 8 coupling example (case 1) d(x,y)=1 • { j’,j’’}∩{1,2}=φ ½ ½ x1 x2 x y1 y2 j’ j’’ ½ ½ y ½ ½ (x2, y2) (x1, y1) (x, y) d(x, y) = d(x1, y1) = d(x2, y2) = 1 (距離は変わらない)
y x y1 x1 y2 x2 coupling example (case 1) d(x,y)=1 • { j’,j’’}∩{1,2}=φ ½ ½ (x2, y2) (x1, y1) (x, y) y1 y2 x2 x y x1 ½ ½