diagnostics

クラス: HamiltonianSampler

マルコフ連鎖モンテカルロの診断情報

構文

tbl = diagnostics(smp,chains)
tbl = diagnostics(smp,chains,'MaxLag',maxlag)


説明

tbl = diagnostics(smp,chains) は、chains 内の連鎖についてマルコフ連鎖モンテカルロの診断情報を返します。

tbl = diagnostics(smp,chains,'MaxLag',maxlag) では、有効な標本サイズの計算に使用する自己相関ラグの最大数を指定します。

入力引数

smp — ハミルトニアン モンテカルロ サンプラー

HamiltonianSampler オブジェクト

ハミルトニアン モンテカルロ サンプラー。HamiltonianSampler オブジェクトとして指定します。

サンプラーを作成するには、関数 hmcSampler を使用します。

chains — MCMC 連鎖

行列 | cell 配列

MCMC 連鎖。次のいずれかを指定します。

• 行列。各行は標本、各列はパラメーターです。

• 行列の cell 配列。連鎖 chains{i} は、各行が標本、各列がパラメーターである行列です。

パラメーターの個数 (つまり行列の列数) は、smp サンプラーの StartPoint プロパティの要素数に等しくなければなりません。

maxlag — 自己相関ラグの最大数

100 (既定値) | 正の整数

有効な標本サイズを計算するための自己相関ラグの最大数。正の整数を指定します。

有効な標本サイズの計算では、標本数が maxlag より多い chains 内の各連鎖について 1,2,...,maxlag 個のラグが使用されます。

標本数が maxlag 個以下である連鎖については、N_i - 1 個のラグが計算で使用されます。N_i は連鎖 i の標本数です。

例: 'MaxLag',50

出力引数

tbl — MCMC の診断情報

table

MCMC の診断情報。chains 内のすべての連鎖を使用して計算され、次の列があるテーブルとして返されます。

列	説明
Name	変数名
Mean	事後平均推定値
MCSE	モンテカルロ標準誤差の推定値 (事後平均推定値の標準偏差)
SD	事後標準偏差の推定値
Q5	周辺事後分布の 5 番目の分位数の推定値
Q95	周辺事後分布の 95 番目の分位数の推定値
ESS	事後平均推定値の有効な標本サイズ有効な標本サイズが大きくなると、結果がより正確になります。標本が独立している場合、有効な標本サイズは標本数に等しくなります。
RHat	Gelman-Rubin 収束統計量。目安として、RHat の値が 1.1 未満の場合は、目的の分布に連鎖が収束したと解釈できます。RHat が 1.1 より大きい変数がある場合は、より多くのモンテカルロ標本を抽出してください。

例

マルコフ連鎖モンテカルロの診断情報の計算

ハミルトニアン モンテカルロ (HMC) サンプラーを使用して MCMC 連鎖を作成し、MCMC の診断情報を計算します。

はじめに、多変量正規対数確率密度およびその勾配を返す関数を MATLAB® パス上に保存します。この例では、normalDistGrad という名前の関数を使用します。定義は例の終わりにあります。次に、関数 hmcSampler の入力引数 logpdf を定義する引数を使用して、この関数を呼び出します。

means = [1;-2;2];
standevs = [1;2;0.5];
logpdf = @(theta)normalDistGrad(theta,means,standevs);

開始点を選択します。HMC サンプラーを作成し、そのパラメーターを調整します。

startpoint = randn(3,1);
smp = hmcSampler(logpdf,startpoint);
smp = tuneSampler(smp);

いくつかの独立した連鎖を使用して、事後密度から標本を抽出します。各連鎖について、ランダムに分布している異なる開始点を選択します。マルコフ連鎖の先頭から破棄するバーンイン標本の数と、バーンイン後に生成する標本の数を指定します。

NumChains  = 4;
chains     = cell(NumChains,1);
Burnin     = 500;
NumSamples = 1000;
for c = 1:NumChains
    chains{c} = drawSamples(smp,'Burnin',Burnin,'NumSamples',NumSamples,...
        'Start',randn(size(startpoint)));
end

MCMC の診断情報を計算し、結果を表示します。means 内の真の平均をテーブル MCMCdiagnostics の Mean という名前の列と比較します。真の事後平均は、推定事後平均からモンテカルロ標準誤差 (MCSE) の数倍以内の範囲にあります。HMC サンプラーは真の平均を正確に復元しています。同様に、SD 列の推定標準偏差は standev の真の標準偏差に非常に近くなっています。

MCMCdiagnostics = diagnostics(smp,chains)

MCMCdiagnostics=3×8 table
     Name      Mean        MCSE        SD          Q5        Q95       ESS      RHat
    ______    _______    ________    _______    ________    ______    ______    ____

    {'x1'}     1.0038    0.016474    0.96164    -0.58601     2.563    3407.4     1  
    {'x2'}    -2.0435    0.034933      1.999     -5.3476    1.1851    3274.5     1  
    {'x3'}     1.9957    0.008209    0.49693      1.2036    2.8249    3664.5     1  

means

means = 3×1

     1
    -2
     2

standevs

standevs = 3×1

    1.0000
    2.0000
    0.5000

関数 normalDistGrad は、startpoint と同じ長さの列ベクトルまたはスカラーとして指定された、Mu に格納されている平均および Sigma に格納されている標準偏差を使用して、多変量正規確率密度の対数を返します。2 番目の出力引数は、対応する勾配です。

function [lpdf,glpdf] = normalDistGrad(X,Mu,Sigma)
Z = (X - Mu)./Sigma;
lpdf = sum(-log(Sigma) - .5*log(2*pi) - .5*(Z.^2));
glpdf = -Z./Sigma;
end

ヒント

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

バージョン履歴

R2017a で導入

参考

関数

• hmcSampler

クラス

• HamiltonianSampler

トピック

• ハミルトニアン モンテカルロの使用によるベイズ線形回帰