読者です 読者をやめる 読者になる 読者になる

どせいたんさき。

ナスダヨー

Metropolis-Hastings 法で乱数生成をやってみた

目的

任意の分布に従う乱数を Metropolis-Hastings 法で生成する. Metropolis-Hastings 法がどのようにして乱数列を生成するのかを視覚的に理解したい.

手法

Octave を用いて以下のコードを実行.のこぎりの刃型の分布関数に従う乱数を 10,000 個生成する.初期の分布としては標準偏差が 3.0 の正規分布を用いた.for 文ではこの分布を Metropolis-Hastings 法にしたがって更新していく.mh(x0,func,xgen) は数列を更新するための関数である.定義は以下の gist を参照.本当は分布が十分に収束したかどうかを調べる必要があるが今回は省略.とりあえず初期の挙動だけが分かればいいということで 200 ステップだけ実行した.

x0 = 3*randn(1,1e4);

func = @(x) max(0, (1-0.5*x).*(x>-1));
xgen = @(x) 3*randn(size(x));

for n=2:200; X(n,:) = mh(X(n-1,:),func,xgen); endfor

結果

とりあえず最初の 40 ステップだけを抜き出してヒストグラムの変化をアニメーションにしてみた.初期状態では数列は正規分布でばらついているが,ステップを進めるとヒストグラムはターゲットである三角形の分布に収束していく.