1 / 82

Random Generation of B m × J Contingency Tables ver. 5.0

Random Generation of B m × J Contingency Tables ver. 5.0. 東京大学 松井知己 東海大学 松井泰子 東京理科大学 小野陽子. Random Generation of {1 ,2 } m × { 1,2,..., n } Contingency Tables. 東京大学 松井知己 東海大学 松井泰子 東京理科大学 小野陽子. Path Coupling 法を用いた 多元分割表生成のためのマルコフ連鎖設計. 東京大学 松井知己 東海大学 松井泰子 東京理科大学 小野陽子. 本研究の目的.

Download Presentation

Random Generation of B m × J Contingency Tables ver. 5.0

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript


  1. Random Generation of Bm×J Contingency Tablesver. 5.0 東京大学 松井知己 東海大学 松井泰子 東京理科大学 小野陽子

  2. Random Generation of {1,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. 列 (column) contingency table • {1,2}2×{1,2,3,4,5,6}の例

  7. 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

  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

  9. 分布関数 • 周辺和を固定したとき、 • T =( nijk)という分割表が出現する確率が、 • (1)一様分布すると仮定 ←Diaconis-Efron検定 • Pr[T]=1/ |Ω| • Ω:与えられた表と周辺和が等しい分割表の集合. • (2)超幾何分布すると仮定 ←対数線形モデル • Δ • 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元分割表 • 行と列の独立性を検定したい. • total: 周辺和

  12. 行と列の独立性 • 行と列の独立性を検定したい. H0: pij =pi・p・j pij

  13. 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

  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元分割表の最近の研究(主にアルゴリズム分野) • 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: セル内の数値の総和

  17. 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の中の更に特殊ケース)

  18. 4. 一様分布を定常分布に持つMarkov chain の提案 an extension of Dyer and Greenhill’s chain[2000]

  19. 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 推移できないときは, 列選択から再試行

  20. -θ 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

  21. -θ 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

  22. x1 x2 … x6 x7 x8 … x13 … x18 x19 … x24 geometric interpretation • 分割表の集合 • =24次元空間中の多面体中の整数格子点 x24 … x4 x3 x2 x1 O

  23. 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

  24. θ 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

  25. properties of new Markov chain • 提案する Markov chain M1の特徴 • Markov chain M1は, エルゴード的. • finite • a-periodic • 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. Markov chainの無記憶性 →yから出発したsample path と同じに見える! coupling sample path • 2つのサンプルパスから, 新しいサンプルパスを作る. 状態が最初に一致した時点で 他方のsample path に乗り換える 状態空間 x y 0 時刻

  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 chainM を 異なる初期状態 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 couplingmethod[Bubley and Dyer 1997] Markov chain の mixing time を算定する手法 coupling method[Aldous 1983]の改訂版 ●アルゴリズムの分野では, Markov chain の mixing time を算定する基本的な方法として定着しつつある

  37. 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である必要は無い.)

  38. β 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)

  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 • P(t)(y):Markov chain Mが, 初期状態 x から • x • t 回目の推移で y に到達する確率. • τx(ε)= min{t:∀s≧t, DTV(π, P(s))≦ε} • x • τ(ε)= max { τx(ε) :x∈Ω}:mixing time

  41. 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∈Ω} 中に収束する.

  42. 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] ≦εでもよい)

  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) ] • =π(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)

  44. 1. mixing time 実際に 有向グラフと coupling を作る.

  45. 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

  46. 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

  47. 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が同じ列から選ぶ

  48. 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が同じ列から選ぶ

  49. 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 (距離は変わらない)

  50. 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 ½ ½

More Related