【Acceptance—Rejection Sampling】

import seaborn as sns
import random
import math
import matplotlib.pyplot as plt
import numpy as np

sns.set_style('darkgrid')
plt.rcParams['figure.figsize'] = (12, 8)

power = 4
t = 0.4
# 被采样函数f(x)的分母
sum_ = (math.pow(1-t, power+1) - math.pow(-t, power+1)) / (power + 1)
# 采样点
x = np.linspace(0, 1, 101)
# 计算满足c*g(x) >= f(x)的c的最小值
#观察图形可知，在自变量x∈[0，1]的范围内
# x取1时c取最小值，此时c=f(1)/g(1)=7.36

# 目标概率密度函数f(x)
f = lambda x: (x-t)**power/sum_
# 辅助概率密度函数g(x)
g = lambda x: 1
# x取1时c取最小值，此时c=f(1)/g(1)=7.36
c = f(1)/g(1)
# print(c)
cc = [c for xi in x]
plt.plot(x, cc, '--', label='c*g(x)')
y = [f(xi) for xi in x]
plt.plot(x, y, label='f(x)')


# 采样1000000个点
samples = []
for i in range(1000000):
    # -- 接受拒绝采样方法 --
    Y = random.uniform(0, 1)
    U = random.uniform(0, 1)
    if U*c*g(Y) <= f(Y):
        samples.append(Y)
    # --------------------

plt.hist(samples, bins=100, stacked=True, density=True, label='sampling')
plt.legend()
plt.show()