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
def reverse_sample():
u = np.random.uniform(0, 1)
return cdf_reverse(u)
def reject_sample():
x = np.random.uniform(0, 1)
k = 2
u = np.random.uniform(0, 1)
if u <= pdf(x) / k:
return x
else:
return None
def mcmc_sample(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():
samples = []
for _ in range(10000):
val = reverse_sample()
if val:
samples.append(val)
return samples
def get_reject_sample():
samples = []
for _ in range(10000):
val = reject_sample()
if val:
samples.append(val)
return samples
def get_mcmc_sample():
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)
print("理论均值:", 33 / 72 - 1 / 24)
print("采样均值:", mean(samples))
if __name__ == '__main__':
main()
【推荐】国内首个AI IDE,深度理解中文开发场景,立即下载体验Trae
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步