◆◆ ◆

前言

“So much of life, it seems to me, is determined by pure randomness.” – Sidney Poitier

在之前的博客中，或多或少都提到过一些随机、伪随机、随机数等等，但基本上只是直接使用，没有探寻背后的一些原理，刚好最近偶然看到python标准库中如何生成服从正态分布随机数的源码，所以本文就简单聊聊如何生成正态分布

◆◆ ◆

知识回顾

为了防止你对接下来的内容一头雾水，我觉得还是有必要回顾一下我们曾经学过的高数和概率统计知识

均匀分布

不过这些约束怎么来的本文就暂不讨论了

我们以Java中Random的实现为例，为了便于理解，这里对代码进行了部分调整，只截取了相关的片段

//这两个东西在本文后面讲正态分布的时候会涉及到

private double nextNextGaussian;

private boolean haveNextNextGaussian = false;

//随机数序列生成器种子

private final AtomicLong seed;

//线性同余发生器乘数a

private final static long multiplier = 0x5DEECE66DL;

//线性同余发生器加数c

private final static long addend = 0xBL;

//线性同余发生器模数m

private final static long mask = (1L << 48) - 1;

public Random(long seed) {

this.seed = new AtomicLong(0L);

setSeed(seed);

}

synchronized public void setSeed(long seed) {

//实现线性同余算法

seed = (seed ^ multiplier) & mask;

//atomic set具有写入(分配)volatile 变量的内存效果(即具有内存可见性)

this.seed.set(seed);

//暂且忽略。。

haveNextNextGaussian = false;

}

protected int next(int bits) {

long oldseed, nextseed;

AtomicLong seed = this.seed;

//ompareAndSet 如果当前值==预期值，则以原子方式将该值设置为给定的更新值(利用了CPU的硬件原语CAS指令)

do {

//atomic get具有读取volatile 变量的内存效果

oldseed = seed.get();

//实现线性同余算法

nextseed = (oldseed * multiplier + addend) & mask;

} while (!seed.compareAndSet(oldseed, nextseed));

//截取bits位整型

return (int)(nextseed >>> (48 - bits));

}

有人可能会有疑惑，这个代码中的实现nextseed = (oldseed * multiplier + addend) & mask好像和递推公式不一样啊？那个模运算为什么变成了与运算？

注意，x&[(1L << 48)–1]与x(mod 2^48)是等价的。为什么呢？从二进制的角度来考虑这个问题就很清楚了。一个数x除以2^n，在二进制中相当于将x右移n位，商和余数分别在小数点左侧和右侧。如23/8=2，23%8=7，23的二进制表示为10111，除以2^3=8相当于右移3位，得到10.111，左侧为商10也就是2，右侧为余数111也就是7。也就是说如果一个数对2^n取余，那么只需要得到该数的低n位即可，很自然的想到如果能得到这样

一个二进制数，与原来的数做与运算即可，而2^n-1恰好可以得到这样的一个二进制数。如2^3-1=7→111，2^7-1=127→1111111。

现在回头看看nextseed = (oldseed * multiplier + addend) & mask这行代码可以理解了吧？

关于均匀分布的更多细节可以参考JDK源码分析——从java.util.Random源码分析线性同余算法

概率分布函数和概率密度函数

直接引用教材上的黑体字吧

看下面这张图也非常清楚，概率分布函数F(x)是cumulative distribution function(CDF)，概率密度函数f(x)是probability density function(PDF)

刚才说到的均匀分布概率密度函数f(x)如下

说人话就是，n个相互独立同分布的随机变量之和的分布近似于正态分布，nn越大，近似程度越好

# -*- coding: utf-8 -*-

import matplotlib.pyplot as plt

import numpy as np

def getExponential(SampleSize,p_lambda):

result = -np.log(1-np.random.uniform(0,1,SampleSize))/p_lambda

return result

# 生成10000个数，观察它们的分布情况

SampleSize = 10000

es = getExponential(SampleSize, 1)

plt.hist(es,np.linspace(0,5,50),facecolor="green")

plt.show()

得到如下结果

对比维基百科里面标准的指数分布

雅可比矩阵与雅可比行列式

这个东西在高数课本中有，只怪当初学习不用功……

以咱们熟悉的平面直角坐标与极坐标转换为例吧，

回顾当时学习二重积分时，利用极坐标变换时的式子如下

现在知道那个多出来的r是怎么回事了吧？雅可比行列式相当于两个不同坐标系中微小区域面积的缩放倍数

拒绝采样(Rejection Sampling)

