棄却法の例

前回の記事で紹介した棄却法を実際にやってみましょう。

今回の例はベータ分布です。(とりうる値の範囲が有限のものの方が適用しやすいので)
とりあえず、$B(2,3)$でやってみましょう。
確率密度関数は区間$x\in[0,1]$の範囲では次の式で表されます。(それ以外の$x$に対しては$0$です)。
$$
f(x) = \frac{x(1-x)^2}{B(2, 3)} = 12x(1-x)^2.
$$
この関数は$x=1/3$で最大値$f(1/3)=16/9$をとります。
(単純な式なので微分してすぐに確認できます。)

さて、実際にプログラムで実行してみたのがこちらです。
念のためですが、今回の目的は前回の記事で紹介したアルゴリズムで目的とする乱数が得られることを確認することです。
単にベータ分布に従う乱数が必要な場合は、scipyのrvs関数を使いましょう。

# ベータ分布の確率密度関数
f = beta.freeze(a=2, b=3)


def beta_rejection_sampling():
    while True:
        u, v = uniform().rvs(2)
        x = u  # 今回は 0<= x <= 1 なのでu をそのまま使用
        y = 16/9 * v
        if y <= f.pdf(x):
            return x


# 棄却法で10000データ生成する
data = [beta_rejection_sampling() for _ in range(10000)]

# 生成された乱数のヒストグラムと、確率密度関数を可視化
fig = plt.figure(facecolor="w")
ax = fig.add_subplot(111)
ax.hist(data, bins=100, density=True, label="棄却法により生成")
ax.set_xlim([0, 1])
x = np.linspace(0, 1, 100)
ax.plot(x, f.pdf(x), label="確率密度関数")
ax.legend()
plt.show()

出力がこちら。

しっかり機能していますね。