前言

本篇文章简要介绍统计模拟（Simulation），也就是如何生成符合特定分布的随机数据的方法

如何获得随机数

在考虑如何生成符合特定概率分布的随机数据之前，我们先考虑最基本的一个问题——如何获得随机数，我们正式定义“获得随机数”这件事为
::: align-center
生成随机且服从U(0,1)均匀分布的数字序列
:::
均匀分布很容易解释，至于为何在[0,1)之间，我们在下一章中会详细说明

下面介绍两种经典的算法

平方取中法

平方取中法考虑从一个N位数开始，将其平方，取中间的N位数作为下一个随机数，并重复此过程

最开始的这个N位数X_n是人给出的，相同的这个数每次生成的随机数一定是唯一且相同的，它就是我们熟悉的随机种子（Seed），一般选一个大的质数

之后将取到的随机数前加0.便可以将其映射到[0,1)区间

这是最早的随机数生成算法之一，但是周期通常较短，质量不高

线性同余法

线性同余法考虑人为给定一个乘子a，增量c，模数m以及初始随机种子X_n，此后通过递推公式
::: align-center
X_{n+1}=(aX_n+c)\bmod{n}
:::
不断迭代获得一组随机数，之后计算X_{n}/m将其映射到(0,1]区间即可
在参数的选择上，模数m通常选择尽可能大的整数，乘子a应与m互质以使序列具有满周期，c通常选择与m互质的数

线性同余法曾是上世纪中叶左右的主流算法，计算简单，但选择不好的参数会导致序列相关性差、周期短，现代算法在其上进行改进（比如著名的MT19937梅森旋转算法）取得了更优的质量

如何获得符合分布的数据

在获得了[0,1)区间上的随机数后，我们便可以考虑“如何将随机数转换为符合概率分布的数据”这件事了，注意到我们的随机数均满足位于[0,1)区间内，而这正好是概率分布的区间，于是我们便成功建立了“随机数”与“概率”之间的关系

这里笔者还希望给读者建立一个“扔飞镖”的直觉。

我们想象将概率密度函数PDF/PMF画在一张二维坐标系的纸上，比如正态分布画在纸上就是一条钟形曲线，它和横坐标轴组成一个钟形区域。

现在我们考虑用一个简单的游戏解释这背后的原理
我们考虑均匀地向这个钟形区域扎飞镖🎯（因为，我们要保证游戏是“公平”的），读者如果实操的话，一定会发现最后扔出的飞镖很难落在钟形曲线的两侧（因为两侧的区域太“窄”了），而更倾向于落在钟形曲线的中间，换句话说，当我们保证我们的飞镖是均匀落在钟形区域的时候，它的实际落点会受到正态分布PDF的限制，我们记落点的横坐标为x，那么落点横坐标x的分布便一定符合正态分布。

上面的方法就是经典的蒙特卡洛方法，而在实际问题中，我们一般考虑如何精确的获得x，毕竟我们还没到真的去扎飞镖的地步

以此为引子，我们先介绍最直观的逆变换采样法，它就是蒙特卡洛方法族中的一种精确蒙特卡洛算法

逆变换采样（Inverse Transformation）

一个最简单的想法便是我们考虑对于给定的，已经生成好的随机数u\in[0,1]，设目标概率分布的CDF为F(x)，如果我们找到x，使得
::: align-center
F(x)=u
:::
那么这个x便是符合我们要求概率分布的随机变量X的一个样本
上式通常又写为x=F^{-1}(u)，其中F^{-1}(u)是给定概率分布CDF的反函数

举个例子
::: align-center
生成符合U(a,b)分布的一组数据
:::
我们假设提前获得的随机数为u，均匀分布的CDF经求得为F(x)=\frac{x-a}{b-a}，我们考虑解方程
::: align-center
\frac{x-a}{b-a}=u\Rightarrow x=a+(b-a)u
:::
则F^{-1}(u)=a+(b-a)u，将我们的随机数u代入便可以得到符合U(a,b)的随机数了

