吉布斯采样（Gibbs采样）

目录

MCMC(一)蒙特卡罗方法 https://www.cnblogs.com/emanlee/p/12356492.html

MCMC(二)马尔科夫链 https://www.cnblogs.com/emanlee/p/12357341.html

MCMC(三)MCMC采样和M-H采样 https://www.cnblogs.com/emanlee/p/12358022.html

MCMC(四)Gibbs采样  https://www.cnblogs.com/emanlee/p/12358194.html

 

 

 

 

 

 

 

 

 

 

 

 

 

import math
import random
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.stats import multivariate_normal

samplesource = multivariate_normal(mean=[5,-1], cov=[[1,1],[1,4]])

def p_ygivenx(x, m1, m2, s1, s2):
    return (random.normalvariate(m2 + rho * s2 / s1 * (x - m1), math.sqrt((1 - rho ** 2) * (s2**2))))

def p_xgiveny(y, m1, m2, s1, s2):
    return (random.normalvariate(m1 + rho * s1 / s2 * (y - m2), math.sqrt((1 - rho ** 2) * (s1**2))))

N = 5000
K = 20
x_res = []
y_res = []
z_res = []
m1 = 5
m2 = -1
s1 = 1
s2 = 2

rho = 0.5
y = m2

for i in range(N):
    for j in range(K):
        x = p_xgiveny(y, m1, m2, s1, s2)
        y = p_ygivenx(x, m1, m2, s1, s2)
        z = samplesource.pdf([x,y])
        x_res.append(x)
        y_res.append(y)
        z_res.append(z)

num_bins = 50
plt.hist(x_res, num_bins, normed=1, facecolor='green', alpha=0.5)
plt.hist(y_res, num_bins, normed=1, facecolor='red', alpha=0.5)
plt.title('Histogram')
plt.show()

 

 

 

然后我们看看样本集生成的二维正态分布，代码如下：

fig = plt.figure()
ax = Axes3D(fig, rect=[0, 0, 1, 1], elev=30, azim=20)
ax.scatter(x_res, y_res, z_res,marker='o')
plt.show()

 

 

 

 

 

 

from

https://www.cnblogs.com/pinard/p/6645766.html

 

REF

https://blog.csdn.net/arcers/article/details/88732639

https://www.cnblogs.com/xbinworld/p/4266146.html

https://wenku.baidu.com/view/d57107bff9c75fbfc77da26925c52cc58bd690ee.html

https://www.plob.org/article/5061.html

http://blog.sciencenet.cn/blog-255662-542389.html

https://blog.csdn.net/xc_xc_xc/article/details/76257344