数理解析研究所講究録 1325 巻 2003 年 215-220215
分割表の近似数え上げスキームの提案$m\cross n$

東京大学大学院情報理工学系研究科数理情報学専攻
来嶋秀治(Shuji Kijima)

松井知己(Tomomi Matsui)

Department of Mathematical Informatics,
Graduate School of Information Science and Technology,
University of Tokyo, Bunkyo ku , Tokyo 113-8656, Japan概要

分割表に対する新しい近似数え上げ手法を提案する。 この手法は、 Dyerand Greenhill が 2 行分割表に対して提案した手法を改良し、拡張したものである。 我々は近似解の精度を保証するだけでなく、 得られた推定量の厳密解からの偏りの大きさについても議論
本報告では$m\mathrm{x}n$する。

2 元分割表は正の整数からなる行和と列和を持ち、 表中の値として非負整数をとる表 (行列) である。 2 元分割表は医療統計の分野などで統計データを扱うのに用いられる。 与えられた行和およ完全であること
ひ列和を溝たす 2 元分割表の個数を厳密に数える問題は、 行数が 2 の時でさえ$\#\mathrm{P}$が知られている [4]。

2000 年に Dyer and Groenhill によって MCMC (マルコフ連鎖モンテカルロ) 法を用いた 2$\mathrm{x}n$分割表に分割表の個数を求める近似解法が提案された [3]。 この方法は非常に直感的だが、適用した場合、 精度や偏りの理論的考察が困難となる。
実際、我々は Dyer and Greenhill の手法と我々の手法を実装し、 2 行分割表に対して厳密解とそれぞれの解法で得られた近似解の平均値との比較を行った。 図 1 は行和 $(42, 38)$ 、列和 $(10, \ldots, 10)$$m\mathrm{x}n$の周辺和をもつ $2\cross 8$ 分割表に対して、 それぞれの解法を 2,000 回ずつ行って得られた近似解のヒストグラムである。 この問題の厳密解は 9, 162, 736 であるが、 Dyer and Greenhill の手法による近似解の平均は大きく偏ったのに対して、我々の手法による近似解にはほとんど偏りは認めるられなかった。
本報告では、 分割表の個数に関する性質を考慮することで Dyer and Greenhill の方法を改良す分割表に拡張することが可能である。 本報告では、 我々の手法にる。 さらに我々の方法は$m\mathrm{x}n$よって得られる推定量の精度とその期待値の厳密解からの偏りの大きさについて議論する。

2近似解法
整数 (非負整数、正整数) 全体の集合をそれぞれは正整数と$\mathbb{Z}(\mathbb{Z}_{+\text{、}}\mathbb{Z}_{++})$$(r_{1}, \ldots, r_{m})\in \mathbb{Z}_{+}^{m}$} こ対して、 $\sum_{i=1}^{m}r_{t}=\sum_{j=1}^{n}Sj=N$をもち、 非負整数を表値にとる $m$ 行 列の 2 元分割表全体$s=(s_{1}, \ldots, s_{n})\in \mathbb{Z}_{+}^{m}$を満たすとする。 行和$r$

およひ列和

$s$

で表すことにする。ベクトル $r=$

$N\in \mathbb{Z}_{++}$

$n$

216

$\Xi\xi \mathrm{a}$図 1: The histogramの集合$\Sigma_{r,s}$.

$\Sigma_{r,s=}\mathrm{d}\mathrm{e}\mathrm{f}$を、.

$\{X\in \mathbb{Z}_{+}^{m\mathrm{x}n}|\sum_{j=1}^{n}X_{ij}=r_{1}$$(1\leq\forall i\leq m)$,$(1\leq\forall j\leq n)\}$$\sum_{i=1}^{m}X_{ij}=s_{j}$は行と列をは分割表 $X$ の 行 列の値を表す。明らかに分割表の個数
したがって、一般性を失うことなく $m\leq n$ とする。部分集合置き換えても変わらなで定義する。但し、$i$$X_{ij}$$|\Sigma r,s|$$j$$\Omega_{i}\subset\Sigma_{r,s}$$\mathrm{V}^{\mathrm{a}_{\text{。}}}$

