蒙特卡洛方法与 MCMC 采样

一、蒙特卡洛方法

  1. 蒙特卡洛方法Monte Carlo 可以通过采用随机投点法来求解不规则图形的面积。

    求解结果并不是一个精确值,而是一个近似值。当投点的数量越来越大时,该近似值也越接近真实值。

  2. 蒙特卡洛方法也可以用于根据概率分布来随机采样的任务。

1.1 布丰投针问题

  1. 布丰投针问题是1777年法国科学家布丰提出的一种计算圆周率的方法:随机投针法。其步骤为:

    • 首先取一张白纸,在上面绘制许多条间距为 的平行线。

    • 取一根长度为 的针,随机地向纸上投掷 次,观测针与直线相交的次数,记做

    • 计算针与直线相交的概率 。可以证明这个概率 。因此有:

  2. 由于向纸上投针是完全随机的,因此用二维随机变量 来确定针在纸上的具体位置。其中:

    • 表示针的中点到平行线的距离,它是 之间的均匀分布。
    • 表示针与平行线的夹角,它是 之间的均匀分布。

    时,针与直线相交。

    由于 相互独立,因此有概率密度函数:

    因此,针与直线相交的概率为:

    根据 即可得证。

  3. 布丰投针问题中,蒙特卡洛方法是利用随机投点法来求解面积 。因为曲线的积分就是面积,这里的曲线就是概率密度函数

1.2 蒙特卡洛积分

  1. 对于函数 ,其在区间 上的积分 可以采用两种方法来求解:投点法、期望法。

  2. 投点法求积分:对函数 ,对其求积分等价于求它的曲线下方的面积。

    此时定义一个常数 ,使得 ,该常数在区间 上的面积就是矩形面积

    随机向矩形框中随机的、均匀的投点,设落在函数 下方的点为绿色,落在 之间的点为红色。则有:落在 下方的点的概率等于 的面积比上矩形框的面积

    具体做法是:从 之间的均匀分布中采样 ,从 之见的均匀分布中采样 构成一个随机点。

    • ,则说明该随机点在函数 下方,染成绿色。
    • ,则说明该随机点在函数 上方,染成红色。

    假设绿色点有 个,红色点有 个,总的点数为 ,因此有:

  3. 期望法求积分:假设需要求解积分 ,则任意选择一个概率密度函数 ,其中 满足条件:

    令:

    则有: ,它刚好是一个期望:设随机变量 服从分布 ,则

    则期望法求积分的步骤是:

    • 任选一个满足条件的概率分布
    • 根据 ,生成一组服从分布 的随机数
    • 计算均值 ,并将 作为 的近似。
  4. 在期望法求积分中,如果 均为有限值,则 可以取均匀分布的概率密度函数:

    此时

    其物理意义为: 为在区间 上函数的平均高度,乘以区间宽度 就是平均面积。

  5. 对于期望 ,如果 或者 的表达式比较复杂,则也可以转化为另一个期望的计算。

    选择一个比较简单的概率密度函数 ,根据:

    ,则原始期望转换为求另一个期望 。此时可以使用期望法求积分的策略计算。

1.3 蒙特卡洛采样

  1. 采样问题的主要任务是:根据概率分布 ,生成一组服从分布 的随机数

    • 如果 就是均匀分布,则均匀分布的采样非常简单。

    • 如果 是非均匀分布,则可以通过均匀分布的采样来实现。其步骤是:

      • 首先根据均匀分布 随机生成一个样本

      • 为概率分布 的累计分布函数:

        ,计算得到 ,其中 为反函数,则 为对 的采样。

  2. 通过均匀分布的采样的原理:假设随机变量 满足 ,则 的概率分布为:

    因为 上面的均匀分布,因此 为概率分布 的累计分布函数,因此 。因此上式刚好等于 ,即: 服从概率分布

    这其中有两个关键计算:

    • 根据 ,计算累计分布函数
    • 根据 计算反函数

    如果累计分布函数无法计算,或者反函数难以求解,则该方法无法进行。

  3. 对于复杂的概率分布 ,难以通过均匀分布来实现采样。此时可以使用接受-拒绝采样 策略。

    • 首先选定一个容易采样的概率分布 ,选择一个常数 ,使得在定义域的所有位置都满足

    • 然后根据概率分布 随机生成一个样本

    • 计算 ,以概率 接受该样本。

      具体做法是:根据均匀分布 随机生成一个点 。如果 ,则接受该样本;否则拒绝该样本。

      或者换一个做法:根据均匀分布 生成一个随机点,如果该点落在灰色区间()则拒绝该样本;如果该点落在白色区间()则接受该样本。

  4. 接受-拒绝采样 在高维的情况下会出现两个问题:

    • 合适的 分布比较难以找到。
    • 难以确定一个合理的 值。

    这两个问题会导致拒绝率很高,无效计算太多。

二、马尔可夫链

  1. 马尔可夫链是满足马尔可夫性质的随机过程。

    马尔可夫链 描述了一个状态序列,其中每个状态值取决于前一个状态。 为随机变量,称为时刻 的状态,其取值范围称作状态空间。

    马尔可夫链的数学定义为:

