ギブスサンプリング

統計学において、ギブスサンプリング（またはギブスサンプラー）は、統計力学ではヒートバスアルゴリズムとも呼ばれ、同時分布からの直接サンプリングが困難な場合に、条件付き分布からのサンプリングの方が実用的であるような、指定された多変量確率分布からのサンプリングを行うためのマルコフ連鎖モンテカルロ（MCMC）アルゴリズムです。このシーケンスは、同時分布を近似する（例えば、分布のヒストグラムを生成する）、変数の1つまたは変数のサブセットの周辺分布を近似する（例えば、未知のパラメータや潜在変数）、または積分を計算する（例えば、変数の1つの期待値を計算する）ために使用できます。通常、変数の中には値が既知の観測値に対応するものがあり、それらはサンプリングする必要はありません。

ギブスサンプリングは、統計的推論、特にベイズ推論の手法として広く用いられています。これはランダム化アルゴリズム（つまり、乱数を用いるアルゴリズム）であり、期待値最大化アルゴリズム（EMアルゴリズム）などの決定論的アルゴリズムに代わる統計的推論手法です。

他のMCMCアルゴリズムと同様に、ギブスサンプリングはマルコフ連鎖を生成するサンプル群であり、各サンプルは近傍のサンプルと相関関係にある。そのため、独立したサンプルが必要な場合は注意が必要である。連鎖の最初（バーンイン期間）のサンプルは、目的とする分布を正確に表していない可能性があり、通常は破棄される。

ギブスサンプリングは、サンプリングアルゴリズムと統計物理学との類似性にちなんで、物理学者ジョサイア・ウィラード・ギブスにちなんで名付けられました。このアルゴリズムは、ギブスの死後約80年後の1984年にスチュアートとドナルド・ゲーマン兄弟によって記述され^{[ 1 ]}、周辺確率分布、特に事後分布の計算のために統計学コミュニティで普及しました^{[ 2 ] 。}

ギブスサンプリングは、基本バージョンではメトロポリス・ヘイスティングスアルゴリズムの特殊なケースです。しかし、拡張バージョン（下記参照）では、各変数（場合によっては各変数グループ）を順番にサンプリングすることで、多数の変数からサンプリングを行うための一般的なフレームワークとみなすことができ、メトロポリス・ヘイスティングスアルゴリズム（またはスライスサンプリングなどの手法）を組み込んで、1つ以上のサンプリングステップを実行できます。

ギブスサンプリングは、同時分布が明示的にわかっていない場合や、直接サンプリングするのが難しい場合に適用できますが、各変数の条件付き分布はわかっており、サンプリングが容易（または少なくとも容易）な場合にも適用できます。ギブスサンプリングアルゴリズムは、他の変数の現在の値を条件として、各変数の分布からインスタンスを順番に生成します。サンプルのシーケンスはマルコフ連鎖を構成し、そのマルコフ連鎖の定常分布がまさに求められる同時分布であることが示されます。^{[ 3 ]}

ギブスサンプリングは、ベイジアンネットワークの事後分布をサンプリングするのに特に適しています。なぜなら、ベイジアンネットワークは通常、条件付き分布の集合として定義されるからです。

ギブスサンプリングは、その基本的な形態では、メトロポリス・ヘイスティングスアルゴリズムの特殊なケースです。ギブスサンプリングのポイントは、多変量分布が与えられた場合、同時分布で積分して周辺化するよりも、条件付き分布からサンプリングする方が簡単であるということです。 を得たいとしましょう。${\displaystyle k}$サンプル${\displaystyle n}$次元ランダムベクトル${\displaystyle \mathbf {X} =(X_{1},\dots ,X_{n})}$我々は反復的に進める。

