ギブスサンプリング

統計学と統計物理学において、ギブスサンプリング（英: Gibbs sampling, Gibbs sampler）は、直接サンプリングが難しい確率分布の代わりにそれを近似するサンプル列を生成するMCMC法（Markov chain Monte Carlo algorithm）の1つである。この生成された数列は、同時分布や周辺分布や期待値などの積分計算を近似するために用いられる。通常は観測として与えられている変数に関してはサンプリングをする必要はない。ギブスサンプリングは統計的推定やベイズ推定の手法として頻繁に用いられている。ランダムアルゴリズムであり、変分ベイズ法（英語版）（variational Bayes）やEMアルゴリズム（expectation-maximization algorithm）のような統計的推定法のための決定論的な方法の代替法である。

他のMCMC法と同様に、ギブスサンプリングはサンプルのマルコフ連鎖を生成する。得られるサンプル列がマルコフ連鎖であるため、例えば100番目毎にサンプルを選ぶといったサンプルが十分に独立とみなせるように気をつけるべきである。それに加え、サンプル列の始めの方の値は目的の分布を精確には表していないため、初期値を与えたすぐ後はburn-in期間としてサンプルを捨てるべきである。

導出

ギブスサンプリングはメトロポリス・ヘイスティングス法の1つである。同時分布より周辺化された条件付き確率分布から、与えられた確率分布に従ったサンプルをサンプリングする。同時確率${\displaystyle (x_{1},\dotsc ,x_{n})}$から確率変数${\displaystyle {\boldsymbol {X}}=\{x_{1},\dotsc ,x_{n}\}}$ を${\displaystyle \left.k\right.}$サンプル得たい。${\displaystyle n}$次元の${\displaystyle i}$番目のサンプルを${\displaystyle {\boldsymbol {X}}^{(i)}=\{x_{1}^{(i)},\dotsc ,x_{n}^{(i)}\}}$とする。

手順は以下の通りである。

1. 初期値${\displaystyle {\boldsymbol {X}}^{(0)}}$を設定する。
2. 条件付き確率${\displaystyle p(x_{j}|x_{1}^{(i)},\dotsc ,x_{j-1}^{(i)},x_{j+1}^{(i-1)},\dotsc ,x_{n}^{(i-1)})}$から${\displaystyle i}$番目のサンプル${\displaystyle x_{j}^{(i)}}$をサンプリングする。
3. ${\displaystyle i}$が${\displaystyle k}$以下であれば${\displaystyle i=i+1}$して2へ。それ以外だと終了。

サンプルは他の変数に条件付けされた分布からサンプリングされるが、他の変数には新しいサンプルを使用すべきである。つまりは、1ステップ毎に1変数をサンプルし、新しいサンプルに入れ替えていく。このサンプル集合は、全ての変数に関する同時分布を近似する。

それに加え、全てのサンプルに関して平均をとることで期待値を近似することができる。

以下は注意点である。

• 初期値の設定にEMアルゴリズムなどを用いたり、ランダムに決めてもいい。しかし初期値の設定に敏感にならなくても、burn-inの期間を設ければ問題ではない。
• サンプリング初期の値を捨てる期間（burn-in period）を設けることがよく行われる。また、期待値を計算するときには${\displaystyle n}$番目毎の値しか用いないといったこともよく行われる。この理由には2つある。1つは生成されるサンプル列はマルコフ連鎖でありサンプル間にある程度の相関が存在するため独立なサンプルではない。もう1つはマルコフ連鎖の定常分布は目的の分布になるが初期値から定常分布に到達するまでには時間がかかる。 自己相関の量や${\displaystyle n}$をアルゴリズムで決定することもできるが、それは黒魔術的である。
• 焼きなまし法(Simulated Annealing)はサンプリングの初期によく用いられる。しかしサンプル空間の移動が遅くサンプルの相関が強くなってしまう。その他の自己相関を減少させる方法には collapsed Gibbs samplingやblocked Gibbs sampling、ordered overrelaxationなどがある。

条件付き分布と同時分布の関係

他のすべての変数が与えられたときのある変数に関する条件付き確率は同時確率に比例する。

${\displaystyle p(x_{j}|x_{1},\dotsc ,x_{j-1},x_{j+1},\dotsc ,x_{n})={\frac {p(x_{1},\dotsc ,x_{n})}{p(x_{1},\dotsc ,x_{j-1},x_{j+1},\dotsc ,x_{n})}}\propto p(x_{1},\dotsc ,x_{n})}$

ここで比例するとは、分母が${\displaystyle x_{j}}$の関数ではなく、${\displaystyle x_{j}}$がどんな値であれ分母が定数であることである。周辺化定数は同時確率を${\displaystyle x_{j}}$に関して周辺化した値である。