を $\Omega_{i}=\{X\in\Sigma_{r,s}|X_{in}\geq\lceil s_{n}/m\rceil\}$ で定義する。 また、添え字

は $r_{k}= \max\{r_{1}, \ldots, r_{m}\}$ を満た$k$すとする。 この時、$\mathrm{d}\mathrm{e}\mathrm{f}=$$\tilde{r}$$\mathrm{d}\mathrm{e}\mathrm{f}=$$\tilde{s}$.
.

$(r_{1}, \ldots, r_{k}-\lceil s_{n}/m\rceil, \ldots, r_{m})$$(s_{1}, \ldots, s_{n-1}, \lfloor\frac{m-1}{m}s_{n}\rfloor)$$\{$$(s_{1}, \ldots, s_{n-1})$,,,$s_{n}>1$$s_{n}=1$,,がわかればが成り立つ。 もし、 $\rho=|\Omega_{k}|/|\Sigma r,s|$ とを定義する。 明らかにを求めることは困難である。 いを計算できる。 しかし、一般に およひ上で一様標本抽出が可能ならば、 の値の推定にモンテカルロ法を適用できる。すなわち
ま、$M$ 回の標本抽出を行い、 $U$ 個の分割表が \Omega \sim こ含まれるならば、 $(U+1)/(M+1)$ をの推定量と$|\Sigma|\tilde{r},\tilde{s}$$|\Omega_{k}|=|\Sigma_{\tilde{r},\tilde{s}}|$$|\Sigma r,s|=|\Sigma|\overline{r},\tilde{s}/\rho$$\Sigma_{r,s}$$\rho$$|\Sigma_{\overline{r},\overline{s}}|$$\rho$$\rho$の数え上げ問題に帰着させることが出来る。 この手続きを繰り返して、分割表のサイズが 2 行 2 列になるまで小さくする。 2を求めることができる。
分割表の個数は定数時間で求めることができるので、 再帰的にする。 こうして、 もとの問題に比べてサイズの小さな分割表$\Sigma_{\overline{r},\overline{s}}$$\mathrm{x}2$$|\Sigma r,s|$217
3

分割表の個数に関する諸定理
我々のスキームでは、 問題のサイズを小さくする目的で部分集合を定義した。 この定義には$\Omega_{k}$で二つの重要な意味がある。 一つは行の添え字 の決め方で、 もう一つはサイズ縮小の値ある。 これは二つの背反な要求から来る。 まず、 多項式時間のアルゴリズムを得るために、効率的に問題のサイズ縮小を行う必要がある。 また、比の値 が小さすぎると の推定量が誤差に敏感に$\lceil s_{n}/m\rceil$$k$$\rho$$\rho$なるので$\rho$はある程度十分大きくなければならな

定理 1 添え字 $k\in\{1, \ldots, m\}$ は$\mathrm{A}\mathrm{a}_{\text{。}}$

この目的で我々は次の定理を得た。
を満たすとする。 この時、

$r_{k}=\mathrm{m}\mathrm{f}\mathrm{f}\mathrm{i}\{r_{1}, \ldots, r_{m}\}$$| \Omega_{k}|\geq\frac{1}{m}|\Sigma r,s|$が成り立つ。
定理 1 から、 我々のスキームに現れる理論比$\rho$

は $\rho\geq 1/m$ を満たすことが分かる。 以下、 定理 1を示す手順について説明する。
まず、分割表の個数について以下の 2 つの補題を示すことができる。
補 ffi 1 ベクトノレ $r,$
$N$$r’\in \mathbb{Z}_{+}^{2}$

およひを満たすとする。 この時、集合補題 2 ベクトノレ $r,$$r’\in \mathbb{Z}_{+}^{m}$$\sum_{i=1}^{m}r_{i}’=\sum_{j=1}^{n}s_{j}=N$

およひは $|r_{1}-r_{2}|\leq|r_{1}’-r_{2}’|,$ $\sum_{i=1}^{2}r_{i}=\sum_{1=1}^{2}.r_{\dot{l}}’=\sum_{j=1}^{n}s_{j}=$$s\in \mathbb{Z}_{+}^{n}$$|\Sigma_{r,s}|$$s\in \mathbb{Z}_{+}^{n}$と $|\Sigma r’,s|$ [こついて $|\Sigma r,s|\geq|\Sigma r’,s|$ が成り立つ。
は $|r_{1}-r_{2}|\leq|r_{1}’-r_{2}’|,$ $r_{i}=r_{i}’(i=3, \ldots, m),$を満たすとする。 この時、集合$|\Sigma r,s|$と$|\Sigma r’,s|$について$\sum_{i=1}^{m}r_{\dot{l}}=$$|\Sigma_{r,s}|\geq|\Sigma_{r\prime,s}|$が成り立つ。
いま、 点 $y=(y_{1}, \ldots, y_{m})$ は$\mathbb{R}^{m}$上の点とする。 $N^{-}= \sum_{j=1}^{n-1}s_{j}=N-s_{n}$ として、$\mathbb{R}^{m}$

空間上に超平面$S=\mathrm{d}\mathrm{e}\mathrm{f}$. $\{y\in \mathbb{R}^{m}|\sum_{i=1}^{m}y_{i}=N^{-}\}$を定義する。 また、 $r=(r_{1}, \ldots, r_{m})$ こ対し、 $(m-1)$ 次元単体

$S^{+}$と$S^{+}=\{y\in S\mathrm{d}\mathrm{e}\mathrm{f}.|y_{i}\geq 0, i=1, \ldots,m\}$$S_{0}=\{y\in S\mathrm{d}\mathrm{e}\mathrm{f}.|r_{i}-s_{n}\leq y_{i}\leq r:, i=1, \ldots, m\}$

で定める。
-

いま、

$x\in \mathbb{R}^{m}$を.

$x=r-y\mathrm{d}\mathrm{e}\mathrm{f}$$S_{0}$

をそれぞれ、(1),,と定める。 任意の $x$ [こ対して、 $x_{i}=X_{in}$(2)$(i=1, \ldots, m)$となる分割表 $X\in\Sigma_{r,s}$ が存在するための必要十分条件は、$(i=1, \ldots, m)$ ,$0 \leq X:\leq\min\{r:, s_{n}\}$$\sum_{\dot{\iota}=1}^{m}x_{i}=s_{n}$,である。 この条件は $x=r-y$ より、領域$x\in \mathbb{Z}^{m}$$S^{+},$

任意の $k,$ $l\in\{1, \ldots, m\}$ , k\neq \sim こ対して、$S_{0}$$\mathrm{R}^{m}$,を用いて$\exists y\in S^{+}\cap S0\cap \mathbb{Z}^{m}$

空間上の超平面$H_{kl}$を$H_{kl}=\{y\in \mathbb{R}^{m}|-y_{k}+y_{l}=-r_{k}+r_{l}\}$で定義する。 また、 $i\in\{1, \ldots, m\}$ に対して領域$S_{:}$(3)を. $= \max\{r_{l}-y_{l}|l=1, \ldots, m\}\}$$S_{i}=\{y\in S_{0}\mathrm{d}\mathrm{e}\mathrm{f}.|r_{\dot{\mathrm{t}}}-y_{1}$と表される。

(4)

218と定義する。 さらに、 任意の点 $y\in S$ に対して関数 $q(y)$ を$q(y)=\{$$|\Sigma y,s-|$$(y\in \mathbb{Z}_{+}^{m}\cap S0\cap S^{+})$,(5)0(otherwise),と定義する。 但し、 $s^{-}=(s_{1}, \ldots, s_{n-1})$ である。 明らか [こ、 $| \Sigma r,s|=\sum_{y\in s_{0}}q(y)$ が成り立つ。

$\sim$の時、 $q(y)$ について以下の補題が成り立つ。

補題 3 行の添え字 $k\in\{1, \ldots, m\}$ は $rk= \max\{r_{1}, \ldots, r_{m}\}$ を溝たすとする。また、$l\in\{1, \ldots, m\},$に関して対称な点 $y^{*}=h_{kl}(y)$ は $q(y)\geq q(y^{*})$ を漢たす。
に対し $y\in S_{k}$ の$k$$l\neq$$H_{kl}$

任意の $i\in\{1, \ldots, m\}$ に対し、関数が威り立つ。 いま、 次の補題を得る。

.

$Q_{i}= \sum_{y\in s_{\mathrm{s}}}q(y)\mathrm{d}\mathrm{e}\mathrm{f}$補題 4 行の添え字$k$を定義する。明らかに、 $| \Sigma r,s|\leq\sum_{i=1}^{m}Q_{i}$は $r_{k}= \max\{r_{1}, \ldots, r_{m}\}$ を満たすとする。 この時、任意の $i\in\{1, \ldots, m\}$ に対して、 $Q_{k}\geq Q_{i}$ が成り立つ。

補題 4 から、 行の添え字が成り立つ。 また、$k$が $r_{k}= \max\{r_{1}, \ldots, r_{m}\}$ を満たすとき、 $| \Sigma r,s|\leq\sum_{i=1}^{m}Q_{i}\leq mQ_{k}$の定義より$Q_{k}$$Q_{k}$$=$$|\{X\in\Sigma r,s\}|\forall iX_{kn}\geq X_{in}\}|$$\leq$$| \{X\in\Sigma r,s|X_{kn}\geq\frac{1}{m}s_{n}\}|$$=|\Omega_{k}|$が成り立つ。 すなわち $|\Sigma r,s|\leq m|\Omega_{k}|$ となり、 定理 1 を得る。

4

推定量の誤差と偏り
ここでは、我々の近似解法で得られる近似解の誤差と偏りについて議論する。提案した近似解法中で $m\mathrm{x}n$ 分割表の一様標本抽出を仮定した。 しかし、一般に $m\mathrm{x}n$ 分割表の一様標本抽出は困難であるため、我々の提案するスキームでは近似一様標本抽出を仮定する。集合
いだの総分布距離$\mathrm{d}_{\mathrm{T}\mathrm{V}}(\pi, \nu)$上の–様分布をで表す。 任意の正数標本は$\Sigma_{r,s}$

$\pi$

上の分布

$\nu$

を、

$\mathrm{d}_{\mathrm{T}\mathrm{V}}(\pi, \nu)=\frac{1}{2}\sum_{x\in\Omega}|\pi(x)-\nu(x)|$$\epsilon<1$を満たすとする。 但し、$R$

は$Z$$\nu$

上の分布$\pi$と$\nu$で定義する。 いま、集合に対して、 ある近似的なに従って抽出されるものとし、分布$\Omega$$–arrow$のあ$\Sigma_{r,s}$様標本抽出法が存在して、は総分布距離$\mathrm{d}_{\mathrm{T}\mathrm{V}}(\pi, \nu)\leq\epsilon/(6mR)$を得るために必要な問題のサイズ縮小の回数とする。 この標本抽出法を用いて標本数 $\mathrm{A}f=108mR^{2}\epsilon^{-2}\ln(2R/\delta)$ のモンテカルロ法を行い、 前節で提案した手法における比$\rho$の推定量$\hat{\rho}$を計算する。 逐次この推定比を用いて$|\Sigma r,s|$の推定量 $Z$ を求める。 すなわ

ち我々のスキームは多項式回の問題縮小と多項式個の標本抽出で終丁する。 この時、 推定 5関して次の -.. つの定理が成り #. つ。
定理 2 推定量 $Z$ は$\mathrm{P}\mathrm{r}[(1-\epsilon)|\Sigma r,s|\leq Z\leq(1+\epsilon)|\Sigma r,s|]\geq 1-\delta$を満たす。
定理 2 は $[4, 3]$ と同様な以下の手順で示すことができる。

1.

$1\leq\forall i\leq R,\hat{\rho}_{i}=\mathrm{E}\mathrm{d}\mathrm{e}\mathrm{f}.[U_{i}/M]\geq 1/m-\epsilon/6mR$

,

$Z$に

218 2192.

$| \rho_{i}-\hat{\rho}_{i}|\leq\frac{\xi \mathrm{i}}{6R-\epsilon}\hat{\rho}_{l}$,3.

$\mathrm{P}\mathrm{r}[|Z_{l}-\hat{\rho}_{i}|>\frac{\epsilon}{6R-\epsilon}\mathrm{A}^{\cdot}i]\leq\frac{\delta}{R}$4.

$1-\delta$

以上の確率で、 $|(Z_{1}\cdots Z_{R})^{-1}-(\rho_{1}\cdots\rho_{R})^{-1}|\leq\epsilon(\rho_{1}\cdots\rho_{R})^{-1}$ .

定理 3 推定量 $Z$ は$\frac{|\mathrm{E}[Z]-\Sigma r,s|}{|\Sigma r,s|}\leq\frac{\epsilon}{4}+\mathrm{e}^{-90R^{3}\epsilon^{-2}\ln(2R/\delta)}<(\frac{1}{4}+\frac{1}{10^{27}})\epsilon$

を漢たす。
証明:

$Z_{1},$$\ldots,$

は独立なので$Z_{R}$$\mathrm{E}[Z]=\sigma \mathrm{E}[\prod_{i=1}^{R}\frac{1}{Z_{i}}]=\sigma\prod_{\dot{\iota}=1}^{R}\mathrm{E}[\frac{1}{Z_{i}}]$’が成り立つ。 いま、$\mathrm{E}[\frac{1}{Z_{i}}]$$=$$\sum_{U.=0}^{M}\frac{M+1}{U_{i}+1}(\begin{array}{l}MU_{i}\end{array})\hat{\rho}_{i}^{U}(:1-\hat{\rho}_{i})^{M-U_{\mathrm{i}}}$$=$$\frac{1}{\hat{\rho}_{i}}\sum_{U\dot{.}=0}^{M}(\begin{array}{ll}M +1U_{\dot{\mathrm{t}}}+1 \end{array}) \hat{\rho}_{i}^{U_{\dot{\mathrm{t}}}+1}(1-\hat{\rho}_{i})^{M-U}\dot{\cdot}=\frac{1}{\hat{\rho}_{i}}\{1-(1-\hat{\rho}_{i})^{M+1}\}$,である。 したがって等式$\mathrm{E}[Z|$

$=$

(6)

$\mathrm{e}\tau\prod_{i}\frac{1}{\hat{\rho}_{i}}\{1-(1-\hat{\rho}_{i})^{M+1}\}$

が成り立つ。 式 (6) より、 $||\Sigma_{r,s}|-\mathrm{E}[Z]|$ は.

$|| \Sigma_{r,s}|-\mathrm{E}[Z]|---|\sigma\prod$$\frac{1}{\rho_{i}}-\sigma\prod_{i}\frac{1}{\hat{\rho}_{i}}\{1-(1-\hat{\rho}_{i})^{M+1}\}|$$\leq_{-}|\sigma\prod_{i}\frac{1}{\rho_{i}}-\sigma\prod_{i}\frac{1}{\hat{\rho}_{i}}||\{1--(1-\hat{\rho}_{i})^{M+1}\}|+|\sigma\prod_{i}\frac{1}{\rho_{i}}(1-\hat{\rho}_{i})^{M+1}|$’を満たす。 従って、$| \sigma\prod_{i}\frac{1}{\rho_{\mathrm{i}}}-\sigma\prod_{i}\frac{1}{\hat{\rho}_{i}}|=\sigma\prod_{i}\frac{1}{\rho_{i}}|1-\prod_{i}\frac{\rho j}{\hat{\rho}_{i}}|\leq|\Sigma r_{\tau}s|\{(1+\frac{\epsilon}{6R-\epsilon})R- 1\}$,$\leq|\Sigma r_{\mathrm{I}}s|\{\exp(\frac{\epsilon R}{6R-\epsilon})- 1\}\leq|\Sigma r.s|\frac{\epsilon R}{6R-\epsilon-\epsilon R}\leq\frac{\epsilon}{4}|\Sigma_{T.\mathit{8}}|$が成り立つ。 よって$|| \Sigma_{r,s}|-\mathrm{E}[Z]|\leq\frac{\epsilon}{4}|\Sigma r,s|$を得る。 最後に$\prod_{i}(1-\hat{\rho}_{i})^{M+1}$$\leq$$\prod_{i}(1-\hat{\rho}_{i})^{M+1}$

十 $| \Sigma r,s|\prod_{i}(1-\hat{\rho}_{i})^{M+1}$ ,の大きさを$\{1-(\frac{1}{m}-\frac{\epsilon}{6mR})\}^{RM}=\{$ $1-( \frac{1-\frac{\epsilon}{6R}}{m})\}^{RM}$220$\leq$$\leq$$=$$(1- \frac{5}{6m})^{RM}\leq(1-\frac{5}{6m})^{R108mR^{2}c^{-2}\ln(2R/\delta)}$$(\mathrm{e}^{-1})$

命$108mR^{3}\epsilon^{-2}\ln(2R/\delta)$$\mathrm{e}^{-90R^{3}\epsilon^{-2}\ln(2R/\delta)}$.

従って、 偏りの大きさは以下のようになる。
$||\Sigma r,s|-\mathrm{E}[Z]|$$\leq$$\frac{\epsilon}{4}|\Sigma r,s|+\mathrm{e}^{-90R^{\theta}\epsilon^{-2}\ln(2R/\delta)}|\Sigma r,s|$$\leq$$| \Sigma r,s|(\frac{\epsilon}{4}+\mathrm{e}^{-90R^{3}\epsilon^{-2}\ln(2R/\delta)})$.

口

この結果から、提案したアルゴリズムによって得られる近似解と厳密解との偏りの大きさは、主に標本分布と一様分布とのずれに依存し、モンテカルロ法を採用したことによる偏りはほとんど無いと言える。

5

結論と課題

分割表数え上げ問題に対して多項式回の問題縮小と多項式個の標本抽出で誤差の大きさと偏りの幅を確率論的に押さえられた近似解を与える。 2002 年 Crym et al.
が行数が定数の 2 元分割表に対する heat bath マルコフ連鎖力. rapid mixing であることを示した
分割表数[1]。 もちろんこれを我々のスキームに適用すると、我々のスキームも行数固定の分割表に対するえ上げ問題に対する多項式時間確率論的近似解法となる。 しかし、 –般の多項式時間近似一様標本抽出法の存在性は現時点では未解決である。
我々のスキームは、$m\mathrm{x}n$$\mathrm{a}$$m\mathrm{x}n$$m\mathrm{x}n$参考文献[1] M. Cryan, M. Dyer, L. A. Goldberg, M. Jerrum, and R. Martin, ‘Rapidly mixing Markovchains for samPling contingency tables with constant number of rows,” Prvceedings of theAnnual Symposium on Foundations of Computer Science (FOCS) (2002), PP. 711-720.
$\mathrm{t}$$\mathit{4}\mathit{3}rd$[2] M. Dyer, A. Frieze, and R. Kannan, “A random polynomial time algorithm for approximatthe volume of convex bodies,” Journal of the $ACM,$ $38$ (1991), PP. 1-17.
$\mathrm{i}\mathrm{n}\mathrm{g}$[3] M. Dyer and C. Greenhill, “Polynomial-time counting and sampling of two rowed contin-gency tables,” Theoretical Computer Sciences, 246 (2000), PP. $265\cdot-278$ .

[4] M. Dyer, R. Kannan, and J. Mount, “Sampling contingency tables,” Random Structures
and Algorithms, 10 (1997), PP. 487-506.
[5] D. Hernek, “Random generation of 213(1998), pp. 71-79.

$\mathrm{x}n$contingency tables,” Random Struct. Algorithms,[6] W. H\"offding, “Probability inequalities for sums of bounded random variables,” J. Amer.
Stat. Assoc. 58(1963), PP. 13-30.

[7] S. Kijima and T. Matsui, “Approximate counting scheme forTech. Rep., The Univ. of Tokyo, METR 03-01, January (2003).

$m\mathrm{x}n$contingency tables,”