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

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

Inverse Transform Method

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

ITM算法描述

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

ITM算法说明

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

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

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

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

\[P(U \in (F(x_a),F(x_b))) = F(x_b) - F(x_a)\]

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

\[P(U \in (F(x_a),F(x_b))) = P(F^{-1}(U) \in (F^{-1}(F(x_a)),F^{-1}(F(x_b)))) = P(X \in (x_a,x_b)) = F(x_b) - F(x_a)\]

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

ITM实现示例

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

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

\[F(X)=1-e^{-\lambda X}\]

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

\[F^{-1}(Y)=-\frac{1}{\lambda}\ln(1-Y)\]

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

\[F^{-1}(Y)=-\frac{1}{\lambda}\ln(Y)\]

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

import random
import math
def exponential_rand(lam):
    if lam