Main Content

mhsample

メトロポリス・ヘイスティングス標本

ページ内をすべて折りたたむ

構文

smpl = mhsample(start,nsamples,'pdf',pdf,'proppdf',proppdf, 'proprnd',proprnd)
smpl = mhsample(...,'symmetric',sym)
smpl = mhsample(...,'burnin',K)
smpl = mhsample(...,'thin',m)
smpl = mhsample(...,'nchain',n)
[smpl,accept] = mhsample(...)

説明

smpl = mhsample(start,nsamples,'pdf',pdf,'proppdf',proppdf, 'proprnd',proprnd) は、メトロポリス・ヘイスティングス アルゴリズムを使って、ターゲットの定常分布 pdf から nsamples の無作為標本を抽出します。

start は、マルコフ連鎖の開始値を含む行ベクトルです。nsamples は、生成する標本数を指定する整数です。pdf、proppdf、および proprnd は、@ を使用して作成された関数ハンドルです。proppdf は、提案分布の密度を定義し、proprnd は提案分布に対する乱数発生器を定義します。pdf と proprnd は、start と同じタイプとサイズをもつ 1 つの引数を入力として受け入れます。proppdf は、start と同じタイプとサイズをもつ 2 つの引数を入力として受け入れます。

smpl は、標本を含む列ベクトルまたは行列です。対数密度関数が選ばれた場合、'pdf' と 'proppdf' は、'logpdf' と 'logproppdf' に置き換えることが可能です。メトロポリス・ヘイスティングス アルゴリズムで使われる密度関数は、必ずしも正規化されるわけではありません。

提案分布 q(x,y) は、y が現在の点である場合、次の点として選択する x に対して確率密度を与えます。q(x|y) として記述される場合もあります。

proppdf または logproppdf が q(x,y) = q(y,x) を満たす場合、すなわち提案分布が対称の場合、mhsample は、ランダム ウォーク メトロポリス・ヘイスティングス サンプリングを実行します。proppdf または logproppdf が q(x,y) = q(x) を満たす場合、すなわち提案分布が現在の値に依存しない場合、mhsample は、独立メトロポリス・ヘイスティングス サンプリングを実行します。

smpl = mhsample(...,'symmetric',sym) は、メトロポリス・ヘイスティングス アルゴリズムを使って、ターゲットの定常分布 pdf から nsamples の無作為標本を抽出します。sym は、提案の分布が対称であるかどうかを示す論理値です。既定値は false で、非対称な提案分布に対応します。sym が true、すなわち提案分布が対称の場合、proppdf と logproppdf はオプションです。

smpl = mhsample(...,'burnin',K) は、開始点と生成されたシーケンスで省略された k 番目の点の間の値をもつマルコフ連鎖を生成します。k 番目を超える値は保持されます。k は既定の設定の幅 0 をもつ非負の整数です。

smpl = mhsample(...,'thin',m) は、生成されたシーケンスで省略された m 値からの m-1 をもつマルコフ連鎖を生成します。m は、既定値 1 をもつ正の整数です。

smpl = mhsample(...,'nchain',n) は、メトロポリス・ヘイスティングス アルゴリズムを使用して、n マルコフ連鎖を生成します。n は、既定値 1 をもつ正の整数です。smpl は、標本を含んでいる行列です。最後の次元は、独立した連鎖のインデックスを含みます。

[smpl,accept] = mhsample(...) は、提案分布の受容比である accept も返します。accept は、単一連鎖が生成された場合はスカラーで、複数連鎖が生成された場合はベクトルです。

例

すべて折りたたむ

ライブ スクリプトを開く

独立メトロポリス・ヘイスティングス サンプリングを使用して、ガンマ分布の 2 次積率を推定します。

rng default;  % For reproducibility
alpha = 2.43;
beta = 1;
pdf = @(x)gampdf(x,alpha,beta); % Target distribution
proppdf = @(x,y)gampdf(x,floor(alpha),floor(alpha)/alpha);
proprnd = @(x)sum(...
              exprnd(floor(alpha)/alpha,floor(alpha),1));
nsamples = 5000;
smpl = mhsample(1,nsamples,'pdf',pdf,'proprnd',proprnd,...
                'proppdf',proppdf);

結果をプロットします。

xxhat = cumsum(smpl.^2)./(1:nsamples)';
figure;
plot(1:nsamples,xxhat)

Figure contains an axes object. The axes object contains an object of type line.

ライブ スクリプトを開く

ランダム ウォーク メトロポリス・ヘイスティングス サンプリングを使用して、標準正規分布から標本データを生成します。

rng default  % For reproducibility
delta = .5;
pdf = @(x) normpdf(x);
proppdf = @(x,y) unifpdf(y-x,-delta,delta);
proprnd = @(x) x + rand*2*delta - delta;   
nsamples = 15000;
x = mhsample(1,nsamples,'pdf',pdf,'proprnd',proprnd,'symmetric',1);

標本データをプロットします。

figure;
h = histfit(x,50);
h(1).FaceColor = [.8 .8 1];

Figure contains an axes object. The axes object contains 2 objects of type bar, line.

バージョン履歴

R2006a で導入

参考

|

トピック