• 初期値から始めましょう$\displaystyle \mathbf {X} ^{(0)}}$。
• サンプルが与えられた場合${\displaystyle \mathbf {X} ^{(i)}=\left(x_{1}^{(i)},\dots ,x_{n}^{(i)}\right)}$次のサンプルを取得するために${\displaystyle \mathbf {X} ^{(i+1)}=\left(x_{1}^{(i+1)},x_{2}^{(i+1)},\dots ,x_{n}^{(i+1)}\right)}$各コンポーネントをサンプリングできます${\displaystyle x_{j}^{(i+1)}}$分布から${\displaystyle X_{j}}$これまでにサンプリングされたすべてのコンポーネントを条件とします。${\displaystyle X_{\ell }^{(i+1)}}$すべての人々のために${\displaystyle \ell \leq j-1}$、そして${\displaystyle X_{\ell }^{(i)}}$のために${\displaystyle j+1\leq \ell \leq n}$言い換えれば、我々はサンプリングする${\displaystyle x_{j}^{(i+1)}}$分布によると${\displaystyle P\left(X_{j}=\cdot |X_{1}=x_{1}^{(i+1)},\dots ,X_{j-1}=x_{j-1}^{(i+1)},X_{j+1}=x_{j+1}^{(i)},\dots ,X_{n}=x_{n}^{(i)}\right)}$。

このようなサンプリングを実施する場合、以下の重要な事実が当てはまります。

• これらのサンプルは、すべての変数の同時分布を近似的に表している。
• 変数の任意の部分集合の周辺分布は、その部分集合の変数のサンプルのみを考慮し、残りの変数を無視することで近似できる。
• 任意の変数の期待値は、すべてのサンプルについて平均をとることによって近似できる。

サンプリングを実施する際：

• 変数の初期値は、ランダムに決定することも、期待値最大化などの他のアルゴリズムによって決定することもできます。
• 最初にサンプリングする変数の初期値を決定する必要は実際にはありません。
• 最初はいくつかのサンプルを無視し（いわゆるバーンイン期間）、その後はすべてのサンプルのみを考慮するのが一般的です。${\displaystyle n}$期待値を計算する際の平均値の計算には、 番目のサンプルを使用します。たとえば、最初の 1,000 個のサンプルは無視され、その後 100 番目ごとのサンプルが平均化され、残りはすべて破棄されます。その理由は、(1)マルコフ連鎖の定常分布は変数に関する望ましい同時分布ですが、その定常分布に到達するには時間がかかる場合があること、(2) 連続するサンプルは互いに独立ではなく、ある程度の相関を持つマルコフ連鎖を形成することです。場合によっては、アルゴリズムを使用してサンプル間の自己相関の量と値を決定することができます。${\displaystyle n}$（実際に使用されるサンプル間の期間）はここから計算されますが、実際にはかなりの「黒魔術」が関わっています。
• シミュレーテッドアニーリング法は、サンプリングプロセスの初期段階における「ランダムウォーク」挙動（すなわち、望ましいように素早く移動するのではなく、サンプル空間内をゆっくりと移動し、サンプル間の自己相関が高くなる傾向）を低減するためによく用いられます。自己相関を低減する他の手法としては、コラプテッドギブスサンプリング、ブロックドギブスサンプリング、オーダードオーバーリラクゼーションなどがあります（下記参照）。

ギブスサンプリングは、統計的推論（例えば、特定の日に特定の店舗で買い物をする人の数、有権者が最も投票する可能性の高い候補者など、パラメータの最適値を決定する）によく用いられます。この手法の考え方は、観測データごとに個別の変数を作成し、それらの変数からサンプリングするのではなく、対象となる変数を観測値に固定することで、観測データをサンプリングプロセスに組み込むというものです。残りの変数の分布は、観測データに基づいて条件付けられた事後分布となります。

目的のパラメータの最も可能性の高い値（最頻値）は、最も頻繁に出現するサンプル値を選択することで簡単に選択できます。これは、本質的にパラメータの最大事後推定に相当します。（パラメータは通常連続であるため、最頻値の意味のある推定値を得るには、サンプリングされた値を有限個の範囲または「ビン」のいずれかに「ビン化」する必要がある場合がよくあります。）しかし、より一般的には、サンプリングされた値の期待値（平均または平均値）が選択されます。これは、ベイズサンプリングから得られる分布全体に関する追加データを利用するベイズ推定器であり、期待値最大化（EM）などの最大化アルゴリズムは、分布から単一の点しか返すことができません。たとえば、単峰性分布の場合、平均（期待値）は通常、最頻値（最も頻繁に出現する値）と似ていますが、分布が一方の方向に歪んでいる場合、平均はその方向に移動します。これにより、その方向の余分な確率質量が効果的に考慮されます。 （分布が多峰性の場合、期待値は意味のある値を返すとは限らず、通常はいずれかの最頻値を用いる方が適切です。）

