Rejection Sampling을 통한 Beta 분포 확률변수 생성

Rejection sampling을 사용하여 Beta (α, β); α=6, β=4 분포로부터의 확률변수 값 1000개를 모의 생성하여 평균, 표준편차 및 히스토그램을 나타내는 파이썬 코드 실습

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import beta

# hyper parameters
num_samples = 1000
a, b = 6, 4
x = np.linspace(0, 1, num_samples)

# target / proposal distribution
target_dist = beta.pdf(x, a, b)
proposal_dist = [np.random.uniform(0,1)*3 for _ in range(num_samples)]  # 3은 proposal_dist의 상수

# 판별을 위한 list
accepted, rejected = [], []

# random point가 bata분포 안에 들어가는지 아닌지 판별
for i in range(num_samples):
    accepted.append((x[i], proposal_dist[i])) if proposal_dist[i] <= target_dist[i] else rejected.append((x[i], proposal_dist[i]))

# 판별 결과에 대한 display
plt.plot(x, target_dist, color='blue')
plt.plot([x[0] for x in accepted], [x[1] for x in accepted], 'ro', color='g')
plt.plot([x[0] for x in rejected], [x[1] for x in rejected], 'ro', color='r')
plt.show()

# beta 분포에 해당하는 값들의 평균, 표준편차, 히스토그램
beta_avg = np.mean([x[1] for x in accepted])
beta_std = np.std([x[1] for x in accepted])
print(f"Avg : {beta_avg:.4f}, Std : {beta_std:.4}")
plt.hist([x[1] for x in accepted], bins=50)
plt.show()