30分钟了解蒙特卡洛方法

30分钟入门蒙特卡洛模拟

本文主要讲解三部分：

• 背景介绍
• 蒙特卡洛方法介绍
• 结果展示

---

背景介绍

  这一小节我们简要介绍一下引出蒙特卡洛方法的实际场景。
  机器学习/深度学习中的图像叠加文字识别需要大量的训练样本，自动生成样本（使用程序在背景图片上叠加文字）是一种样本的获取方式。但色彩值（为了兼顾各方向的同学，原谅我用一个这么不专业的词汇，此值可以是RGB到[0,1]区间的映射，让它能代表颜色的性质）的选择很重要，为了防止（控制）发生叠加文字与背景图片的色彩值相近的情况发生，叠加文字的色彩值最好服从我们指定的概率分布。这样就需要根据指定的概率分布来产生色彩值——蒙特卡洛方法擅长解决的问题。

蒙特卡洛方法介绍

  蒙特卡洛方法的应用场景很多，横跨物理、金融、计算机。拿计算机科学来举例，自然语言处理中的LDA模型，hinton较早提出的深度学习模型DBN都用到了蒙特卡洛方法。此文第一部分简要介绍了实际问题，简而言之蒙特卡洛方法就是生成样本，即蒙特卡洛采样。即根据某已知分布的概率密度函数f(x)，产生服从此分布的样本X。
  下面首先介绍一种最简单最易理解的蒙特卡洛方法——Accept-Rejection method（下文称接受拒绝采样），然后给出这个方法的直观解释，最后证明方法的正确性。

1. 接受拒绝采样介绍

   • 需求：
     根某已知分布的概率密度函数f(x)，产生服从此分布的样本X

   • 准备工作：

     1. 需要一个辅助的“建议分布”G（已知其概率密度函数g(y)）来产生候选样本。

        由于需要用此分布来产生候选样本，因此需要我们能够从G抽样。即我们必须能够生成服从此概率分布的样本Y。比如可以选择均匀分布、正态分布，大部分语言都有生成其样本的实现。

     2. 还需要另一个辅助的均匀分布U(0,1)。

        当然也要求我们能从U(0,1)抽样。

     3. 计算一个常数值c。

        c为满足c∗g(x)≥f(x)的最小值

   • 开始样本生成：

     1. 从建议分布G抽样，得到样本Y。
     2. 从分布U(0,1)抽样，得到样本U。
     3. 如果 U≤f(Y)c∗g(Y)，则令X=Y（接受Y），否则继续执行步骤1（拒绝）。

2. 接受拒绝采样的直观解释

     为了方便理解，我们使用U(0,1)来作为“建议分布”G，这样g(x)=1∀x∈[0,1]。如下图所示，我们每次生成的两个样本Y与U，对应下图中矩形内的一点P(Y,U∗c∗g(Y))。接受条件U≤f(Y)c∗g(Y)，即U∗c∗g(Y)≤f(Y)的几何意义是点P在f(x)下方，不接受Y的几何意义是点P在f(x)的上方。在f(x)下方的点(o形状)满足接受条件，上方的点(+形状)不满足接受条件。
   

3. 接受拒绝采样方法有效性证明

     （注：本文不着重公式推导，所以简写了一种方法给大家提供一个思路，有不明白的可以留言。需要着重说明的是，对于以下证明过程，看不懂也没关系，因为不会妨碍你直观形象理解采样方法。对于工程人员来讲能直观形象理解公式比证明更重要。）
     实际上是证明在U≤f(Y)c∗g(Y)的条件下，Y的概率分布服从F(y)，即证明P(Y≤y|U≤f(Y)c∗g(Y))=F(y).证明过程如下
   P(Y≤y|U≤f(Y)c∗g(Y))
   =P(U≤f(Y)c∗g(Y)|Y≤y)∗P(Y≤y)P(U≤f(Y)c∗g(Y))

   =P(U≤f(Y)c∗g(Y)|Y≤y)∗P(Y≤y)1c

   =c∗P(U≤f(Y)c∗g(Y),Y≤y)

   =c∗∫y−∞∫f(Y)c∗g(Y)0g(Y)dUdY

   =∫y−∞f(Y)dY

   =F(y)

   证毕。

应用结果展示

  （假设目标概率密度函数为f(x)=(x−0.4)4∫10(x−0.4)4dx，据此分布生成样本。（为什么用这个目标函数？因为它气质好，不是，是效果好，并且非常简单）

• 准备工作：
  
  1. 建议分布函数G选择U(0,1)，概率密度函数为g(y)=1∀y∈[0,1]。（因为这个分布最简单，求常数值c方便）
  2. 辅助的均匀分布U(0,1)
  3. 计算常数值，首先c=max{0.44/sum,0.64/sum}=0.64/sum (其中sum=∫10(x−0.4)4dx，后面可以化简掉)
• 生成代码如下：

import randomimport mathimport matplotlib.pyplot as pltclass RndGenerator:   """docstring for RndGenerator"""   @staticmethod   def AceeptReject(split_val):   power = 4   scale = max(math.pow(split_val, power), math.pow(1 - split_val, power))   while True:      x = random.uniform(0, 1)      y = random.uniform(0, 1)      if y*scale <= math.pow(x - split_val, power):         return x   if __name__ == "__main__":      count = 101             power2 = 4      t = 0.4      sum_ = (math.pow(1-t, power2 + 1) - math.pow(-t, power2 + 1)) / (power2 + 1)      x = [xi*1.0/(count-1) for xi in range(count)]      #常数值      c = max(math.pow(-t, power2), math.pow(1-t, power2))/sum_      cc = [c for xi in x]      p1, = plt.plot(x, cc, '--')      #目标概率密度函数的值f(x)      xx = [math.pow(xi - t, power2)/sum_ for xi in x]      p2, = plt.plot(x, xx)      plt.legend([p1, p2], ["g(x)=c*f(x)", "f(x)"], loc=6)      #采样10000个点      samples = []      for  i in range(10000):         samples.append(RndGenerator.AceeptReject(t))      p3 = plt.hist(samples, bins=40, normed=True)      plt.show()

• 效果说明：
  下图所示，生成的10000个样本的直方图与目标概率密度函数f(x)吻合。
  

留坑

  关于蒙特卡洛采样的效率问题，大家可以思考下接受一次采样需要循环多少次？等式P(U≤f(Y)c∗g(Y))=1c的意义是什么？这个等式与采样效率有什么关系？建议分布G与采样效率有什么关系?

总结

  这篇文章是蒙特卡洛方法入门介绍，接受-拒绝采样是蒙特卡洛方法最简单的一种，也是所有方法的基础。蒙特卡洛方法的精髓就是随机采样，拟合分布。有了样本点之后就可以做很多事情了，最重要的就是近似复杂积分。此方法在深度学习领域内应用甚广，值得大家深入学习。
  最后年假将至，预祝大家新年快乐。