gibbs采样代码python多维gibbs采样实例
更新时间：2023-10-03前言
在贝叶斯统计中，目标是寻找一个后验分布，即给定观测值后参数的分布。然而，由于许多模型的后验分布无法解析计算，因此需要使用MCMC方法来获得后验分布的近似。 Gibbs采样是一种非常常见且受欢迎的MCMC方法之一。 Gibbs采样方法通过逐个条件子分量的抽样替代所有分量同时抽样的策略，这对于许多问题来说是非常有用的。在这篇文章中，我们将介绍如何使用Python编写Gibbs采样算法。同时，我们将提供一个多维Gibbs采样的例子来说明如何使用该算法。代码示例 1
Gibbs采样算法的基本原理是从后验分布的条件分布中抽样，而不是从全体后验分布中抽样。
def gibbs_sampling(posterior, initial_value, num_samples):
"""
Gibbs采样算法实现函数，输出指定数量的样本。
参数：
posterior : 函数类型，用于计算后验分布的对数，参数为关于待抽样参数的数组。
initial_value : 初始值列表。
num_samples : 所需样本数。
返回：
samples : num_samples数量的样本列表。
"""
# 获取参数数量
num_params = len(initial_value)
# 初始化样本空数组
samples = np.zeros((num_samples, num_params))
# 初始化起始值
current_sample = initial_value
# 开始采样
for i in range(num_samples):
# 在后验分布的条件概率分布中抽样
for j in range(num_params):
# 从条件概率分布中抽样
current_sample[j] = np.random.normal(loc=current_sample[:j].sum()+current_sample[j+1:].sum(),
scale=posterior(current_sample[:j]+current_sample[j+1:]+[0])/posterior(current_sample))
# 将采样的样本存储
samples[i,:] = current_sample
return samples
代码示例 2
下面我们来看一下如何使用Gibbs采样算法来进行多维抽样。我们将使用Gibbs采样算法来从二维“双月亮”分布中抽样。先导入一些必要的包：
import matplotlib.pyplot as plt import numpy as np from scipy.stats import multivariate_normal, norm
然后，我们可以定义双月亮分布的概率密度函数：
def bimodal_pdf(x, y):
return .5*(multivariate_normal.pdf([x, y], mean=[1,0], cov=[[2,1],[1,2]])
+ multivariate_normal.pdf([x, y], mean=[-1,0], cov=[[2,1],[1,2]]))
接下来，我们可以定义一个gibbs_sampler函数来使用Gibbs采样返回指定数量的样本：
def gibbs_sampler(num_samples=100, x_start=0, y_start=0):
samples = np.zeros((num_samples,2))
x, y = x_start, y_start
for i in range(num_samples):
x = norm.rvs(loc=y, scale=1)
y = norm.rvs(loc=x, scale=1)
samples[i,:] = [x,y]
return samples
最后，我们调用gibbs_sampler函数，采样num_samples个样本，绘制结果图表：
num_samples = 5000
samples = gibbs_sampler(num_samples=num_samples, x_start=-2, y_start=3)
plt.scatter(samples[:,0], samples[:,1], alpha=.4, s=3)
plt.xlim(-5,5)
plt.ylim(-5,5)
plt.title('Gibbs Sampler for Bimodal Distribution')
plt.show()
代码示例 3
除此之外，我们也可以用多个Gibbs采样函数实现多维数据的采样。下面是一个多元高斯分布的Gibbs采样实例：
def gaussian_sampler(num_samples=10000, mu=[0,0], sigma=[[1,0],[0,1]]):
"""
用Gibbs采样进行多元高斯分布的抽样。
参数：
mu : 指定高斯分布的均值。
sigma : 指定高斯分布的协方差矩阵。
返回：
返回一个两维的样本数组。
"""
x_initial = np.zeros(2)
x = np.zeros(2)
samples = np.zeros((num_samples, 2))
for i in range(num_samples):
x[0] = np.random.normal(mu[0] + sigma[0][1]*x_initial[1]/sigma[1][1],
np.sqrt(sigma[0][0] - sigma[0][1]*sigma[1][0]/sigma[1][1]))
x[1] = np.random.normal(mu[1] + sigma[1][0]*x[0]/sigma[0][0],
np.sqrt(sigma[1][1] - sigma[1][0]*sigma[0][1]/sigma[0][0]))
samples[i,:] = x
x_initial = x
return samples