${\displaystyle p(x_{1},\dotsc ,x_{j-1},x_{j+1},\dotsc ,x_{n})=\int p(x_{1},\dotsc ,x_{n})dx_{j}}$

実用上、確率変数${\displaystyle x_{j}}$の条件付き分布を求めるためのもっとも簡単な方法はグラフィカルモデルの変数のうち${\displaystyle x_{j}}$の関数ではない値を独立とみなして因子化すればいい。

そのほかには、3つの場合が考えられる。

1. 分布が離散分布である場合、${\displaystyle x_{j}}$の取りうる値すべてに関して確率を計算し、足しあわせることで周辺化定数を計算すればいい。
2. 分布が連続で既知の形を持つ場合、1次元の周辺化なので周辺化定数が計算可能である。
3. その他の場合、周辺化定数は無視することができる。大抵のサンプリング法は周辺化定数がなくても問題ではない。

理論

サンプル${\displaystyle \left.X\right.}$を次元${\displaystyle \left.d\right.}$のパラメータベクトル${\displaystyle \theta \in \Theta \,\!}$に基づいた事前分布${\displaystyle g(\theta _{1},\dotsc ,\theta _{d})}$からサンプリングする。

${\displaystyle \left.d\right.}$が大きい場合、数値積分を行ない${\displaystyle \left.\theta _{i}\right.}$の周辺化分布を計算することは困難を伴う。

その周辺化分布の計算を${\displaystyle \left.\Theta \right.}$空間上のマルコフ連鎖を生成することで代替する。以下の2ステップを繰り返す。

1. ${\displaystyle j}$を選ぶ${\displaystyle 1\leq j\leq d}$
2. ${\displaystyle g(\theta _{1},\dotsc ,\theta _{j-1},\,\cdot \,,\theta _{j+1},\dotsc ,\theta _{d})}$の分布にしたがい、新しい変数${\displaystyle \left.\theta _{j}\right.}$を選ぶ。

この2ステップは分布${\displaystyle \left.g\right.}$にしたがう可逆なマルコフ連鎖を生成する。

これは以下の通り証明される。すべての${\displaystyle i\neq j}$ に関して${\displaystyle \left.x_{i}=y_{i}\right.}$であれば、${\displaystyle x\sim _{j}y}$と表す。${\displaystyle x\sim _{j}y}$ は同値関係である。${\displaystyle \left.p_{xy}\right.}$は${\displaystyle x\in \Theta }$から${\displaystyle y\in \Theta }$へ遷移する確率の分布を表す。

よって遷移確率は

${\displaystyle p_{xy}={\begin{cases}{\dfrac {1}{d}}{\dfrac {g(y)}{\sum _{z\in \Theta :z\sim _{j}x}g(z)}}&x\sim _{j}y\\0&{\text{otherwise}}\end{cases}}}$

したがって、

${\displaystyle g(x)p_{xy}={\frac {1}{d}}{\frac {g(x)g(y)}{\sum _{z\in \Theta :z\sim _{j}x}g(z)}}={\frac {1}{d}}{\frac {g(y)g(x)}{\sum _{z\in \Theta :z\sim _{j}y}g(z)}}=g(y)p_{yx}}$

このように詳細つりあい条件は満たされる。

実用的には${\displaystyle \left.j\right.}$はランダムに選ばれず、並びの順番で選ばれる。 また正確に言うと、これは非定常マルコフ連鎖を生成するが、各ステップは可逆であり、全体のプロセスは求めたい定常分布を与える。

手法の拡張

ギブスサンプリングの拡張について説明する。これらの拡張はサンプル間の自己相関を減少させることで、計算コストを減らすことを目標に考えられている。

Blocked Gibbs sampler

blocked Gibbs sampler は個々の変数を考えるのではなく、変数を複数のグループに分割して条件づき分布を考える。 例えば、隠れマルコフモデルではforward-backward algorithmを使い、隠れ変数に関するマルコフ連鎖を生成する。

Collapsed Gibbs sampler

collapsed Gibbs samplerは周辺化分布の変数を積分消去する。たとえば、A、B、Cの3つの変数があるとする。ギブスサンプリングではp(A|B,C)、p(B|A,C)、p(C|A,B)からサンプリングすることになる。collapsed Gibbs samplerではAをサンプリングする時には例えばBを積分消去した周辺化分布p(A|C)からサンプリングを行う。Aに関してBが共役事前分布であったり、指数分布族であれば計算が容易にできる。詳しくはcompound distribution、Liu (1994)を参考。^[1]

ソフトウェア

• OpenBUGS (Bayesian inference Using Gibbs Sampling) ：複雑な統計モデルのベイズ推定にMCMC法を用いている
• JAGS (Just another Gibbs sampler) MCMC法を用いた階層ベイズモデルを推定できる
• Church ：任意の分布をギブスサンプリングする

参照