与えられた確率分布に従う乱数の生成法

逆関数法

累積分布関数の性質

連続確率変数$X$の累積分布関数を$F(x)=P(X \leq x)$とすると
$$P(F(X) \leq x)=P(X \leq F^{-1}(x))=F(F^{-1}(X))=x$$よって$F(X)$は$U(0,1)$に従う。逆関数法はこの応用である。

逆関数法

連続確率変数$X$の累積分布関数$F(x)$の逆関数$F^{-1}(x)$が求まる場合は、区間$[0,1]$上の一様乱数を$U$として、$Y=F^{-1}(U)$とすればよい。このとき$Y$は$X$の確率分布に従う乱数となる。
なぜなら、$Y$の累積分布関数を$G$としたときに
$$G(y)=P(Y \leq y)=P(F^{-1}(U) \leq y) = P(U \leq F(y))=F(y)$$
が成り立つからである。

逆関数法の例

$X \sim \rm{Exp}(\lambda)$ の場合、$F(x)=1-e^{-\lambda x}=u$を解くと$x=\frac{-1}{\lambda} \log (1-u)$となるので、$Y=\frac{-1}{\lambda} \log (1-U)$とすればよい。

逆関数法のコードの例

import numpy as np
import matplotlib.pyplot as plt

# 逆関数法で指数分布からサンプリングする関数
def inverse_transform_sampling(lambda_param, size):
    u = np.random.rand(size)  # 0から1の一様分布からのサンプリング
    x = -np.log(1 - u) / lambda_param  # 逆関数法で指数分布に変換
    return x

# パラメータの設定
lambda_param = 0.5  # 指数分布のパラメータ
sample_size = 1000

# 指数分布からサンプリング
samples = inverse_transform_sampling(lambda_param, sample_size)

# ヒストグラムの表示
plt.hist(samples, bins=50, density=True, alpha=0.7, color='blue', edgecolor='black')

# 理論的な確率密度関数の描画
x = np.linspace(0, max(samples), 1000)
pdf = lambda_param * np.exp(-lambda_param * x)
plt.plot(x, pdf, 'r-', lw=2)

plt.title('Inverse Transform Sampling: Exponential Distribution')
plt.xlabel('Value')
plt.ylabel('Probability Density')
plt.show()

棄却法

逆関数が複雑な時は棄却法を利用する。

$F$の確率密度関数を$f$として、その最大値を$f_{\rm{max}}$とする。このとき、確率$f(U)/f_{\rm{max}}$で$U$を保持し、$1-f(U)/f_{\rm{max}}$で棄却する。

すなわち$V$を区間$[0,1]$上の一様分布に従う確率変数としたときに$V \leq f(U)/f_{\rm{max}}$が成り立つならば$U$を保持するということである。

これを必要なサンプル数が得られるまで繰り返す。

棄却法のコードの例

import numpy as np
import matplotlib.pyplot as plt

# 目標確率分布（標準正規分布）の確率密度関数
def target_distribution(x):
    return np.exp(-0.5 * x**2) / np.sqrt(2 * np.pi)

# 提案確率分布（一様分布）の確率密度関数
def proposal_distribution(x):
    return 1.0 / np.sqrt(2 * np.pi)

# 棄却法によるサンプリング
def rejection_sampling(sample_size):
    samples = []
    for _ in range(sample_size):
        # 提案確率分布からサンプリング
        x_proposal = np.random.uniform(-3, 3)
        # 一様分布からサンプリング
        u = np.random.uniform(0, proposal_distribution(x_proposal))
        # 棄却条件を確認
        if u < target_distribution(x_proposal) / proposal_distribution(x_proposal):
            samples.append(x_proposal)
    return np.array(samples)

# パラメータの設定
sample_size = 1000

# 棄却法によるサンプリング
samples = rejection_sampling(sample_size)

# ヒストグラムの表示
plt.hist(samples, bins=50, density=True, alpha=0.7, color='blue', edgecolor='black')

# 目標確率分布の確率密度関数の描画
x = np.linspace(-3, 3, 1000)
plt.plot(x, target_distribution(x), 'r-', lw=2)

plt.title('Rejection Sampling: Normal Distribution')
plt.xlabel('Value')
plt.ylabel('Probability Density')
plt.show()