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)
ランダム ウォーク メトロポリス・ヘイスティングス サンプリングを使用して、標準正規分布から標本データを生成します。
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];
バージョン履歴
R2006a で導入
