前面文章MCMC和Gibbs Sampling介绍了MCMC采用算法,这里是其编程实例

In [2]:
import random
import math
import numpy as np
from scipy.stats import norm
import scipy.special as ss
import matplotlib.pyplot as plt

%matplotlib inline


def mcmc(N_hops,p):
    states = []
    cur = random.uniform(0,1)
    for i in range(0,N_hops):
        states.append(cur)
        nextf = norm.rvs(loc=cur)
        alpha = min(p(nextf)/p(cur),1) # 计算接受概率
        u = random.uniform(0, 1)
        if u < alpha:
            cur = nextf
    return states[-1000:] # 返回进入平稳分布后的1000个状态

M-H采样正态分布

In [3]:
# 高斯分布
def norm_dist_prob(theta):
    y = norm.pdf(theta, loc=3, scale=2)
    return y

pi = mcmc(100000,norm_dist_prob)
plt.scatter(pi, norm.pdf(pi, loc=3, scale=2),label="Real Distribution")
plt.hist(pi, 50, density=1, facecolor='red', alpha=0.7,label="Simulated_MCMC")
plt.legend()
plt.show()

M-H采样Beta分布

In [4]:
# Beta分布概率密度函数
def beta(x):
    a=0.5
    b=0.6
    return (1.0 / ss.beta(a,b)) * x**(a-1) * (1-x)**(b-1)

Ly = []
Lx = []
i_list = np.mgrid[0:1:100j]
for i in i_list:
    Lx.append(i)
    Ly.append(beta(i))
    
plt.plot(Lx, Ly, label="Real Distribution")
plt.hist(mcmc(100000,beta),density=1,bins=25, histtype='step',label="Simulated_MCMC")
plt.legend()
plt.show()
/usr/local/lib/python3.6/site-packages/ipykernel_launcher.py:5: RuntimeWarning: divide by zero encountered in double_scalars
  """
/usr/local/lib/python3.6/site-packages/ipykernel_launcher.py:5: RuntimeWarning: invalid value encountered in double_scalars
  """

Gibbs Sampling

In [5]:
from mpl_toolkits.mplot3d import Axes3D
from scipy.stats import multivariate_normal

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

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

# 定义二维高斯分布
samplesource = multivariate_normal(mean=[5,-1], cov=[[1,0.5],[0.5,2]])

N = 5000*20
m1 = 5
m2 = -1
s1 = 1
s2 = 2
rho = 0.5
y = m2
x_res = []
y_res = []
z_res = []

for _ in range(N):
    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, density=1, facecolor='green', alpha=0.5)
plt.hist(y_res, num_bins, density=1, facecolor='red', alpha=0.5)
plt.title('Histogram')
plt.show()
In [6]:
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()
In [ ]:
 
posted @ 2019-01-01 19:56:46
评论加载中...

发表评论