生成正态分布的随机数

2018-02-09


正态分布是非常常见的一种随机分布。本文就来讲讲几个生成正态分布随机数的算法。

下文都是以生成标准正态分布,也就是数学期望为0,标准差为1的正态分布为例。

相信大家都在概率论课上学过了,标准正态分布的概率密度函数\( f(x) \)是:

\[ f(x) = {{1 \over \sqrt{2 \pi}} exp(-{x^2 \over 2})} \]


Reject Sampling

Reject Sampling (拒绝采样) 应该是坠简单的,坠好理解的生成正态分布的随机数的算法了(事实上不仅仅是适用于正态分布)。我们首先在给定的区间里生成一个随机数(均匀分布的)。根据概率密度函数,我们可以得到这一随机分布在这一取值点取值的概率。比如生成的随机数是0.1,根据概率密度函数,我们可以认为在这一点取值的概率是0.3969525474770118。那么我们再根据这一概率再进行一次随机,让程序有如上概率接受这一值并输出,否则就拒绝这个值,重新进行采样。代码如下:

def normal_probability(x):
    return (1 / math.sqrt(2 * math.pi)) * math.exp(x * x / -2)

def reject_sampling():
    while True:
        num = random.uniform(-500, 500)
        if (random.uniform(0, 1) < normal_probability(num)):
            return num

当然,弊端也很明显,迭代次数会非常多。区间越大,效率也越低。


Box Muller

Box Muller就是一个比较简洁的算法了。它只需要一些初等函数就能快速生成。在区间(0, 1)上生成俩随机数U和V,然后根据下列公式获得正态分布的随机数:

\[ X = \sqrt{−2\log(U)}\, \cos(2\pi V) \quad \text{or} \quad X = \sqrt{−2\log(U)}\, \sin(2\pi V)\ \]

早期这是计算机里通用的生成正态分布随机数的方法。使用极坐标的方法可以对这一公式进行证明。手打有点累,可参考:https://math.stackexchange.com/questions/1110168/proof-of-the-box-muller-method