任意の確率分布に従う標本生成：棄却法

累積密度関数の逆関数から標本を発生させる方法は、一次元標本に限ってしまうが、多変量の多次元に拡張できる。
確率変数$${x: x \in {\bf R}}$$の確率密度関数$${f(x)}$$であり、密度関数の上限$${M: f(x)\le M}$$が与えられているとする。
$${{\bf R}}$$内で一様連続分布に従う乱数$${x}$$と、$${(0,M)}$$内の一様連続分布に従う乱数$${u}$$を生成する。
$${x: u\leq f(x); accepted}$$
$${x: u>f(x); rejected}$$
これを$${a=1, b=5}$$のベータ分布に従う標本作成に適応したのが以下のコードとなる。
$${f(x)=\displaystyle{\frac{(1-x)^{4}}{B(1,5)}=5(1-x)^{4}}}$$より、

import numpy as np
from numpy import random as rand
import matplotlib.pyplot as plt
rand.seed(0)
randnum=1_000_000

b=5
a=1.0
b1=rand.beta(a, b, size=randnum)
x =rand.random(randnum)
u = rand.uniform(0,5,randnum)
bb=x[u<5*(1-x)**4]
plt.hist(bb, density=True, bins='auto',alpha=0.6)
plt.hist(b1, density=True, bins='auto',alpha=1.0,histtype="step",label='numpy.rand.beta')
plt.title('Beta distribution two methods')
plt.ylabel('Probability')
plt.legend()
plt.xlabel('Data');