下面总结这么几点：

1、蒙特卡洛数值积分

2、均匀分布，Box-Muller 变换

3、Monte Carlo principle

4、接受-拒绝抽样（Acceptance-Rejection sampling)

5、重要性抽样(Importance sampling)

6、马尔科夫链，马尔科夫稳态

7、MCMC——Metropolis-Hasting算法

8、MCMC——Gibbs Sampling算法

1.

MCMC概述

从名字我们可以看出，MCMC由两个MC组成，即蒙特卡罗方法（Monte Carlo Simulation，简称MC）和马尔科夫链（Markov Chain ，也简称MC）。要弄懂MCMC的原理我们首先得搞清楚

蒙特卡罗方法

和

马尔科夫链的原理

。

gibbs简析：

gibbs采样

需要知道样本中一个属性在其它所有属性下的条件概率，然后利用这个条件概率来分布产生各个属性的样本值。

gibbs采样属于随机模拟抽样算法中的一种（一类近似求解的方法）。随机模拟的核心是对一个分布进行抽样，常用的

抽样算法

包括：

1. 接受-拒绝抽样；
2. 重要性抽样；
3. MCMC（马尔科夫链蒙特卡洛方法）方法，它包括两个非常著名的采样算法（
   
   metropolis-hasting算法
   
   和它的特例
   
   Gibbs采样算法
   
   ）（补充：MCMC方法最早由Metropolis（1954）给出，后来Metropolis的算法由Hastings改进，合称为M-H算法。M-H算法是MCMC的基础方法。由M-H算法演化出了许多新的抽样方法，包括目前在MCMC中最常用的Gibbs抽样也可以看做M-H算法的一个特例）。

2、蒙特卡洛数值积分

3、均匀分布，Box-Muller 变换

4、Monte Carlo principle

4.1. 随机投点法

4.2. 平均值法

import random

r = 1.0  # 圆半径，假设为1
n = 300000  # 总投点数
count = 0  # 落到圆内投点数

# 投点x,y的范围
x_min, x_max = -r, r
y_min, y_max = -r, r

for i in range(0, n):

    # 在 [min, max] 范围内随机生成实数
    x = random.uniform(x_min, x_max)
    y = random.uniform(y_min, y_max)

    # 落到圆内投点数+1
    if x * x + y * y <= r:
        count += 1

pi = (count / float(n)) * 4
print("pi is ", pi)

# pi is 3.142053333333333

import random

n = 300000  # 总投点数
count = 0  # 位于曲线之下的投点数

# 投点x,y的范围
x_min, x_max = 1.0, 2.0
y_min, y_max = 0.0, 8.0

for i in range(0, n):

    # 在[min, max]范围内随机生成实数
    x = random.uniform(x_min, x_max)
    y = random.uniform(y_min, y_max)

    # 位于曲线之下的投点数+1
    if x * x * x >= y:
        count += 1

# 所求的积分值即为曲线下方的面积与正方形面积8的比
result = count / float(n) * 8
print("result is ", result)

# result is  3.7480533333333335

5. 接受-拒绝采样

下面我们来证明下接受-拒绝方法采样得到的样本服从 π(x) 分布。

证明：accept x 服从 π(x) 分布，即 p(x|accept) =π(x)。

import random
import math
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np

%matplotlib inline
sns.set_style('darkgrid')
plt.rcParams['figure.figsize'] = (12, 8)


def AceeptReject(split_val):
    global c
    global power
    while True:
        x = random.uniform(0, 1)
        y = random.uniform(0, 1)
        if y*c <= math.pow(x - split_val, power):
            return x

power = 4
t = 0.4  
sum_ = (math.pow(1-t, power + 1) - math.pow(-t, power + 1)) / (power + 1)  #求积分
x = np.linspace(0, 1, 100)
#常数值c
c = 0.6**4/sum_
cc = [c for xi in x]
plt.plot(x, cc, '--',label='c*f(x)')
#目标概率密度函数的值f(x)
y = [math.pow(xi - t, power)/sum_ for xi in x]
plt.plot(x, y,label='f(x)')
#采样10000个点
samples = []
for  i in range(10000):
    samples.append(AceeptReject(t))
plt.hist(samples, bins=50, normed=True,label='sampling')
plt.legend()
plt.show()

6.重要性采样

要性采样是蒙特卡洛方法中的一个重要策略。该方法不改变统计量，只改变概率分布，可以用来降低方差。

7.马尔科夫链

马尔科夫链的细致平稳条件

8、MCMC

8.1 M-H采样

import random
import math
from scipy.stats import norm
import matplotlib.pyplot as plt
%matplotlib inline

def norm_dist_prob(theta):
    y = norm.pdf(theta, loc=3, scale=2)
    return y

T = 5000
pi = [0 for i in range(T)]
sigma = 1
t = 0
while t < T-1:
    t = t + 1
    pi_star = norm.rvs(loc=pi[t - 1], scale=sigma, size=1, random_state=None)
    alpha = min(1, (norm_dist_prob(pi_star[0]) / norm_dist_prob(pi[t - 1])))

    u = random.uniform(0, 1)
    if u < alpha:
        pi[t] = pi_star[0]
    else:
        pi[t] = pi[t - 1]


plt.scatter(pi, norm.pdf(pi, loc=3, scale=2))
num_bins = 50
plt.hist(pi, num_bins, normed=1, facecolor='red', alpha=0.7)
plt.show()

M-H采样总结

8.2 Gibbs采样

二维Gibbs采样


Gibbs抽样是MCMC的一个特例，它交替的固定某一维度，然后通过其他维度的值来抽样该维度的值，注意，gibbs采样只对z是高维（2维以上）（Gibbs sampling is applicable in situations where Z has at least two dimensions）情况有效。






二维Gibbs采样实例

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 xrange(N):
    for j in xrange(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()

多维Gibbs采样

1、

数学八卦