2.1 马尔可夫链示例

  1. 社会学家把人按照经济状况分成三类:下层、中层、上层。用状态 1,2,3 代表着三个阶层。社会学家发现:决定一个人的收入阶层的最重要因素就是其父母的收入阶层。

    • 如果一个人的收入属于下层,则他的孩子属于下层的概率是 0.65,属于中层的概率是 0.28,属于上层的概率是 0.07 。
    • 如果一个人的收入属于中层,则他的孩子属于下层的概率是 0.15,属于中层的概率是 0.67,属于上层的概率是 0.18 。
    • 如果一个人的收入属于上层,则他的孩子属于下层的概率是 0.12,属于中层的概率是 0.36,属于上层的概率是 0.52 。

    从父代到子代,收入阶层的变化的转移概率如下:

     子代阶层1子代阶层2子代阶层3
    父代阶层10.650.280.07
    父代阶层20.150.670.18
    父代阶层30.120.360.52
  2. 使用矩阵的表示方式,转移概率矩阵记作:

    假设当前这一代人在下层、中层、上层的人的比例是概率分布 ,则:

    • 他们的子女在下层、中层、上层的人的概率分布是
    • 他们的孙子代的分布比例将是
    • ....
    • 代子孙在下层、中层、上层的人的比例是
  3. 假设初始概率分布为 ,给出前 14 代人的分布状况:

    可以看到从第 9 代开始,阶层分布就趋向于稳定不变。

  4. 如果换一个初始概率分布为 ,给出前 14 代人的分布状况:

    可以发现到第 8 代又收敛了。

  5. 两次不同的初始概率分布,最终都收敛到概率分布 。 这说明收敛的行为和初始概率分布 无关,而是由概率转移矩阵 决定的。

    计算

    可以看到 :

    发现当 足够大的时候, 矩阵 收敛且每一行都稳定收敛到

2.2 平稳分布

  1. 马尔可夫链定理:如果一个非周期马尔可夫链具有转移概率矩阵 ,且它的任何两个状态是联通的,则有:

    其中:

    • 为所有可能的状态。

    • 是转移概率矩阵 的第 行第 列的元素,表示状态 转移到状态 的概率。

    • 概率分布 是方程 的唯一解,其中

      称概率分布 为马尔可夫链的平稳分布。

  2. 注意,在马尔可夫链定理中:

    • 马尔可夫链的状态不要求有限,可以是无穷多个。

    • 非周期性在实际任务中都是满足的。

    • 两个状态的连通指的是:状态 可以通过有限的 步转移到达 (并不要求从状态 可以直接一步转移到状态 )。

      马尔可夫链的任何两个状态是联通的含义是:存在一个 ,使得矩阵 中的任何一个元素的数值都大于零。

  3. 从初始概率分布 出发,在马尔可夫链上做状态转移,记时刻 的状态 服从的概率分布为 , 记作 ,则有:

    假设到达第 步时,马尔可夫链收敛,则有:

    所以 都是同分布的随机变量(当然它们并不独立)。

    如果从一个具体的初始状态 开始,然后沿着马尔可夫链按照概率转移矩阵做调整,则得到一个转移序列

    根据马尔可夫链的收敛行为,当 较大时, 将是平稳分布 的样本。

  4. 定理:如果非周期马尔可夫链的转移矩阵 和某个分布 满足: ,则 是马尔可夫链的平稳分布。

    这被称作马尔可夫链的细致平稳条件detailed balance condition ,其证明如下:

    .

三、MCMC 采样

  1. 概率图模型中最常用的采样技术是马尔可夫链蒙特卡罗方法Markov Chain Monte Carlo:MCMC

    给定连续随机变量 的概率密度函数 ,则 在区间 中的概率可以计算为:

    如果函数 , 则可以计算 的期望:

    • 如果 不是一个单变量,而是一个高维的多元变量 ,且服从一个非常复杂的分布,则对于上式的求积分非常困难。为此,MCMC先构造出服从 分布的独立同分布随机变量 , 再得到 的无偏估计:

    • 如果概率密度函数 也很复杂,则构造服从 分布的独立同分布随机变量也很困难。MCMC 方法就是通过构造平稳分布为 的马尔可夫链来产生样本。

3.1 MCMC 算法

  1. MCMC 算法的基本思想是:先设法构造一条马尔可夫链,使其收敛到平稳分布恰好为 。然后通过这条马尔可夫链来产生符合 分布的样本。最后通过这些样本来进行估计。

    这里马尔可夫链的构造至关重要,不同的构造方法将产生不同的MCMC算法。Metropolis-Hastings:MH算法是MCMC的重要代表。

  2. 假设已经提供了一条马尔可夫链,其转移矩阵为 。目标是另一个马尔科夫链,使转移矩阵为 、平稳分布是

    通常 ,即 并不满足细致平稳条件不成立。但是可以改造已有的马尔可夫链,使得细致平稳条件成立。

    引入一个函数 ,使其满足: 。如:取 ,则有:

    令: ,则有 。其中 构成了转移矩阵 。而 恰好满足细致平稳条件,因此它对应的马尔可夫链的平稳分布就是

  3. 在改造 的过程中引入的 称作接受率。其物理意义为:在原来的马尔可夫链上,从状态 的概率跳转到状态 的时候,以 的概率接受这个转移。

    • 如果接受率 太小,则改造马尔可夫链过程中非常容易原地踏步,拒绝大量的跳转。这样使得马尔可夫链遍历所有的状态空间需要花费太长的时间,收敛到平稳分布 的速度太慢。

    • 根据推导 ,如果将系数从 1 提高到 ,则有:

      于是: