这篇文章的问题来源于博客的随机初始化单位向量，实际上那篇文章也是在用均匀分布去构造另一种分布。在搜集那篇文章的答案的时候，碰到了“如何用均匀分布构造正态分布”的问题，觉得很有代表性，就记录了下来。

Inverse Transform Method - 逆变换方法

先抛开构造具体的正态分布这个议题，我们先来看看用均匀分布去构造其他分布的一种较为通用的方法。当我们知道了要构造的目标分布的密度函数$f_{X}(x)$，如果你可以求出它的累积分布函数（Cumulative Distribution Function-CFD）$F_{X}(x)$，那么我们就可以用0到1之间的一个均匀分布变量$U(0,1)$去构造这个目标分布。先求出累积分布函数的逆函数，再将均匀分布作为变量输入该逆函数，即$F^{-1}_{X}(U)$，输出就遵循要构造的目标分布了。

这个证明也非常简单，要证明两个随机变量是同分布的，最直接的办法就是证明他们的累积分布函数是相同的。假设随机变量$X$是目标分布，那么就要证明随机变量$Y=F_{X}^{-1}(U)$的累积分布函数也是$F_{X}$。变量$Y$的累积分布函数$F_{Y}$可以写为

$$
F_{Y}(x) = P(Y<x) \\
F_{Y}(x) = P(F^{-1}_{X}(U)<x)
$$

对于不等式$F_{X}^{-1}(U)<x$因为累积分布函数（注意不是概率密度函数）一定是单调递增的，所以这个不等式两边代入$F_{X}(x)$，不等号的方向还是一样的，即

$$
F_{X}(F_{X}^{-1}(U)) < F_{X}(x)\\
U < F_X(x)
$$

所以$F_Y(x)$就可以写为

$$
F_Y(x) = P(U<F_X(x))
$$

这里的$U$是0到1上的均匀分布，那么$U$小于$F_{X}(x)$的概率就为$(F_X(x)-0)*1$，即表明
$$
F_Y(x) = F_{X}(x)
$$
也就是说这两个随机变量的的累积分布函数是一样的，即这两个随机变量是同分布的。

举个例子，假设要构造一个随机变量服从概率密度函数$f_{X}(x) = x^2(x\in(0,\sqrt[3]{3}))$，步骤如下

  • 求出累积分布函数$F_{X}(x) = \frac{1}{3}x^3(x\in(0,\sqrt[3]3))$
  • 求出累积分布函数的逆函数为$F_{X}^{-1} = (3x)^{-3}$
  • 用一个0到1的均匀分布的随机变量$U(0,1)$，则变量$(3U)^{-3}$就会服从目标分布

代码段如下

1
2
3
4
std::uniform_real_distribution<> values {0.0, 1.0};
std::random_device rd; 
std::default_random_engine rng {rd()}; 
double X = pow((3 * values(rng)), -3);

构造正态分布变量

在我搜集到的资料当中，一般有三种方法，分别是The Box–Muller transform，以及它的改良版Marsaglia polar method，还有一种The Ziggurat algorithm。这里主要介绍前两种，因为跟前面介绍的逆变换方法有关。

The Box–Muller transform

标准正态分布的概率密度函数为$f(x) = \frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}$，这个函数直接去积分求它的累积分布函数是比较麻烦的，但是我们可以同时构造两个互相独立的标准正态分布变量。比如要构造两个独立的标准正态分布变量$X$和$Y$，那么这两个随机变量的联合概率密度函数$f(x,y)$即为
$$
f(x,y) = \frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}} *\frac{1}{\sqrt{2\pi}} e^{-\frac{y^2}{2}}\\
f(x,y) = \frac{1}{2\pi}e^{-\frac{x^2+y^2}{2}}
$$
这个联合分布的累积分布函数是好求的，那么它的积分（$I^2$）就为
$$
\iint_{-\infty} ^{+\infty}\frac{1}{2\pi} e^{-\frac{x^2+y^2}{2}} dxdy
$$
这是在笛卡尔直角坐标系下的表达，可以转成极坐标下的表达（注意$dxdy$是和$rdrd\theta$等价的，不要忘了那个$r$，忘记的可以看看雅可比矩阵）
$$
I^{2} = \iint \frac{1}{2\pi}e^{-\frac{r^2}{2}}*rdrd\theta
$$
我们可以看出，这个被积分函数是跟$\theta$没关系的，也就是说$\theta$和$r$是也是独立的，上述分布可以写成
$$
I^2 = \int_{0}^{2\pi}\frac{1}{2\pi}d\theta * \int_{0}^{+\infty}re^{-\frac{r^2}{2}}dr
$$
左边关于$\theta$很好解决，右边的关于$r$单独拿出来，这个原函数也是可以求出来的，这两个变量的累积分布函数即为
$$
F_R(r) = 1-e^{-\frac{r^2}{2}} \quad r\in(0,+\infty) \\
F_{\Theta}(\theta) = \frac{\theta}{2\pi} \quad \theta\in(0,2\pi)
$$
根据前面说的逆变换方法，这两个变量是可以通过两个0到1的随机变量$U_1$和$U_2$构造出来的，求出上述两个分布函数的逆函数，再代入$U_1$和$U_2$，即可以得到$R$跟$\Theta$两个随机变量
$$
R = \sqrt{-2\ln(1-U_1)} \\
\Theta = 2\pi U_2
$$
注意我们是要构造$X$和$Y$这两个随机变量，极坐标系跟笛卡尔直角坐标系有如下转换关系
$$
X = R * \cos(\Theta) \\
Y = R * \sin(\Theta)
$$
即
$$
X = \sqrt{-2\ln(1-U_1)}*\cos(2\pi U_2) \\
Y = \sqrt{-2\ln(1-U_1)}*\sin(2\pi U_2)
$$
其中因为$U_1$是均匀分布的，也就是说$1-U_1$也是均匀分布，所以可以直接用$U_1$替代，即简化为
$$
X = \sqrt{-2\ln U_1}*\cos(2\pi U_2) \\
Y = \sqrt{-2\ln U_1}*\sin(2\pi U_2)
$$