这个方法有的时候也称接收-拒绝采样，使用场景是有些函数p(x)太复杂在程序中没法直接采样，那么可以设定一个程序可抽样的分布q(x)比如正态分布等等，然后按照一定的方法拒绝某些样本，达到接近p(x)分布的目的：

具体操作如下，设定一个方便抽样的函数q(x)，以及一个常量kk，使得p(x)总在kq(x)的下方。(参考上图)

x轴方向：从q(x)分布抽样得到a

y轴方向：从均匀分布(0,kq(a))中抽样得到u

如果刚好落到灰色区域：u>p(a)，拒绝；否则接受这次抽样

重复以上过程

证明过程就不细说了，知道怎么用就行了，感兴趣的可以看看这个文档

Acceptance-Rejection Method

不过在高维的情况下，拒绝采样会出现两个问题，第一是合适的q分布比较难以找到，第二是很难确定一个合理的k值。这两个问题会造成图中灰色区域的面积变大，从而导致拒绝率很高，无用计算增加。

◆◆ ◆

暴力生成正态分布

根据中心极限定理，生成正态分布就非常简单粗暴了，直接生成n个独立同分布的随机变量，求和即可。注意，无论你使用什么分布，当n趋近于无穷大时，它们和的分布都会趋近正态分布！

以最简单的均匀分布为例，看代码

# -*- coding: utf-8 -*-

import matplotlib.pyplot as plt

import numpy as np

def getNormal(SampleSize,n):

xsum = []

for i in range(SampleSize):

# 利用中心极限定理，[0,1)均匀分布期望为0.5，方差为1/12

tsum = (np.mean(np.random.uniform(0,1,n))-0.5)*np.sqrt(12*n)

xsum.append(tsum)

return xsum

# 生成10000个数，观察它们的分布情况

SampleSize = 10000

# 观察n选不同值时，对最终结果的影响

N = [1,2,10,1000]

plt.figure(figsize=(10,20))

subi = 220

for index,n in enumerate(N):

subi += 1

plt.subplot(subi)

normalsum = getNormal(SampleSize, n)

# 绘制直方图

plt.hist(normalsum,np.linspace(-4,4,80),facecolor="green",label="n={0}".format(n))

plt.ylim([0,450])

plt.legend()

plt.show()

得到结果如下图所示

可以看到，n=1时其实就是均匀分布，n=2时有正态分布的样子了，但不够平滑，随着n逐渐增大，直方图轮廓越来越接近正态分布了~因此利用中心极限定理暴力生成服从正态分布的随机数是可行的

但是这样生成正态分布速度是非常慢的，因为要生成若干个同分布随机变量，然后求和、计算，效率非常低。那有没有其他办法呢？

当然有！利用反变换法

◆◆ ◆

反变换法生成正态分布

正态分布的概率分布函数(CDF)如下图所示，

在y轴上产生服从(0,1)均匀分布的随机数，水平向右投影到曲线上，然后垂直向下投影到x轴，这样在x轴上就得到了正态分布。

当然正态分布的概率分布函数不方便直接用数学形式写出，求反函数也无从说起，不过好在scipy中有相应的函数，我们直接使用即可

# -*- coding: utf-8 -*-

import matplotlib.pyplot as plt

import numpy as np

from scipy.stats import norm

def getNormal(SampleSize):

iid = np.random.uniform(0,1,SampleSize)

result = norm.ppf(iid)

return result

SampleSize = 10000000

normal = getNormal(SampleSize)

plt.hist(normal,np.linspace(-4,4,81),facecolor="green")

plt.show()

结果如下图所示，

以上两个方法虽然方便也容易理解，但是效率实在太低，并不实用，那么在实际中到底是如何生成正态分布的呢？

◆◆ ◆

Box–Muller算法

说来也巧，某天闲的无聊突然很好奇python是如何生成服从正态分布的随机数的，于是就看了看random.py的代码，具体实现的代码就不贴了，大家自己去看，注释中有下面几行

# When x and y are two variables from [0, 1), uniformly

# distributed, then

#

#    cos(2*pi*x)*sqrt(-2*log(1-y))

#    sin(2*pi*x)*sqrt(-2*log(1-y))

#

# are two *independent* variables with normal distribution

写程序实现一下

# -*- coding: utf-8 -*-

import matplotlib.pyplot as plt

import numpy as np

def getNormal(SampleSize):

iid = np.random.uniform(0,1,SampleSize)

normal1 = np.cos(2*np.pi*iid[0:SampleSize/2-1])*np.sqrt(-2*np.log(iid[SampleSize/2:SampleSize-1]))

normal2 = np.sin(2*np.pi*iid[0:SampleSize/2-1])*np.sqrt(-2*np.log(iid[SampleSize/2:SampleSize-1]))

