Main Content

マルコフ連鎖サンプラーを使用する標本分布の表現

より複雑な確率分布については、一般的な疑似乱数発生法で説明されている方法より高度な方法で標本を生成する必要がある可能性があります。このような分布は、データのベイズ解析やマルコフ連鎖モンテカルロ (MCMC) シミュレーションの大規模な組み合わせ問題などで発生します。代わりの方法は、状態の分布が目標に収束する最初のバーンイン期間の後に、鎖の状態を使って乱数を生成し、目標の標本分布に等しい定常分布をもつマルコフ連鎖を作成することです。

メトロポリス・ヘイスティングス アルゴリズムの使用

メトロポリス・ヘイスティングス アルゴリズムは、定数のみから求められる分布から標本を取り出します。乱数は、提案関数と等しいかまたは比例する、確率密度関数をもつ分布から生成されます。

乱数を生成するには、次の操作を行います。

  1. 初期値 x(t) を仮定します。

  2. 提案分布 q(y|x(t)) から標本 y(t) を抽出します。

  3. r(x(t),y(t)) という確率で y(t) を次の標本 x(t + 1) として採択し、1 - r(x(t),y(t)) という確率で x(t) を次の標本 x(t + 1) として保持します。ここで、

    r(x,y)=min{f(y)f(x)q(x|y)q(y|x),1}

  4. t を t + 1 としてインクリメントし、望ましい標本数が得られるまで手順 2 と 3 を繰り返します。

関数 mhsample でメトロポリス・ヘイスティングス法を使用して、乱数を生成します。メトロポリス・ヘイスティングス アルゴリズムを使って、質の良い標本を効率よく生成するには、適切な提案分布を選択することが重要です。有効な提案分布を見つけることが難しい場合は、スライス サンプリング (slicesample) またはハミルトニアン モンテカルロ (hmcSampler) を代わりに使用します。

スライス サンプリングの使い方

効率的なメトロポリス・ヘイスティングス提案分布を見つけることが難しい例の場合、スライス サンプリング アルゴリズムは明示的な指定を必要としません。Slice Sampling アルゴリズムは、垂直方向と水平方向のステップの列を使って、密度関数の下の領域から標本を抽出します。はじめに、0 ~密度関数 f (x) で、無作為に高さを選択します。次に、選択された高さよりも高い密度の水平方向の "スライス" からサンプリングすることによって、新しい x の値を無作為に選択します。同様の Slice Sampling アルゴリズムは、多変量分布に対しても使用されます。

密度関数に比例する関数 f (x) が与えられる場合は、以下の操作で乱数を生成します。

  1. f (x) の定義域内で、初期値 x(t) を仮定します。

  2. (0, f(x(t))) から実数 y を等間隔に抽出することで、水平方向の "スライス" を S = {x: y < f(x)} として定義します。

  3. すべて、またはほとんどの "スライス" S を含む x(t) 周辺の区間 I = (L, R) を検出します。

  4. この区間内に新たな点 x(t + 1) を描きます。

  5. t → t + 1 と増加させ、目的とする標本数が得られるまで手順 2 から 4 を繰り返します。

密度の "スライス" である、区間 I = (L,R) を効率的な数値処理で見つけることができれば、スライス サンプリングは、任意の形の密度関数をもつ分布から乱数を生成できます。

関数 slicesample によりスライス サンプリング法を使用し、乱数を生成します。

ハミルトニアン モンテカルロの使用

メトリポリス・ヘイスティングとスライス サンプリングは、特に中次元問題と高次元問題で、混合が遅く定常分布に収束するまで時間がかかる MCMC 連鎖を生成する可能性があります。このような状況におけるサンプリングを高速化するには、勾配に基づくハミルトニアン モンテカルロ (HMC) サンプラーを使用します。

HMC サンプリングを使用するには、(最大で加法定数までの) log f(x) とその勾配を指定しなければなりません。数値勾配を使用できますが、サンプリングが遅くなります。すべてのサンプリング変数は制約なしでなければなりません。つまり、すべての実数 x について log f(x) とその勾配が明確に定義されます。制約がある変数をサンプリングするには、それらを制約なしの変数に変換してから HMC サンプラーを使用します。

HMC サンプリング アルゴリズムでは、ランダムな "運動量ベクトル" z を導入し、z の同時密度と "位置ベクトル" x を P(x,z) = f(x)g(z) として定義します。目標は、この同時分布からサンプリングしてから、z の値を無視すること、x の周辺分布が必要な密度 f(x) をもつことです。

HMC アルゴリズムでは、共分散行列 M ("質量行列") をもつガウス密度を z に割り当てます。

g(z)exp(12zTM1z)

そして、"エネルギー関数" を次のように定義します。

E(x,z)=logf(x)+12zTM1z=U(x)+K(z)

U(x) = - log f(x) は "位置エネルギー"、K(z) = zTM-1z/2 は "運動エネルギー" です。同時密度は P(x,z) ∝ exp {-E(x,z)} によって与えられます。

無作為標本を生成するため、HMC アルゴリズムでは以下を行います。

  1. 位置ベクトルの初期値を x と仮定します。

  2. 運動量ベクトルの標本 z ∼ g(z) を生成します。

  3. "運動方程式" を使用して、いくらかの架空の時間 τ について状態 (x, z) を新しい状態 (x’,z’) に変化させます。

    dzdτ=Ux

    dxdτ=Kz

    運動方程式を厳密に解くことができるのであれば、エネルギー (したがって密度) は一定のままです (E(x,z) = E(x’,z’))。実際には、(通常は、いわゆるリープフロッグ法を使用して) 運動方程式を数値的に解かなければならず、エネルギーは保存されません。

  4. x’ を pacc = min(1, exp{E(x,z) – E(x’,z’)}) の確率で次の標本として受け入れ、x を 1 - pacc の確率で次の標本として保持します。

  5. 必要な数の標本を生成するまで、手順 2 ~ 4 を繰り返します。

HMC サンプリングを使用するには、関数 hmcSampler を使用してサンプラーを作成します。サンプラーの作成後、MAP (最大事後確率) 点推定の計算、サンプラーの調整、標本の抽出、収束診断のチェックを行うことができます。このワークフローの例については、ハミルトニアン モンテカルロの使用によるベイズ線形回帰を参照してください。

参考

関数