Marsaglia polar method

上述的方法还是要计算三角函数，在Marsaglia polar method里就可以略去这个过程。

其实主要问题在这两个三角函数，这是我们不希望有的，那么有没有一种变量的分布刚好是$\cos(2\pi U_2)$或者$\sin(2\pi U_2)$这样的分布呢，可能你的第一想法还是用逆变换方法去做，但是会马上意识到，这样做会陷入一种三角函数的循环。那我们来试一试，假设你要构造的变量是$K=\cos(2\pi U_2)$，那么它的累积分布函数为
$$
F_K(k) = P(\cos(2\pi U_2)<k)
$$
如图所示

则上述等式可写为
$$
F_K(k) = \frac{2\pi - 2\cos^{-1}(k)}{2\pi}
$$
这个函数的反函数依然是三角函数，并没有产生什么实质性的变化，所以继续用逆函数变换方法去做是取消不掉三角函数的运算的。

Marsaglia polar method是先产生两个在-1到1之间的均匀分布的变量$U_3 U_4$，这样的话在坐标系中就是一个边长为2的正方形，紧接着再剔除那些落在单位圆之外的点，这样就得到了一个在单位圆内均匀分布的点。在这个单位圆内，存在一个均匀分布的变量（注意这个时候的横纵坐标由于筛选，并不服从均匀分布了），那就是极角，这一点很好理解。这个均匀分布极角有个很好的特性，就是它的三角函数值是可以用横纵坐标，也就是之前随机出来的两个变量表示出来。注意这时候的$U_3 U_4$只是为了表示变量值，他们已经失去了均匀分布的特性
$$
\cos(2\pi U_2) = \cos(\Theta) = \frac{U_3}{\sqrt{U_3^2+U_4^2}}\\
\sin(2\pi U_2) = \sin(\Theta) = \frac{U_4}{\sqrt{U_3^2+U_4^2}}
$$
其实到这里我们已经解决了要计算三角函数的问题，但是我们用了三个均匀分布的变量，即$U_1 U_3 U_4$（$U_1$是0到1的均匀分布，$U_3U_4$是-1到1的均匀分布），这不是我们想要的，那么上面构造出来的单位圆里面，还存不存在一个在0到1之间均匀分布而且相对于$\Theta$独立的变量呢？其实还是有的，就是点的到原点距离的平方$S = U_3^2 + U_4^2$，且$S$跟极角$\Theta$也是独立的。

同样的考虑积分（因为在单位圆里均匀分布，所以被积函数是一个常数$\frac{1}{\pi}$，单位圆面积的倒数），根据雅可比矩阵
$$
\iint \frac{1}{\pi}dxdy = \frac{1}{\pi}\iint rdrd\theta\\
\iint \frac{1}{\pi}dxdy = \frac{1}{\pi}\iint \frac{1}{2}dr^2d\theta \\
\iint \frac{1}{\pi}dxdy = \frac{1}{2\pi}\int_{0}^{2\pi}d\theta\int_{0}^{1}dr^2
$$
也就是说不光是$\Theta$，半径的平方$r^2$，即$S$也是均匀分布的，而且是在0到1上的均匀分布，而且还是独立的（因为上面的概率密度函数可以分开相乘）。

那么在Box–Muller方法中得出的两组正态分布的变量就可以写为
$$
X = \sqrt{-2\ln S}*\frac{U_3}{\sqrt{S}} \\
Y = \sqrt{-2\ln S}*\frac{U_4}{\sqrt{S}}
$$
c++标准函数库里std::normal_distribution用的就是Marsaglia polar method，下面是具体的代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
if (_M_saved_available)//因为一次性可以生成两个正态分布的变量，所以可以利用起来
{
  _M_saved_available = false;
  __ret = _M_saved;
}
else
{
  result_type __x, __y, __r2;
  do
  {
    __x = result_type(2.0) * __aurng() - 1.0;
    __y = result_type(2.0) * __aurng() - 1.0;
    __r2 = __x * __x + __y * __y;
  }
  while (__r2 > 1.0 || __r2 == 0.0);//筛去单位圆以外以及等于0的

  const result_type __mult = std::sqrt(-2 * std::log(__r2) / __r2);
  _M_saved = __x * __mult;
  _M_saved_available = true;
  __ret = __y * __mult;
}