逆概率采样-接受拒绝采样-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()