ギブスサンプリングをpythonで実装してみた.

正規分布をギブスサンプリングによって生成するコードを書いてみました.
数式による説明は余裕ができたら書きます.(更新します)

コードを実行すると以下のようなアニメーションが再生される.正規乱数を生成しているので実行するごとに違う結果が出る



参考文献
計算統計 2 マルコフ連鎖モンテカルロ法とその周辺 (統計科学のフロンティア 12)


import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

# ボックスミューラー法による実装.

b = 0.8 # 相関項
x = 10.0
y = 3.0
times = 200
im = []
ims = []

fig = plt.figure()

for i in range(1, times):
    u1 = np.random.rand()
    u2 = np.random.rand()
    x_old = x
    x = np.sqrt(-2*np.log(u1))*np.cos(2*np.pi*u2) + b*y
    im = im + plt.plot([x_old, x], [y, y], color="blue")
    y_old = y
    y = np.sqrt(-2*np.log(u1))*np.sin(2*np.pi*u2) + b*x
    im = im + plt.plot([x, x], [y_old, y], color="blue")
    ims.append(im)

ani = animation.ArtistAnimation(fig, ims, interval=400)

# アニメーションのgifを作成する.(しない場合は必要ない)
ani.save("output.gif", writer="imagemagick")
# アニメーションを再生する
plt.show()

コメントを残す

メールアドレスが公開されることはありません。

このサイトはスパムを低減するために Akismet を使っています。コメントデータの処理方法の詳細はこちらをご覧ください