定义

输入:目标概率分布的密度函数$p(x)$,函数$f(x)$

输出:$p(x)$的随机样本$x_{m+1},x_{m+2},\cdots ,x_n$,函数样本均值$f_{mn}$;

参数:收敛步数$m$,迭代步数$n$。

(1)初始化。给出初始样本$x^0=(x_1^{(0)},x_2^{(0)},\cdots ,x_k^{(0)}{)}^T$。

(2)对i循环执行

设第$(i-1)$次迭代结束时的样本为$x^{(i-1)}=(x_1^{(i-1)},x_2^{(i-1)},\cdots ,x_k^{(i-1)}{)}^T$,则第$i$次迭代进行如下几步操作:

$\{\begin{array}{l}(1)由满足条件分布p(x_1|x_2^{(i)},\cdots ,x_k^{(i-1)})抽取x_1^{(i)}\\ ⋮\\ (j)由满足条件分布p(x_j|x_1^{(i)},\cdots ,x_{j-1}^{(i)},x_{j+1}^{(i)},\cdots ,x_k^{(i-1)})抽取x_j^{(i)}\\ ⋮\\ (k)由满足条件分布p(x_k|x_1^{(i)},\cdots ,x_{(k-1)}^{(}i)抽取x_k^{(i)}\end{array}$

得到第$i$次迭代值$x^{(i)}=(x_1^{(i)},x_2^{(i)},\cdots ,x_k^{(i)}{)}^T$
(3)得到样本集合
{ x ( m + 1 ) , x ( m + 2 ) , ⋯ , x ( n ) } \{ x^{(m+1)},x^{(m+2)},\cdots,x^{(n)} \} {x(m+1),x(m+2),⋯,x(n)}
(4)计算
$$f_{mn}=\frac{1}{n-m}\sum _{i=m+1}^{n}f(x^{(i)})$$

import numpy as np
from matplotlib import pyplot as plt
from scipy.stats import gaussian_kdedef gibbs_sampling(dim, conditional_sampler, x0=None, burning_steps=1000, max_steps=10000, epsilon=1e-8, verbose=False):"""给定一个从p（x_j|x_1，x_2，…x_n）采样的条件采样器返回样本x~p的列表，其中p是条件分布的原始分布。x0是x的初始值。如果未指定，则将其设置为零向量。条件采样器将（x，j）作为参数"""x = np.zeros(dim) if x0 is None else x0samples = np.zeros([max_steps - burning_steps, dim])for i in range(max_steps):for j in range(dim):x[j]  = conditional_sampler(x, j)if verbose:print("New value of x is", x_new)if i >= burning_steps:samples[i - burning_steps] = xreturn samplesif __name__ == '__main__':i = 0def demonstrate(dim, p, desc, **args):samples = gibbs_sampling(dim, p, **args)z = gaussian_kde(samples.T)(samples.T)plt.scatter(samples[:, 0], samples[:, 1], c=z, marker='.')plt.plot(samples[: 100, 0], samples[: 100, 1], 'r-')plt.title(desc)plt.savefig(str(i)+".png")plt.show()# example 1:mean = np.array([2, 3])covariance = np.array([[1, 0],[0, 1]])covariance_inv = np.linalg.inv(covariance)det_convariance = 1def gaussian_sampler1(x, j):return np.random.normal()demonstrate(2, gaussian_sampler1, "Gaussian distribution with mean of 2 and 3")# example 2:mean = np.array([2, 3])covariance = np.array([[1, 0],[0, 1]])covariance_inv = np.linalg.inv(covariance)det_convariance = 1def gaussian_sampler2(x, j):if j == 0:return np.random.normal(2)else:return np.random.normal(3)demonstrate(2, gaussian_sampler2, "Gaussian distribution with mean of 2 and 3")# example 3:def blocks_sampler(x, j):sample = np.random.random()if sample > .5:sample += 1.return sampledemonstrate(2, blocks_sampler, "Four blocks")# example 4:def blocks_sampler(x, j):sample = np.random.random()if sample > .5:sample += 100.return sampledemonstrate(2, blocks_sampler, "Four blocks with large gap.")