广义逆变换采样

但是，狭义逆变换采样要求F(x)可逆，而对于一般的概率分布（比如正态分布）而言，其CDF一般不可逆/求逆复杂，于是我们引入广义逆变换采样

广义逆变换采样的公式为：
::: align-center
F^{-1}(u)=\inf\{x|F(x)\ge u\}
:::
也就是找到一个x，使得F(x)是第一个大于等于u的值

举个例子，在”骰子点数“这个离散随机变量统计学问题上，它的CDF是一个阶梯函数

• F(1)=\frac{1}{6}
• F(2)=\frac{2}{6}
• ...
• F(6)=1

假若我们要生成一组符合骰子点数分布的随机数据（也就是随机变量X=1,2,3,4,5,6的一组样本），假若我们获得的随机数是0.3，那么显然F(2)是第一个值大于0.3的数，于是我们就可以获得一个样本x=2

而对于连续型随机变量，以标准正态分布N(0,1)为例，其广义逆累计分布函数通常写为\Phi^{-1}(0,1,\alpha)，其值通常表示“第一个概率大于1-\alpha的数”，也就是我们常说的上\alpha分位数，因此广义逆累计分布函数又称为分位函数

接受-拒绝法

接受-拒绝法是另一种在CDF难以求逆时可稳定生成符合分布的随机数的方法，它的核心思想是
::: align-center
既然无法直接逆变换采样，那我们可以考虑使用已知分布进行采样，然后利用目标采样的概率密度对采样结果进行约束
:::
首先介绍目标分布与提议分布

• 目标分布就是要采样的分布，其PDF/PMF记为f(x)，它的直接采样通常十分困难
• 提议分布的PDF/PMF通常记为g(x)，其具有易采性、覆盖性
  • 易采性即容易采样
  • 覆盖性是指可以通过一个常数M，使得Mg(x)可完全“覆盖”住f(x)

它的执行方式如下：

1. 寻找一个提议分布g(x)，使得对于定义域内的几乎所有x，有f(x)\le Mg(x)（覆盖性），同时人容易采样（易采性），比如均匀分布
2. 从提议分布g(x)中独立地抽取一个随机样本X_0，由于提议分布具有易采性，其采样通常可通过简单的方法（比如逆变换采样）获得，比如均匀分布的X_0可以通过X_0=a+(b-a)u获得
3. 计算接受概率r=\frac{f(X_0)}{Mg(X_0)}，由于f(x)\le Mg(x)，这个值通常总在0到1间
4. 在区间[0,1]上生成一个均匀分布的随机数U，如果U\le r，则接受X_0作为来自f(x)的一个有效样本，否则拒绝这个样本，回到步骤1重新采样

它的原理，依然以“扎飞镖”这个场景演示，我们考虑用一个更大的Mg(x)包裹住f(x)，之后我们考虑在x=X_0处画一条竖直的线，这一条线一定穿过f(x)的区域，最终到达g(X_0)

现在我们考虑在这条线上均匀的扎飞镖🎯（还记得吗？我们要保证“公平”），只要飞镖扎中了f(x)的那块区域，那X_0这个点就是符合f(x)分布的一个采样点）

那么读者可能又会问？我们还得保证X_0这个点取的是“公平”的啊！而这就又回到了我们如何给g(x)采样这个问题上了，因为g(x)拥有易采性，所以我们能够轻松的取到符合“公平”的X_0

从中我们也可以窥探到处理数学问题的“升维后降维”思想，因为如果读者仔细思考，接受-拒绝法实际上是考虑把问题升到二维后再降到一维处理，类似的思想也可以被用于处理“二维曲线在某一点的法线”的这个问题上，只要我们把问题升维到高维后考虑函数的梯度，而后再投影到二维，便会发现这个向量就是曲线在该点的法向量，详见笔者于24年3月写的这篇文章

