import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
# 目标分布：标准正态分布
def target_distribution(x):
return np.exp(-0.5 * x**2) / np.sqrt(2 * np.pi)
# 建议分布：以当前状态为中心的正态分布
def proposal_distribution(x, sigma=1):
return np.random.normal(x, sigma)
# Metropolis-Hastings算法
def metropolis_hastings(target, proposal, n_samples, initial_state, sigma=1):
samples = [initial_state]
current_state = initial_state
accepted = 0
for _ in range(n_samples):
# 生成候选样本
candidate = proposal(current_state, sigma)
# 计算接受概率
acceptance_prob = min(1, target(candidate) / target(current_state))
# 接受或拒绝
if np.random.random() < acceptance_prob:
samples.append(candidate)
current_state = candidate
accepted += 1
else:
samples.append(current_state)
acceptance_rate = accepted / n_samples
print(f”接受率: {acceptance_rate:.2f}”)
return samples[1000:] # 舍弃前1000个样本作为burn-in期
# 运行MH算法
samples = metropolis_hastings(target_distribution, proposal_distribution, 10000, 0, sigma=1)
# 绘制结果
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.plot(samples)
plt.title(‘样本序列’)
plt.subplot(1, 2, 2)
sns.histplot(samples, kde=True, stat=’density’)
x = np.linspace(-4, 4, 1000)
plt.plot(x, target_distribution(x), ‘r-‘, label=’目标分布’)
plt.title(‘样本分布’)
plt.legend()
plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
# 二元正态分布的参数
mu = np.array([0, 0])
cov = np.array([[1, 0.8], [0.8, 1]])
# 条件分布的参数
def conditional_params(x, index):
# 对于二元正态分布，条件分布也是正态分布
# X1|X2=x2 ~ N(μ1 + ρ*σ1/σ2*(x2-μ2), σ1²*(1-ρ²))
# X2|X1=x1 ~ N(μ2 + ρ*σ2/σ1*(x1-μ1), σ2²*(1-ρ²))
rho = cov[0, 1] / np.sqrt(cov[0, 0] * cov[1, 1])
if index == 0: # X1|X2
mean = mu[0] + rho * np.sqrt(cov[0, 0] / cov[1, 1]) * (x[1] – mu[1])
var = cov[0, 0] * (1 – rho**2)
else: # X2|X1
mean = mu[1] + rho * np.sqrt(cov[1, 1] / cov[0, 0]) * (x[0] – mu[0])
var = cov[1, 1] * (1 – rho**2)
return mean, var
# Gibbs抽样
def gibbs_sampling(mu, cov, n_samples, initial_state):
samples = [initial_state]
current_state = initial_state
for _ in range(n_samples):
# 更新X1
mean1, var1 = conditional_params(current_state, 0)
x1 = np.random.normal(mean1, np.sqrt(var1))
# 更新X2
mean2, var2 = conditional_params([x1, current_state[1]], 1)
x2 = np.random.normal(mean2, np.sqrt(var2))
current_state = [x1, x2]
samples.append(current_state)
return np.array(samples[1000:]) # 舍弃前1000个样本作为burn-in期
# 运行Gibbs抽样
samples = gibbs_sampling(mu, cov, 10000, [0, 0])
# 绘制结果
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.plot(samples[:, 0], samples[:, 1], ‘b-‘, alpha=0.3)
plt.scatter(samples[::100, 0], samples[::100, 1], c=’red’, s=10)
plt.title(‘样本轨迹’)
plt.subplot(1, 2, 2)
sns.kdeplot(x=samples[:, 0], y=samples[:, 1], cmap=”Blues”, shade=True)
plt.title(‘样本密度’)
plt.tight_layout()
plt.show()

import pymc3 as pm
import numpy as np
import matplotlib.pyplot as plt
# 生成一些模拟数据
np.random.seed(42)
true_mu = 5
true_sigma = 2
data = np.random.normal(true_mu, true_sigma, size=100)
# 使用HMC进行贝叶斯推断
with pm.Model() as model:
# 定义先验分布
mu = pm.Normal(‘mu’, mu=0, sigma=10)
sigma = pm.HalfNormal(‘sigma’, sigma=10)
# 定义似然函数
likelihood = pm.Normal(‘likelihood’, mu=mu, sigma=sigma, observed=data)
# 使用HMC进行采样
trace = pm.sample(1000, tune=1000, step=pm.HamiltonianMC())
# 查看结果
pm.summary(trace)
# 绘制后验分布
pm.plot_posterior(trace, var_names=[‘mu’, ‘sigma’])
plt.show()
# 绘制轨迹图
pm.plot_trace(trace, var_names=[‘mu’, ‘sigma’])
plt.show()