return np.hstack((normal1,normal2))

# 生成10000000个数，观察它们的分布情况

SampleSize = 10000000

es = getNormal(SampleSize)

plt.hist(es,np.linspace(-4,4,80),facecolor="green")

plt.show()

得到的结果如下图所示，

这里抽样次数达到1千万次，1秒左右就完成了，速度比暴力生成正态分布要快的多~

ps:由于Box–Muller算法一次性生成了两个独立且服从正态分布的随机数，所以可以把其中一个保存起来，下次直接使用即可。本文刚开始的那段代码中nextNextGaussian就是用来保存它的~

◆◆ ◆

Ziggurat Algorithm

Box–Muller算法虽然快了许多，但是由于用到了三角函数和对数函数，相对来说还是比较耗时的，如果想要更快一点有没有办法呢？

当然有，这就是Ziggurat算法，不仅可以用于快速生成正态分布，还可以生成指数分布等等。其基本思想就是利用拒绝采样，其高效的秘密在于构造了一个非常精妙的q(x)，看下面这张图

如果为了方便，我们当然可以直接使用一个均匀分布，也就是一个矩形，但是这样的话，矩形与正态分布曲线间的距离很大，容易造成拒绝率很高，无用计算增加，高效也就无从谈起了

再看看上面那张图，我们用多个堆叠在一起的矩形，这样保证阴影部分(被拒绝部分)的始终较小，这样就非常高效了

简单对图作一个解释：

我们用R[i]来表示一个矩形，R[i]的右边界为x[i]，上边界为y[i]。S[i]表示一个分割，当i=0时，S[0]=R[0]+tail，其他情况S[i]==R[i]

除了R[0]以外，其他每个矩形面积相等，设为定值A。R[0]的面积=A-tail的面积。这样保证从任何一个分割中抽样(x,y)的概率相等

当任意选定一个R[i]在其中抽样(x,y)，若x

这里为了方便解释，只用了几个矩形，在程序实现的时候，一般使用128或256个矩形

可以看出，为了提高效率，Ziggurat算法中使用了许多技巧性的东西，这在其C代码实现中更加明显，使用了与运算和字节等各种小技巧，代码就不在这里贴了，感兴趣可以看看下面几个版本，C版本的追求的是极致的速度，每个矩形的边界已经提前计算好了。C#版本中的注释非常详细，Java版的基本与C#一致，但是效率一般。

C

C#

Java

最后对比一下Ziggurat算法与Box-muller算法的效率

◆◆ ◆

总结

本文介绍了多种生成正态分布的方法，其中Box-muller算法应对一般的需求足够了，但是要生成大量服从正态分布的随机数时，Ziggurat算法效率会更高一点~

再说点题外话，作为一名普通的程序员，对于很多东西往往不需要了解的非常深入，说白了“会用就行了”。但是有的时候探寻其背后的原理往往能发现别人领会不到的数学之美，这也是写程序之余的一点乐趣吧~

◆◆ ◆

参考资料

大数定律与中心极限定理

http://ac-cf2bfs1v.clouddn.com/71ErfPP72vdJRcCTIcC9iOVnAPEOVFDveKI1JMO1.pdf

Jacobian 矩陣與行列式

https://ccjou.wordpress.com/2012/11/26/jacobian-%E7%9F%A9%E9%99%A3%E8%88%87%E8%A1%8C%E5%88%97%E5%BC%8F/

从随机过程到马尔科夫链蒙特卡洛方法

http://www.cnblogs.com/daniel-D/p/3388724.html

随机数产生原理

http://ac-cf2bfs1v.clouddn.com/uPcW0ce2E0FIDIt53mLwHGJ5s6xTadE4mqVCpsWd.ppt

Transformed Random Variables

http://www.mathematik.uni-ulm.de/numerik/teaching/ss09/NumFin/Script/chap2_4-2_5.pdf

Box-Muller Transform

Normalityhttp://math.stackexchange.com/questions/1005236/box-muller-transform-normality

The Ziggurat Method for Generating Random Variables

http://ac-cf2bfs1v.clouddn.com/nNPPeDafleeRW3U073UA4mPbYiAxhgb0soziU5Uo.pdf

The Ziggurat Algorithm for Random Gaussian Sampling

http://heliosphan.org/zigguratalgorithm/zigguratalgorithm.html

作者：Bindog

原文链接：http://bindog.github.io/blog/2016/06/04/from-sne-to-tsne-to-largevis

推荐阅读

知识无价，感谢作者分享！扫一扫向作者打赏~

专业大数据竞赛平台

中国数据青年成长之家