接受拒绝法的算法效率可以用接受率\frac{1}{M}表示，因为在给定X_0=x时，接受样本的条件是r\le\frac{f(x)}{Mg(x)}，其中r\sim U(0,1)，且与X_0独立，因此
::: align-center
P(\text{接受}|X_0=x)=\frac{f(x)}{Mg(x)}
:::
又已知取到X_0=x的概率就是g(x)，则由全概率积分公式得到接受率为
::: align-center
\int \frac{f(x)}{Mg(x)}g(x)dx=\frac{1}{M}\int f(x)dx=\frac{1}{M}
:::
在f(x)和g(x)固定时，M越小，接受率越高，算法效率越高，实际计算中我们通常尽可能取M为待采样分布比提议分布的上确界以保证最佳效率
::: align-center
M=\sup \frac{f(x)}{g(x)}
:::

采样序列分布&高维联合独立分布

前面介绍的采样方法均只针对独立样本，而对于天气的变化、文本的生成等来说，这些就属于前后状态可能相互关联的序列数据，或者高维上的联合分布，因为原则上，如果我们把序列数据x_1,x_2,\cdots,x_t看作一个高维单点X=(x_1,x_2,\cdots,x_t)，则理论上依然可以使用接受-拒绝法等方法从这个序列的联合概率分布f(X)中采样，但因为维度过高，我们的提议分布g(X)将变得很难设计，而如果我们采用简单的提议分布，在高维中它的接受率往往极低，计算效率不高

但如果我们将上面的描述简化，考虑一个系统在任何时刻的状态只依赖于其前一个时刻的状态，那我们就得到了一个马尔科夫模型，而如果我们再考虑系统的真实状态是隐藏的、不可直接观测的，但每个隐藏的状态又会生成一个可观测的输出时，我们就得到了隐马尔科夫模型

而如果我们将问题简化为马尔科夫模型，我们对这个高维分布的采样便不再需要估计一个全局的提议分布，我们只需要基于当前的状态，考虑如何向下一个状态移动。这就好比我们想知道一片非常大的森林里树木的分布，使用接受拒绝法相当于我们每次提出一个更大的但我们已知分布的森林然后与这片森林比较，很显然这是不现实的。但如果我们考虑让一个旅行者每次沿着森林里的树走，由一棵树走向另一颗树，中途记录下树木的位置，只要这个旅行者走的次数足够多，我们便可以通过研究它记录的树木的分布得到这个森林树木分布的估计

在数学上，我们可以证明马尔科夫链具有收敛性，不论我们的初始概率分布如何，经足够大的次数迭代后其最终的状态转移矩阵一定会收敛于某一个确定的值

而如果我们考虑精心设计一些算法，使得“旅行者”在状态空间里长期游走后访问各个区域的频率分布正好收敛于待采样分布的概率分布，此时“旅行者”经过的每一个采样点，便就是待采样分布的一个实际的采样。这些算法属于马尔科夫蒙特卡洛方法族（Markov Chain Monte Carlo，MCMC），MCMC通过构造马尔科夫链随机游走避免了需要提出全局提议分布的问题，从而能有效的采样高维空间，特别适用于序列分布&高维联合独立分布的采样

如何从已知的离散状态的马尔科夫链中采样下一个状态

我们先从一个最简单的问题开始，假如我们有一个定义好的马尔科夫链（已知其转移概率矩阵为P），并且已知当前时刻t的状态为X_t=i，如何采样下一个状态X_{t+1}？

假如当前的状态i已知，由于马尔科夫链是已知的，我们便知道了它向下一个阶段各个状态转移的概率分布，在直观上，你找到了转移概率矩阵的第i行，这一行的每一列都是一个概率值，整个第i行构成了一个离散的概率分布

那么针对下一阶段的采样就变成了事实上的，可以由广义逆变换采样解决的问题，我们只需要计算出此离散分布的CDF，抽取随机数u\in [0,1]，通过广义逆变换采样便可以得到下一状态的采样i+1