変数の中には、通常、関心のあるパラメータに対応するものもありますが、変数間の関係を適切に表現するためにモデルに導入される、関心のない（「邪魔な」）変数もあります。サンプリングされた値はすべての変数の同時分布を表しますが、期待値や最頻値を計算する際には、邪魔な変数は単純に無視できます。これは、邪魔な変数を周辺化することと同等です。複数の変数の値が必要な場合は、各変数について個別に期待値を計算します。（ただし、最頻値を計算する場合は、すべての変数をまとめて考慮する必要があります。）

教師あり学習、教師なし学習、半教師あり学習（欠損値を含む学習とも呼ばれる）はすべて、値が既知のすべての変数の値を固定し、残りの変数からサンプリングすることで処理できます。

観測データの場合、例えば一連の観測値の標本平均や標本分散に対応する変数ではなく、各観測値に対して1つの変数が存在します。実際には、一般的に「標本平均」や「標本分散」といった概念に対応する変数は全く存在しません。その代わりに、未知の真の平均と真の分散を表す変数が存在し、これらの変数の標本値はギブスサンプラーの動作によって自動的に決定されます。

一般化線形モデル（すなわち線形回帰の変種）も、ギブスサンプリングで処理できる場合があります。例えば、回帰係数に正規分布の事前分布を設定して、特定の二値（はい／いいえ）選択の確率を決定するプロビット回帰は、追加の変数を追加して共役性を利用できるため、ギブスサンプリングで実装できます。しかし、ロジスティック回帰はこの方法では処理できません。一つの方法として、ロジスティック関数を（通常7～9種類の）正規分布の混合で近似する方法があります。ただし、より一般的には、ギブスサンプリングの代わりにメトロポリス・ヘイスティングス法が使用されます。

サンプルを仮定します${\displaystyle \left.X\right.}$パラメータベクトルに依存する分布から取得される${\displaystyle \theta \in \Theta \,\!}$長さ${\displaystyle \left.d\right.}$事前分布あり${\displaystyle g(\theta _{1},\ldots ,\theta _{d})}$もしかしたら${\displaystyle \left.d\right.}$は非常に大きく、その周辺密度を求めるための数値積分は${\displaystyle \left.\theta _{i}\right.}$計算コストが高くなるだろう。そこで、周辺密度を計算する別の方法として、空間上にマルコフ連鎖を作成する方法がある。${\displaystyle \left.\Theta \right.}$以下の2つの手順を繰り返すことで：

1. ランダムなインデックスを選択${\displaystyle 1\leq j\leq d}$
2. 新しい値を選択してください${\displaystyle \left.\theta _{j}\right.}$によると${\displaystyle g(\theta _{1},\ldots ,\theta _{j-1},\,\cdot \,,\theta _{j+1},\ldots ,\theta _{d})}$

これらの手順は、目的の不変分布を持つ可逆マルコフ連鎖を定義する。${\displaystyle \left.g\right.}$これは次のように証明できます。${\displaystyle x\sim _{j}y}$もし${\displaystyle \left.x_{i}=y_{i}\right.}$すべての人々のために${\displaystyle i\neq j}$そして${\displaystyle \left.p_{xy}\right.}$からジャンプする確率を表す${\displaystyle x\in \Theta }$に${\displaystyle y\in \Theta }$すると、遷移確率は次のようになる。

${\displaystyle p_{xy}={\begin{cases}{\frac {1}{d}}{\frac {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 x\sim _{j}y}$これは同値関係である。したがって、詳細平衡方程式が満たされ、連鎖は可逆であり、不変分布を持つことを意味する。${\displaystyle \left.g\right.}$。

実際には、この指標は${\displaystyle \left.j\right.}$はランダムに選択されるのではなく、チェーンはインデックスを順番に巡回します。一般に、これは非定常マルコフ過程となりますが、個々のステップは可逆であり、全体的なプロセスは依然として望ましい定常分布を持ちます（チェーンが固定された順序の下で全ての状態にアクセスできる限り）。

させて${\displaystyle y}$は、標本分布から生成された観測値を表す。${\displaystyle f(y|\theta )}$そして${\displaystyle \pi (\theta )}$パラメータ空間でサポートされている事前分布である${\displaystyle \Theta }$ベイズ統計学の中心的な目標の一つは、事後確率密度を近似することである。

${\displaystyle \pi (\theta |y)={\frac {f(y|\theta )\cdot \pi (\theta )}{m(y)}}}$

周辺尤度${\displaystyle m(y)=\int _{\Theta }f(y|\theta )\cdot \pi (\theta )d\theta }$すべての${\displaystyle y}$。

ギブスサンプラーを説明するために、パラメータ空間が${\displaystyle \Theta }$分解すると

${\displaystyle \Theta =\prod _{i=1}^{K}\Theta _{i}=\Theta _{1}\times \cdots \Theta _{i}\times \cdots \times \Theta _{K},\quad \quad (K>1)}$、

どこ${\displaystyle \times }$はデカルト積を表します。各コンポーネントのパラメータ空間${\displaystyle \Theta _{i}}$スカラー成分、サブベクトル、または行列の集合である。

集合を定義する${\displaystyle \Theta _{-i}}$それは${\displaystyle \Theta _{i}}$ギブスサンプラーの必須成分は${\displaystyle i}$各 の 番目の完全条件付き事後分布${\displaystyle i=1,\cdots ,K}$

${\displaystyle \pi (\theta _{i}|\theta _{-i},y)=\pi (\theta _{i}|\theta _{1},\cdots ,\theta _{i-1},\theta _{i+1},\cdots ,\theta _{K},y)}$。

以下のアルゴリズムは、一般的なギブスサンプラーの詳細を示しています。

${\displaystyle {\text{Initialize: pick arbitrary starting value}}\,\,\theta ^{(1)}=(\theta _{1}^{(1)},\theta _{2}^{(1)},\cdots ,\theta _{i}^{(1)},\theta _{i+1}^{(1)},\cdots ,\theta _{K}^{(1)})}$

${\displaystyle {\text{Iterate a Cycle:}}\,}$

${\displaystyle \quad \quad {\text{Step 1. draw}}\,\,\theta _{1}^{(s+1)}\sim \pi (\theta _{1}|\theta _{2}^{(s)},\theta _{3}^{(s)},\cdots ,\theta _{K}^{(s)},y)}$

${\displaystyle \quad \quad {\text{Step 2. draw}}\,\,\theta _{2}^{(s+1)}\sim \pi (\theta _{2}|\theta _{1}^{(s+1)},\theta _{3}^{(s)},\cdots ,\theta _{K}^{(s)},y)}$

${\displaystyle \quad \quad \quad \vdots }$

${\displaystyle \quad \quad {\text{Step i. draw}}\,\,\theta _{i}^{(s+1)}\sim \pi (\theta _{i}|\theta _{1}^{(s+1)},\theta _{2}^{(s+1)},\cdots ,\theta _{i-1}^{(s+1)},\theta _{i+1}^{(s)},\cdots ,\theta _{K}^{(s)},y)}$

${\displaystyle \quad \quad {\text{Step i+1. draw}}\,\,\theta _{i+1}^{(s+1)}\sim \pi (\theta _{i+1}|\theta _{1}^{(s+1)},\theta _{2}^{(s+1)},\cdots ,\theta _{i}^{(s+1)},\theta _{i+2}^{(s)},\cdots ,\theta _{K}^{(s)},y)}$

${\displaystyle \quad \quad \quad \vdots }$

${\displaystyle \quad \quad {\text{Step K. draw}}\,\,\theta _{K}^{(s+1)}\sim \pi (\theta _{K}|\theta _{1}^{(s+1)},\theta _{2}^{(s+1)},\cdots ,\theta _{K-1}^{(s+1)},y)}$

${\displaystyle {\text{end Iterate}}}$

ギブスサンプラーは、反復モンテカルロ法によってサイクル内で実行されることに注意してください。${\displaystyle S}$サンプル数${\displaystyle \{\theta ^{(s)}\}_{s=1}^{S}}$上記のアルゴリズムによって描画されたマルコフ連鎖は、不変分布を目標密度として定式化します。${\displaystyle \pi (\theta |y)}$。

さて、それぞれについて${\displaystyle i=1,\cdots ,K}$情報理論上の以下の量を定義する。

${\displaystyle I(\theta _{i};\theta _{-i})={\text{KL}}(\pi (\theta |y)||\pi (\theta _{i}|y)\cdot \pi (\theta _{-i}|y))=\int _{\Theta }\pi (\theta |y)\log {\bigg (}{\frac {\pi (\theta |y)}{\pi (\theta _{i}|y)\cdot \pi (\theta _{-i}|y)}}{\bigg )}d\theta ,}$

${\displaystyle H(\theta _{-i})=-\int _{\Theta _{-i}}\pi (\theta _{-i}|y)\log \pi (\theta _{-i}|y)d\theta _{-i},}$

${\displaystyle H(\theta _{-i}|\theta _{i})=-\int _{\Theta }\pi (\theta |y)\log \pi (\theta _{-i}|\theta _{i},y)d\theta ,}$

すなわち、事後相互情報量、事後微分エントロピー、事後条件付き微分エントロピーである。同様に情報理論的な量も定義できる。${\displaystyle I(\theta _{-i};\theta _{i})}$、${\displaystyle H(\theta _{i})}$、 そして${\displaystyle H(\theta _{i}|\theta _{-i})}$交換することで${\displaystyle i}$そして${\displaystyle -i}$定義された量において。次に、${\displaystyle K}$方程式は成り立つ。^{[ 4 ]}

${\displaystyle I(\theta _{i};\theta _{-i})=H(\theta _{-i})-H(\theta _{-i}|\theta _{i})=H(\theta _{i})-H(\theta _{i}|\theta _{-i})=I(\theta _{-i};\theta _{i}),\quad (i=1,\cdots ,K)}$。

相互情報${\displaystyle I(\theta _{i};\theta _{-i})}$ランダム量の不確実性の減少を定量化する${\displaystyle \theta _{i}}$一度わかったら${\displaystyle \theta _{-i}}$事後的に。それは、以下の場合に限り消滅する。${\displaystyle \theta _{i}}$そして${\displaystyle \theta _{-i}}$は周辺的に独立しており、事後分布となる。相互情報量${\displaystyle I(\theta _{i};\theta _{-i})}$は、から伝達される量として解釈できます。${\displaystyle i}$への第 1 ステップ${\displaystyle i+1}$ギブスサンプラーの単一サイクル内の 番目のステップ。

基本的なギブスサンプラーには数多くのバリエーションが存在する。これらのバリエーションの目的は、計算コストの増加を相殺できる程度にサンプル間の自己相関を低減することである。

ギブスサンプリングが失敗する原因は2つあります。1つ目は、確率の高い状態が孤立して存在し、それらの間に経路がない場合です。例えば、2ビットベクトルの確率分布を考えてみましょう。ベクトル(0,0)と(1,1)はそれぞれ確率が1/2ですが、他の2つのベクトル(0,1)と(1,0)は確率が0です。ギブスサンプリングは、確率の高い2つのベクトルのうちの1つに閉じ込められ、もう一方のベクトルには到達できません。より一般的には、高次元の実数値ベクトルの任意の分布において、ベクトルの2つの特定の要素が完全に相関している（または完全に反相関している）場合、それらの2つの要素は固定されてしまい、ギブスサンプリングでは変更できなくなります。

2 つ目の問題は、すべての状態がゼロでない確率を持ち、高確率状態の島が 1 つしかない場合でも発生する可能性があります。たとえば、100 ビットのベクトルの確率分布を考えてみましょう。すべてゼロのベクトルは確率 1/2 で発生し、他のすべてのベクトルは等しく発生する確率を持ち、したがって確率は1/2です。${\displaystyle {\frac {1}{2(2^{100}-1)}}}$それぞれ。ゼロベクトルの確率を推定したい場合は、真の分布から 100 または 1000 個のサンプルを取るだけで十分です。そうすれば、おそらく 1/2 に非常に近い答えが得られるでしょう。しかし、おそらくそれ以上のサンプルを取る必要があるでしょう。${\displaystyle 2^{100}}$ギブスサンプリングによるサンプルを用いて同じ結果を得る。どんなコンピュータでも、一生かかってもこれは不可能だ。

この問題は、バーンイン期間の長さに関係なく発生します。これは、真の分布ではゼロベクトルが半分の確率で発生し、それらの発生が非ゼロベクトルとランダムに混ざり合っているためです。小さなサンプルでも、ゼロベクトルと非ゼロベクトルの両方が見られます。しかし、ギブスサンプリングでは、長期間（約${\displaystyle 2^{99}}$連続して)、その後は長期間 (約${\displaystyle 2^{99}}$連続して）。したがって、真の分布への収束は非常に遅く、${\displaystyle 2^{99}}$これほど多くのステップを実行することは、妥当な時間内に計算上実現不可能です。ここで収束が遅いのは、次元の呪いの結果であると考えられます。このような問題は、100ビットベクトル全体を一度にブロックサンプリングすることで解決できます。（これは、100ビットベクトルがより大きな変数セットの一部であることを前提としています。このベクトルだけがサンプリングされる場合、ブロックサンプリングはギブスサンプリングを全く行わないことと同等であり、これは仮定上困難です。）

• OpenBUGSソフトウェア（ギブスサンプリングを用いたベイズ推論）は、マルコフ連鎖モンテカルロ法を用いて複雑な統計モデルのベイズ分析を行います。

• JAGS（Just another Gibbs sampler）は、マルコフ連鎖モンテカルロ法を用いてベイズ階層モデルを解析するためのGPLプログラムです。

• Churchは、確率プログラムとして指定された任意の分布に対してギブス推論を実行するためのフリーソフトウェアです。

• PyMCは、一般的な確率的グラフィカルモデルのベイズ学習のためのオープンソースのPythonライブラリです。
• Turingは、確率的プログラミングを用いたベイズ推論のためのオープンソースのJuliaライブラリです。

• ビショップ、クリストファー M. (2006)、パターン認識と機械学習、スプリンガー、ISBN 978-0-387-31073-2
• ボルスタッド、ウィリアム・M. (2010)、『計算ベイズ統計学の理解』、ジョン・ワイリーISBN 978-0-470-04609-8
• Casella, G. ; George, EI (1992). "ギブスサンプラーの説明". The American Statistician . 46 (3): 167. CiteSeerX 10.1.1.554.3993 . doi : 10.2307/2685208 . JSTOR 2685208 .  （基本的な概要と多数の参考文献が含まれています。）
• Gelfand, Alan E.; Smith, Adrian FM (1990)、「周辺密度を計算するためのサンプリングベースのアプローチ」、Journal of the American Statistical Association、85 (410): 398–409、doi : 10.2307/2289776、JSTOR 2289776、MR 1141740  
• Gelman, A.、Carlin JB、Stern HS、Dunson D.、Vehtari A.、Rubin DB (2013)、ベイズデータ分析、第3版。ロンドン：Chapman & Hall。
• Levin, David A.; Peres, Yuval ; Wilmer, Elizabeth L. (2008)、「マルコフ連鎖と混合時間」、アメリカ数学会。
• Robert, CP; Casella, G. (2004), Monte Carlo Statistical Methods (second edition), Springer-Verlag.