蒙特卡洛方法
蒙特卡罗原来是一个赌场的名称，用它作为名字大概是因为蒙特卡罗方法是一种随机模拟的方法，这很像赌博场里面的扔骰子的过程。最早的蒙特卡罗方法都是为了求解一些不太好求解的求和或者积分问题。比如积分：
如果我们很难求解出f(x)的原函数，那么这个积分比较难求解。当然我们可以通过蒙特卡罗方法来模拟求解近似值。如何模拟呢？假设我们函数图像如下图：
则一个简单的近似求解方法是在
之间随机的采样一个点。比如
,然后用
代表在区间上所有的
的值。那么上面的定积分的近似求解为：
当然，用一个值代表区间上所有的的值，这个假设太粗糙。那么我们可以采样区间的
个值：
，用它们的均值来代表区间上所有的的值。这样我们上面的定积分的近似求解为：
虽然上面的方法可以一定程度上求解出近似的解，但是它隐含了一个假定，即x在[a,b]之间是均匀分布的，而绝大部分情况，
在之间不是均匀分布的。如果我们用上面的方法，则模拟求出的结果很可能和真实值相差甚远。
怎么解决这个问题呢？ 如果我们可以得到在的概率分布函数
，那么我们的定积分求和可以这样进行：
上式最右边的这个形式就是蒙特卡罗方法的一般形式。当然这里是连续函数形式的蒙特卡罗方法，但是在离散时一样成立。
可以看出，最上面我们假设在之间是均匀分布的时候，
，带入我们有概率分布的蒙特卡罗积分的上式，可以得到：
也就是说，我们最上面的均匀分布也可以作为一般概率分布函数在均匀分布时候的特例。那么我们现在的问题转到了如何求出的分布的若干和样本上来。
马尔科夫链
马尔科夫链定义本身比较简单，它假设某一时刻状态转移的概率只依赖于它的前一个状态。
如果用精确的数学定义来描述，则假设我们的序列状态是
，那么我们在时刻
的状态的条件概率仅仅依赖于时刻
，即：
既然某一时刻状态转移的概率只依赖于它的前一个状态，那么我们只要能求出系统中任意两个状态之间的转换概率，这个马尔科夫链的模型就定了。
这个马尔科夫链是表示股市模型的，共有三种状态：牛市（Bull market）, 熊市（Bear market）和横盘（Stagnant market）。每一个状态都以一定的概率转化到下一个状态。比如，牛市以0.025的概率转化到横盘的状态。这个状态概率转化图可以以矩阵的形式表示。如果我们定义矩阵阵
某一位置
的值为
,即从状态
转化到状态
的概率，并定义牛市为状态0， 熊市为状态1, 横盘为状态2. 这样我们得到了马尔科夫链模型的状态转移矩阵为：
马尔科夫链模型状态转移矩阵的性质
transfer_matrix = np.array([[0.6,0.2,0.2],[0.3,0.4,0.3],[0,0.3,0.7]],dtype='float32')start_matrix = np.array([[0.5,0.3,0.2]],dtype='float32')value1 = []value2 = []value3 = []for i in range(30):start_matrix = np.dot(start_matrix,transfer_matrix)value1.append(start_matrix[0][0])value2.append(start_matrix[0][1])value3.append(start_matrix[0][2])print(start_matrix)#进行可视化x = np.arange(30)plt.plot(x,value1,label='cheerful')plt.plot(x,value2,label='so-so')plt.plot(x,value3,label='sad')plt.legend()plt.show()
可以发现，从10轮左右开始，我们的状态概率分布就不变了，一直保持在[0.23076934,0.30769244,0.4615386]。
如果一个非周期的马尔科夫链有状态转移矩阵, 并且它的任何两个状态是连通的，那么
与无关，我们有：
是方程
的唯一非负解，其中
上面的性质中需要解释的有：
- 非周期的马尔科夫链：这个主要是指马尔科夫链的状态转化不是循环的，如果是循环的则永远不会收敛。幸运的是我们遇到的马尔科夫链一般都是非周期性的。用数学方式表述则是：对于任意某一状态，
为集合
的最大公约数，如果
，则该状态为非周期的。 - 任何两个状态是连通的：这个指的是从任意一个状态可以通过有限步到达其他的任意一个状态，不会出现条件概率一直为0导致不可达的情况，即
中任意一个元素都大于零。 - 马尔科夫链的状态数可以是有限的，也可以是无限的。因此可以用于连续概率分布和离散概率分布。
- 我们用
表示在马氏链上跳转第 
步后所处的状态，如果 
存在，很容易证明以上定理的第三个结论。由于
两边取极限就可以得到：
由马氏链收敛的定理, 概率分布
将收敛到平稳分布 
。假设到第 
步的时候马氏链收敛，则有 
，所以 
都是同分布的随机变量，当然他们并不独立。如果我们从一个具体的初始状态 
开始, 沿着马氏链按照概率转移矩阵做跳转，那么我们得到一个转移序列 
，由于马氏链的收敛行为， 
都将是平稳分布 
的样本。
采样
概率分布采样
可以通过PDF求得CDF，再通过CDF的反函数来采样。例如
分布
概率密度函数（PDF）：
def gumbel_pdf(x, mu=0, beta=1):z = (x - mu) / betareturn np.exp(-z - np.exp(-z)) / beta
累计密度函数（CDF）：
def gumbel_cdf(x, mu=0, beta=1):z = (x - mu) / betareturn np.exp(-np.exp(-z))
CDF的反函数：
def inv_gumbel_cdf(y, mu=0, beta=1, eps=1e-20):return mu - beta * np.log(-np.log(y + eps))
我们可以利用反函数法和生成服从 Gumbel 分布的随机数:
def sample_gumbel(shape):p = np.random.random(shape)return inv_gumbel_cdf(p)
接受-拒绝采样
对于概率分布不是常见的分布，一个可行的办法是采用接受-拒绝采样来得到该分布的样本。既然
太复杂在程序中没法直接采样，那么我设定一个程序可采样的分布
比如高斯分布，然后按照一定的方法拒绝某些样本，以达到接近分布的目的，其中叫做 proposal distribution。
具体采用过程：设定一个方便采样的常用概率分布函数，以及一个常量
，使得总在
的下方。如上图。
首先，采样得到的一个样本
，采样方法如概率分布采样。然后，从均匀分布
中采样得到一个值
。如果落在了上图中的灰色区域，则拒绝这次抽样，否则接受这个样本。重复以上过程得到个接受的样本
，则最后的蒙特卡罗方法求解结果为：
整个过程中，我们通过一系列的接受拒绝决策来达到用模拟概率分布的目的。
代码实现
假如我们的目标概率密度函数是
，对此分布生成样本。
准备工作：简单起见，proposal distribution选择uniform(0,1)分布，概率密度函数为
- 计算常数值
```python
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(): global c global power global t while True: x = random.uniform(0, 1) y = random.uniform(0, 1) if y*c <= math.pow(x - t, 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=’cf(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()) plt.hist(samples, bins=50, normed=True,label=’sampling’) plt.legend() plt.show()
<a name="Dfqfg"></a>### 基于马尔科夫链采样> [马尔可夫链](https://www.yuque.com/scistack/recommender/byg2xh#15h5l)如果我们得到了某个平稳分布所对应的马尔科夫链状态转移矩阵，我们就很容易采样出这个平稳分布的样本集。<br />现在我们可以开始采样了，首先，基于初始任意简单概率分布比如高斯分布采样得到状态值，基于条件概率分布采样状态值，一直进行下去，当状态转移进行到一定的次数时，比如到次时，我们认为此时的采样集即是符合我们的平稳分布的对应样本集，可以用来做蒙特卡罗模拟求和了。<a name="mcUJx"></a>##### 马尔科夫链的细致平稳条件(Detailed Balance Condition)给定一个概率平稳分布, 很难直接找到对应的马尔科夫链状态转移矩阵。而只要解决这个问题，我们就可以找到一种通用的概率分布采样方法，进而用于蒙特卡罗模拟。<br />在解决从平稳分布  , 找到对应的马尔科夫链状态转移矩阵  之前，我们还需要先看看马尔科夫链的细致平稳条件。定义如下：<br />如果非周期马尔科夫链状态转移矩阵  和概率分布  对所与的  满足：<br /><br />则称概率分布  是状态转移矩阵  的平稳分布(Stationary Distribution)。**这是个充分非必要条件**（分布为二维时是充要条件）。证明如下：<br />因为 ，所以<br />即满足马尔可夫链的收敛性质。也就是说，只要我们找到了可以使概率分布满足细致平稳分布的矩阵即可。这给了我们寻找从平稳分布, 找到对应的马尔科夫链状态转移矩阵的新思路。<br />不过不幸的是，仅仅从细致平稳条件还是很难找到合适的矩阵。比如我们的目标平稳分布是,随机找一个马尔科夫链状态转移矩阵,它是很难满足细致平稳条件的，即：<a name="d4qiz"></a>#### MCMC由于一般情况下，目标平稳分布和某一个马尔科夫链状态转移矩阵不满足细致平稳条件，即，我们引入一个使得上式取等号：<br /><br />怎样才能成立呢，其实很简单，按照对称性：<br /><br />然后我们就可以得到了分布  对让马尔科夫链状态转移矩阵 <br />其中  一般称之为接受率，取值在  之间，可以理解为一个概率值。这很像[接受-拒绝采样](https://www.yuque.com/scistack/recommender/byg2xh#2iRL2)，那里是以一个常用分布通过一定的接受-拒绝概率得到一个非常见分布， 这里是以一个常见的马尔科夫链状态转移矩阵  通过一定的接受-拒绝概率得到目标转移矩阵  ,两者的解决问题思路是类似的。<br /><br />MCMC采样算法如下：1. 输入任意给定的马尔科夫链状态转移矩阵 ，目标平稳分布  ，设定状态转移次数阈值 ，需要的样本数  ;1. 从任意简单概率分布得到初始状态值 ；1. 1. 从条件概率分布  得到样本值 1. 从均匀分布中采样1. 如果, 则接受转移，即1. 否则不接受转移，样本集即为我们需要的平稳分布对应的样本集。<br />上述过程中  说的都是离散的情形，事实上即便这两个分布是连续的，以上算法仍然是有效，于是就得到更一般的连续概率分布  的采样算法，而  就是任意一个连续二元概率分布对应的条件分布。<br />但是这个采样算法还是比较难在实际中应用，因为在第三步中，由于  可能非常的小，比如0.1，导致我们大部分的采样值都被拒绝转移，采样效率很低。有可能我们采样了上百万次马尔可夫链还没有收敛，也就是上面这个  要非常非常的大，这让人难以接受，怎么办呢？这时就轮到我们的M-H采样出场了。<a name="gS2Bh"></a>#### M-H采样M-H采样解决了我们上一节[MCMC采样](#d4qiz)接受率过低的问题。<br />我们可以对  两边进行扩大，此时细致平稳条件也是满足的，我们将等式扩大  倍，使  （精确来说是使得两边最大的扩大为1），这样我们就提高了采样中的跳转接受率，所以我们可以取 <br />M-H采样算法如下：1. 输入任意给定的马尔科夫链状态转移矩阵 ，目标平稳分布  ，设定状态转移次数阈值 ，需要的样本数  ;1. 从任意简单概率分布得到初始状态值 ；1. 1. 从条件概率分布  得到样本值 1. 从均匀分布中采样1. 如果, 则接受转移，即1. 否则不接受转移，很多时候我们马尔科夫转移状态矩阵是对称的，因此接受率可以写成： <a name="3adbi"></a>##### <br /><a name="LBBfK"></a>##### 代码实现假设目标平稳分布是一个均值3，标准差2的正态分布，而选择的马尔可夫链状态转移矩阵  的条件转移概率是以  为均值,方差1的正态分布在位置  的值。```pythonimport randomfrom scipy.stats import normimport matplotlib.pyplot as pltimport seaborn as sns%matplotlib inlinesns.set_style('darkgrid')plt.rcParams['figure.figsize'] = (12, 8)def norm_dist_prob(theta):y = norm.pdf(theta, loc=3, scale=2)return ydef norm_Q_prob(x_prev):return norm.rvs(loc=x_prev, scale=1, size=1, random_state=None)T = 5000pi = [0 for i in range(T)]t = 0while t < T-1:t = t + 1pi_star = norm_Q_prob(pi[t-1]) #状态转移进行随机抽样alpha = min(1, (norm_dist_prob(pi_star[0]) / norm_dist_prob(pi[t - 1]))) #alpha值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),label='Target Distribution')num_bins = 100plt.hist(pi, num_bins, normed=1, facecolor='red', alpha=0.7,label='Samples Distribution')plt.legend()plt.show()
M-H的问题：
- 我们的数据特征非常的多，M-H采样由于接受率计算式
的存在，在高维时需要的计算时间非常的可观，算法效率很低。同时 
一般小于1，有时候辛苦计算出来却被拒绝了。能不能做到不拒绝转移呢？ - 由于特征维度大，很多时候我们甚至很难求出目标的各特征维度联合分布，但是可以方便求出各个特征之间的条件概率分布。这时候我们能不能只有各维度之间条件概率分布的情况下方便的采样呢？
Gibbs采样
重新寻找合适的细致平稳条件
在M-H采样中我们通过引入接受率使细致平稳条件满足。现在我们换一个思路。
从二维的数据分布开始，假设 
是一个二维联合数据分布，观察第一个特征维度相同的两个点 
和 
，容易发现下面两式成立：
由于两式的右边相等，因此我们有：
也就是：
观察上式再观察细致平稳条件的公式，我们发现在 
这条直线上，如果用条件概率分布 
作为马尔科夫链的状态转移概率，则任意两个点之间的转移满足细致平稳条件！同样的道理，在 
这条直线上，如果用条件概率分布 
作为马尔科夫链的状态转移概率，则任意两个点之间的转移也满足细致平稳条件。那是因为假如有一点 
,我们可以得到：
基于上面的发现，我们可以这样构造分布 
的马尔可夫链对应的状态转移矩阵 
:
有了上面这个状态转移矩阵，我们很容易验证平面上的任意两点 
，满足细致平稳条件：
二维Gibbs采样
利用上一节找到的状态转移矩阵，我们就得到了二维Gibbs采样，这个采样需要两个维度之间的条件概率。算法如下：
- 输入平稳分布
，设定状态转移次数阈值 
，需要的样本个数 
; - 随机初始化初始状态值
和 
; - for
:- 从条件概率分布
中采样得到样本 
- 从条件概率分布
中采样得到样本 
- 从条件概率分布
样本集： 
即为我们需要的平稳分布对应的样本集。
整个采样过程中，马氏链的转移只是轮换的沿着坐标轴 ，采样的过程为：
马氏链收敛后，最终得到的样本就是 
的样本，而收敛之前的阶段称为 burn-in period。额外说明一下，我们看到教科书上的 Gibbs Sampling 算法大都是坐标轴轮换采样的，但是这其实是不强制要求的。最一般的情形可以是，在 
时刻，可以在 
轴和 
轴之间随机的选一个坐标轴，然后按条件概率做转移，马氏链也是一样收敛的。轮换两个坐标轴只是一种方便的形式。
多维Gibbs采样
上面的这个算法推广到多维的时候也是成立的。比如一个 
维的概率分布 
，可以通过在 
个坐标轴上轮换采样，来得到新的样本。对于轮换到的任意一个坐标轴 
上的转移，马尔科夫链的状态转移概率为 
，即固定 
个坐标轴，在某一个坐标轴上移动，即每次转移只允许某一个维度上变化，可以解决维数灾难问题。算法如下：
- 输入平稳分布
，设定状态转移次数阈值 
，需要的样本个数 
- 随机初始化初始状态值
- for
:- 从条件概率分布
中采样得到样本 
- 从条件概率分布
中采样得到样本 
- ….
- 从条件概率分布
中采样得到样本 
- ….
- 从条件概率分布
中采样得到样本 
- 从条件概率分布
样本集： 
即为我们需要的平稳分布对应的样本集。当然这些样本并不独立，但是我们此处要求的是采样得到的样本符合给定的概率分布，并不要求独立。
整个采样过程和Lasso回归的坐标轴下降法算法非常类似，只不过Lasso回归是固定 
个特征，对某一个特征求极值。而Gibbs采样是固定 
个特征在某一个特征采样。
同样的，轮换坐标轴不是必须的，我们可以随机选择某一个坐标轴进行状态转移，只不过常用的Gibbs采样的实现都是基于坐标轴轮换的。
**
二维Gibbs采样实例python实现
假设我们要采样的是一个二维正态分布 
，其中：
， 
;
而采样过程中的需要的状态转移条件分布为：
import math
from mpl_toolkits.mplot3d import Axes3D
from scipy.stats import multivariate_normal
import random
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
sns.set_style('darkgrid')
plt.rcParams['figure.figsize'] = (12, 8)
m1 = 5
m2 = -1
s1 = 1
s2 = 2
rho = 0.5
samplesource = multivariate_normal(mean=[m1, m2], cov=[[s1, s1*s2*rho],[s1*s2*rho, s2*s2]])
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))
N = 5000
K = 20
x_res = []
y_res = []
z_res = []
y = m2
for i in range(N):
for j in range(K):
x = p_xgiveny(y, m1, m2, s1, s2) #y给定得到x的采样
y = p_ygivenx(x, m1, m2, s1, s2) #x给定得到y的采样
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,label='x')
plt.hist(y_res, num_bins, normed=1, facecolor='red', alpha=0.5,label='y')
plt.title('Histogram')
plt.legend()
plt.show()
若有收获，就点个赞吧
0 人点赞
