逆機率取樣-接受拒絕取樣-MCMC取樣

import numpy as np
import scipy
from matplotlib import pyplot as plt


def pdf(x):
    if 0 <= x < 0.25:
        return 8 * x
    elif 0.25 <= x < 1:
        return 8 / 3 - 8 / 3 * x
    else:
        return 0


def cdf(x):
    if x < 0:
        return 0
    elif 0 <= x < 0.25:
        return 4 * x * x
    elif 0.25 <= x < 1:
        return 8 / 3 * x - 4 / 3 * x * x - 1 / 3
    else:
        return 1


def cdf_reverse(x):
    if x < 0 or x > 1:
        return None
    elif 0 <= x < 0.25:
        return x ** 0.5 / 2
    else:
        return 1 - (3 * (1 - x)) ** 0.5 / 2


# 逆CDF取樣
def reverse_sample():
    #########################################################
    # p()是目標分佈PDF
    # F^-1()是目標分佈CDF的反函式
    # u()是0-1均勻分佈
    # 演算法：
    #     從u()中隨機抽取u
    #     返回F^-1(u)
    #########################################################
    u = np.random.uniform(0, 1)
    return cdf_reverse(u)


# 接受拒絕取樣
def reject_sample():
    #########################################################
    # p()是目標分佈
    # q()取0-1均勻分佈 並且使得kq()要覆蓋p()
    # u()是0-1均勻分佈
    # 演算法：
    #     從q()中隨機抽取x0，從u()隨機抽取u
    #     如果 u <= p(x0) / kq(x0),則接受x0,反之則拒絕
    #########################################################
    x = np.random.uniform(0, 1)
    k = 2
    u = np.random.uniform(0, 1)
    if u <= pdf(x) / k:
        return x
    else:
        return None


# MCMC MH演算法
def mcmc_sample(x):
    #########################################################
    # p()是目標分佈
    # q()隨機選擇一個分佈，選擇正態分佈
    # u()是0-1均勻分佈
    # 演算法：
    #     初始化x，並從q()隨機抽取x_star
    #     計算alpha=min(1, p(j)*q(j,i) / p(i)q(i,j))
    #     從u()隨機抽取u
    #     如果 u < alpha,則x = x_star,反之則x = x
    #########################################################
    x_star = np.random.normal(x, 1)
    num = pdf(x_star) * scipy.stats.norm(x_star, 1).pdf(x)
    den = pdf(x) * scipy.stats.norm(x, 1).pdf(x_star)
    alpha = min(1, num / den)
    u = np.random.uniform(0, 1)
    if u < alpha:
        return x_star
    else:
        return x


def plot(x, y, samples):
    plt.hist(samples, color="green", density=True)
    plt.plot(x, y, color="red")
    plt.show()


def mean(samples):
    return np.mean(samples)


def get_reverse_sample():
    # 隨機取樣10000次
    samples = []
    for _ in range(10000):
        val = reverse_sample()
        if val:
            samples.append(val)
    return samples


def get_reject_sample():
    # 隨機取樣10000次
    samples = []
    for _ in range(10000):
        val = reject_sample()
        if val:
            samples.append(val)
    return samples


def get_mcmc_sample():
    # 隨機取樣10000次
    samples = []
    x = 0.1
    for _ in range(10000):
        x = mcmc_sample(x)
        samples.append(x)
    return samples[1000:]


def main():
    x = np.arange(0, 1, 0.01)
    y = [pdf(i) for i in x]
    samples = get_mcmc_sample()
    plot(x, y, samples)

    # 理論均值 = 33/72 - 1/24
    print("理論均值:", 33 / 72 - 1 / 24)
    print("取樣均值:", mean(samples))


if __name__ == '__main__':
    main()