生成特定分布随机数的方法

2017 年 9 月 29 日 算法与数据结构 张洋

来源：CodingLabs-张洋

blog.codinglabs.org/articles/methods-for-generating-random-number-distributions.html

生成随机数是程序设计里常见的需求。一般的编程语言都会自带一个随机数生成函数，用于生成服从均匀分布的随机数。不过有时需要生成服从其它分布的随机数，例如高斯分布或指数分布等。有些编程语言已经有比较完善的实现，例如Python的NumPy。这篇文章介绍如何通过均匀分布随机数生成函数生成符合特定概率分布的随机数，主要介绍Inverse Ttransform和Acceptance-Rejection两种基础算法以及一些相关的衍生方法。下文我们均假设已经拥有一个可以生成0到1之间均匀分布的随机数生成函数，关于如何生成均匀分布等更底层的随机数生成理论，请参考其它资料，本文不做讨论。

基础算法

Inverse Transform Method

最简单的生成算法是Inverse Transform Method（下文简称ITM）。如果我们可以给出概率分布的累积分布函数（下文简称CDF）及其逆函数的解析表达式，则可以非常简单便捷的生成指定分布随机数。

ITM算法描述

生成一个服从均匀分布的随机数U∼Uni(0,1)
设F(X)为指定分布的CDF，F−1(Y)是其逆函数。返回X=F−1(U)作为结果

ITM算法说明

这是一个非常简洁高效的算法，下面说明其原理及正确性。

我们通过图示可以更直观的明白算法的原理。下图是某概率分布的CDF：

我们从横轴上标注两点xa和xb，其CDF值分别为F(xa)和F(xb)。

由于U服从0到1之间的均匀分布，因此对于一次U的取样，U落入F(xa)和F(xb)之间的概率为:

而由于CDF都是单调非减函数，因此这个概率同时等于X落入xa和xb之间的概率，即：

而根据CDF的定义，这刚好说明XX服从以F(x)为CDF的分布，因此我们的生成算法是正确的。

ITM实现示例

下面以指数分布为例说明如何应用ITM。

首先我们需要求解CDF的逆函数。我们知道指数分布的CDF为

通过简单的代数运算，可以得到其逆函数为

由于UU服从从0到1的均匀分布蕴含着1−U服从同样的分布，因此在实际实现时可以用Y代替1−Y，得到：

下面给出一个Python的实现示例程序。

import random

import math

def exponential_rand(lam):

if lam <= 0:

return -1

U = random.uniform(0.0, 1.0)

return (-1.0 / lam) * math.log(U)

Acceptance-Rejection Method

一般来说ITM是一种很好的算法，简单且高效，如果可以使用的话，是第一选择。但是ITM有自身的局限性，就是要求必须能给出CDF逆函数的解析表达式，有些时候要做到这点比较困难，这限制了ITM的适用范围。

当无法给出CDF逆函数的解析表达式时，Acceptance-Rejection Method（下文简称ARM）是另外的选择。ARM的适用范围比ITM要大，只要给出概率密度函数（下文简称PDF）的解析表达式即可，而大多数常用分布的PDF是可以查到的。

ARM算法描述

设PDF为f(x)f(x)。首先生成一个均匀分布随机数X∼Uni(xmin,xmax)
独立的生成另一个均匀分布随机数Y∼Uni(ymin,ymax)
如果Y≤f(X)，则返回X，否则回到第1步

ARM算法说明

通过一幅图可以清楚的看到ARM的工作原理。

ARM本质上是一种模拟方法，而非直接数学方法。它每次生成新的随机数后，通过另一个随机数来保证其被接受概率服从指定的PDF。

显然ARM从效率上不如ITM，但是其适应性更广，在无法得到CDF的逆函数时，ARM是不错的选择。

ARM实现示例

下面使用ARM实现一个能产生标准正态分布的随机数生成函数。

首先我们要得到标准正态分布的PDF，其数学表示为：

为了方便，这里我会直接使用SciPy来计算其PDF。

程序如下。

import random

import scipy.stats as ss

def standard_normal_rand():

while True:

X = random.uniform(-3.0,3.0)

Y = random.uniform(0.0, 0.5)

if Y < ss.norm.pdf(X):

return X