蒙特卡洛采样与Gibbs采样
蒙特卡洛采样与Gibbs采样
采样
采样问题指的是给定一个特定的概率分布p(z),得到一批符合这个概率分布的样本点。
采样的方法有很多,MCMC是其中的一类方法,意思是利用Mento Carlo和Markov Chain完成采样。
当然,要完成对各种分布的采样,有一个默认的假设,就是我们已经能够对均匀分布进行采样了(后面就专指范围为0-1的均匀分布),也就是编程中通常会用到的伪随机数发生器,在各大编程语言中通常以random命名的模块/方法出现。
蒙特卡洛思想
蒙特卡罗方法是一种通过在一定范围内均匀随机抽样来得到某个结果的计算方法。方法的大致思路框架是:
示例
val n = 100000 * slices // slice是控制抽样量的参数
val count = sc.parallelize(1 to n, slices).map { i =>
val x = random * 2 - 1 // 长度坐标范围-1~1
val y = random * 2 - 1 // 宽度坐标范围-1~1
if (x*x + y*y < 1) 1 else 0 //区分出落在圆圈中的点。圆圈半径为1,所以点到圆心距离小于1
}.reduce(_ + _) // 统计点的数量
println("Pi is roughly " + 4.0 * count / n)
马尔科夫链
[2000 14000]∗[0.970.050.030.95]=[2640 13360][2000\ 14000] * \left [
0.970.05amp;0.03amp;0.950.97amp;0.030.05amp;0.95
\right ]=[2640\ 13360][2000 14000]∗[0.970.050.030.95]=[2640 13360]
P(Xn∣Xn−1,…,X0)=P(Xn∣Xn−1)P(Xn|Xn−1,…,X0)=P(Xn|Xn−1)P(Xn∣Xn−1,…,X0)=P(Xn∣Xn−1)
而迁移矩阵的每一个元素其实就对应着一个条件概率值。所以后面会同时用P(j|i)和PijP_{ij}Pij这两种写法来代表迁移矩阵的某一项。
- 马尔科夫链收敛于平稳分布π,使得对于迁移矩阵P,π是方程πP=π的唯一解。把这个解用向量方式表示,就是该分布在各个取值上对应的概率,即π=[π(0),π(1),…π(i),…],∑iπ(i)=1[π(0),π(1),…π(i),…],∑iπ(i)=1[π(0),π(1),…π(i),…],∑iπ(i)=1
- 迁移矩阵P经过反复迭代之后本身也会收敛(每列结果相同),才能满足平稳分布的要求,即:
MCMC
π(i)Pij≠π(j)Pjiπ(i)P_{ij}≠π(j)P_{ji}π(i)Pij̸=π(j)Pji
Metropolis Hastings算法
# -*- coding: utf-8 -*-
import random
import numpy as np
import pylab as pl
import scipy.special as ss
# 完整的beta分布概率密度函数
def beta(a, b, i):
e1 = ss.gamma(a + b)
e2 = ss.gamma(a)
e3 = ss.gamma(b)
e4 = i ** (a - 1)
e5 = (1 - i) ** (b - 1)
return (e1/(e2*e3)) * e4 * e5
# beta分布概率密度函数去掉前面的常数项之后的形式
def beta_s(a,b,i):
return i**(a-1)*(1-i)**(b-1)
# mcmc模拟
def beta_mcmc(N_hops,a,b):
states = []
cur = random.uniform(0,1) # 初始化状态
for i in range(0,N_hops):
states.append(cur)
next = random.uniform(0,1) #从原来的迁移矩阵P采样,这里假设P是一个基于均匀分布的迁移矩阵
#计算接受率,因为beta分布的常数项会抵消,所以用不带常数项的形式,能大幅提速。而且P是均匀分布,所以也相互抵消了
ap = min(beta_s(a,b,next)/beta_s(a,b,cur),1)
if random.uniform(0,1) < ap: #随机采样决定是否跳转
cur = next
return states[-10000:] #取最后的一部分状态,保证已经收敛
#可视化
def plot_beta(a, b):
Ly = []
Lx = []
i_list = np.mgrid[0:1:100j]
for i in i_list:
Lx.append(i)
Ly.append(beta(a, b, i))
pl.plot(Lx, Ly, label="Real Distribution: a="+str(a)+", b="+str(b))
pl.hist(beta_mcmc(1000000,a,b),normed=True,bins=50, histtype='step',
label="Simulated_MCMC: a="+str(a)+", b="+str(b))
pl.legend()
pl.show()
pl.rcParams['figure.figsize'] = (8.0, 4.0)
plot_beta(2, 5)
输出的可视化结果为:
Gibbs采样
Gibbs Sampling需要应用在至少二维的数据上,所以下面先以二维的情况为例。
当然在实际应用中,无论是MH算法还是GS算法,一般都要应用在多维数据上。我们把上述二维的情况推广到多维,就是完整的GS算法:
可以看到,多维的情况下,每一次采样还是只能移动一个维度,其余的n−1个维度都作为条件。
相比MH算法,GS算法每次跳转都没有被拒绝的可能,所以一般会说GS算法是MH算法的一个特例,即接受率为1的MH算法。
# -*- coding: utf-8 -*-
import pylab as pl
import numpy as np
import math
# x和y两个维度的均值都是0
sx = 8 # x维度正态分布的标准差
sy = 2 # y维度正态分布的标准差
cor = 0.5 # x和y的相关系数
# x维度的概率密度函数
def pdf_gaussian_x(x):
return (1 / (math.sqrt(2 * math.pi) * sx)) * math.exp(-math.pow(x, 2) / (2 * math.pow(sx, 2)))
# 条件分布 p(x|y)
def pxgiveny(y):
return np.random.normal(y * (sx/sy) * cor, sx * math.sqrt(1 - cor * cor))
# 条件分布 p(y|x)
def pygivenx(x):
return np.random.normal(x * (sy/sx) * cor, sy * math.sqrt(1 - cor * cor))
def gibbs(N_hop):
x_states = []
y_states = []
#状态随机初始化
x = np.random.uniform()
y = np.random.uniform()
for _ in range(N_hop):
x = pxgiveny(y) #根据y采样x
y = pygivenx(x) #根据x采样y
x_states.append(x)
y_states.append(y)
return x_states[-1000:], y_states[-1000:]
def plot_gibbs():
x_sim, _ = gibbs(100000)
x1 = np.arange(-30, 30, 1)
pl.hist(x_sim, normed=True, bins=x1, histtype='step', label="Simulated_Gibbs")
x1 = np.arange(-30, 30, 0.1)
px1 = np.zeros(len(x1))
for i in range(len(x1)):
px1[i] = pdf_gaussian_x(x1[i])
pl.plot(x1, px1, label="Real Distribution")
pl.legend()
pl.show()
plot_gibbs()
作图部分为了兼顾方便和清晰,只基于一个维度做了可视化,结果如下图:
应用场景
基于MCMC的采样方法主要有两类应用场景。第一类是应用于求函数的期望值,或者说是积分问题。另一类是求解最优化问题。
积分问题
把采样用在积分计算上,一个应用场景就是用在贝叶斯估计中。首先根据贝叶斯公式,写出后验概率:
E[θ]=∫θP(θ∣X)dθE[θ]=∫θP(θ|X)dθE[θ]=∫θP(θ∣X)dθ
事实上我们利用这个思路还可以直接进行预测,也就是求P(y|X):
假设Z是一个连续随机变量(离散的情况类似),其概率密度函数为p(z),现在需要计算函数f(z)的期望值,即
最优化
# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
import random
import math
# 目标函数,有两个变量。为了让收敛过程看起来清晰一点,找了一个复杂一点的函数
def f(x):
x1 = x[0]
x2 = x[1]
obj = x1**2 + x2**2 - 0.1*math.cos(6.0 * math.pi * x1) - 0.1*math.cos(6.0 * math.pi * x2)
return obj
x_start = [1, -.5] # 初始点
t_start = 1 #起始温度参数
t_end = 0.02 #结束温度参数
n = 50 #温度变化分几轮做完
frac = (t_end/t_start)**(1.0/(n-1.0)) #每次温度缩小系数
#记录每一次温度变化之后的最终采样结果
x = np.zeros((n+1,2))
x[0] = x_start
# 每轮mcmc采样中新一轮跳转状态
xi = np.zeros(2)
xi = x_start
# 每轮mcmc采样中当前最好状态
xc = np.zeros(2)
xc = x[0]
# 对应函数值的当前最好状态和每轮状态
fc = f(xi)
fs = np.zeros(n+1)
fs[0] = fc
t = t_start
for i in range(n):
for j in range(200): # MCMC跳转200次
# 在上一轮采样结果的附近采样
xi[0] = random.random() - 0.5 + xc[0]
xi[1] = random.random() - 0.5 + xc[1]
# 根据gibbs分布计算接受率
a = math.exp(-(f(xi)-fc)/t)
# 如果采样满足接受率条件,或者本来就取到了更小的值,就接受跳转
if (random.random() < a) or (f(xi) < fc):
xc[0] = xi[0]
xc[1] = xi[1]
fc = f(xc)
# 存储每轮MCMC采样结果
x[i+1][0] = xc[0]
x[i+1][1] = xc[1]
fs[i+1] = fc
# 减小温度参数
t = frac * t
# 可视化
fig = plt.figure()
ax1 = fig.add_subplot(211)
ax1.plot(fs,'r.-')
ax1.legend(['Objective'])
ax2 = fig.add_subplot(212)
ax2.plot(x[:,0],'b.-')
ax2.plot(x[:,1],'g--')
ax2.legend(['x1','x2'])
plt.show()
过可视化输出可以看到,目标函数值虽然在收敛过程中会出现变大的情况,但是最终还是落到最小值,而相应的变量也逐渐收敛:
除了用在通用的最优化方法上面,MCMC也可以应用在具体的问题上,比如在The Markov Chain Monte Carlo Revolution一文中,就提到过一个破译凯撒密码的例子。通过随机生成解密key,根据解密key替换原文后的评估分数决定是否跳转(跳转就是交换任意两个字母的替换规则),最后得到收敛的结果就是正确的解密key。
引用
蒙特卡洛采样与Gibbs采样相关推荐
- MCMC详解2——MCMC采样、M-H采样、Gibbs采样(附代码)
MCMC是一种随机采样方法,用来处理一些复杂运算的近似求解.在HMM.LDA等模型中都有重要应用. 上一篇 MCMC详解1--蒙特卡洛方法 目录 1,马尔科夫链模型 1.1 马尔科夫链转移矩阵的性质 ...
- 马氏链,Metropolis-Hastings采样与Gibbs采样的理解(附有python仿真)
文章目录 马氏链 原理 采样方法 MH采样 原理 代码 Gibbs采样 原理 代码 马氏链 原理 采样方法 所谓的采样方法,主要就是利用了马氏链的性质 πn(x)\pi_n(x)πn(x)为一个离散 ...
- 受限玻尔兹曼机准备知识——MCMC方法和Gibbs采样
先点明几个名词 MCMC方法:马尔可夫链-蒙特卡洛方法 (千万别叫成梅特罗波利斯蒙特卡罗方法了) Metropolis-Hastings采样:梅特罗波利斯-哈斯廷斯采样 Gibbs采样:吉布斯采样 ...
- 深度学习 --- 受限玻尔兹曼机RBM(MCMC和Gibbs采样)
上一节我们详细的讲解了马尔科夫过程和马尔科夫链,提及了转移矩阵P,马尔科夫的平稳性(遍历性),需要好好理解马尔科夫的平稳性,因为本节将根据马尔科夫的平稳性进行讲解,同时也介绍了采样的原理和过程.好,到 ...
- 漫谈MCMC与Gibbs采样(三)—— 有趣的马尔科夫链
Markov Chain 提起马尔科夫链,大家应该都不陌生.我第一次接触这个概念,是在大一的C语言编程课中,当时用马尔科夫链来做文章的随机生成.马尔科夫链的思想非常简单,在数学上可以表述如下: (6) ...
- LDA主题模型系列(二)求解之Gibbs采样
本系列分为三部分: LDA基本概念 LDA求解之Gibbs采样 LDA求解之变分推断EM算法 该篇为第二部分:LDA求解之Gibbs采样 对于Gibbs采样不了解的可以参考这里 本文只包括思路,具体的 ...
- 漫谈MCMC与Gibbs采样(一)—— 采样背后的逻辑
采样背后的逻辑 我们为什么需要采样 我们先来思考一个最基本的问题,我们为什么需要采样?之所以讨论这一点,是因为深刻理解我们要解决的问题,可以帮助我们更好地掌握解决问题的方法.尤其是当我们在复杂的模型中 ...
- 马尔可夫蒙特卡洛(MCMC)-从平稳分布,细致平衡到Metropolis-Hastings和Gibbs采样
马尔可夫链性质 马尔可夫假设: 某一时刻状态转移概率只依赖于前一个状态 马尔科夫链状态转移矩阵的性质:平稳分布 与初始的概率分布无关,马尔可夫链在有限次转移之后总能收敛到一个稳定的概率分布,称为平稳分 ...
- MCMC(三)蒙特卡洛之Gibbs采样
对于给定的概率分布p(x),我们希望能有便捷的方式生成它对应的样本.由于马氏链能收敛到平稳分布,于是一个很的漂亮想法是:如果我们 能构造一个转移矩阵为P的马氏链,使得该马氏链的平稳分布恰好是P(x), ...
最新文章
- Android-SharedPreferences
- 语义分割领域开山之作:Google提出用神经网络搜索实现语义分割
- Java提高篇——Java实现多重继承
- cmd长ping记录日志和时间_四个网络命令ping、arp、tracert、route的详细用法
- 复杂网络研究:让世界变得简单
- gdb学习(二)[第二版]
- Wireshark coloring rules tips
- bbs.FISHC.com//python_文件
- 927. 三等分(每日一难phase2--day26)
- 音视频编解码学习详解h264 ,mpeg4 ,aac 等音视频格式
- springboot JWT Token 自动续期的解决方案
- 【公益译文】网络威胁信息共享指南
- python:实现Lempel-Ziv压缩算法(附完整源码)
- P1138 第k小整数
- C++ STL使用实例
- [__ob__: Observer]
- redis系列七LUR清除算法
- 独立同分布(I.I.D.)是什么?(转载)
- 暑假java培训班,分享面经!
- 生产者和消费者模型介绍