花式生成正态分布

发表于编程 2016-11-29 0:12 字数: 7121 阅读量: 1,629

◆ ◆ ◆

前言

“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<x[i+1]，y必然在曲线下方，满足条件，接受x；若x[i+1]<x<x[i]，则还需要进一步判断。同样这里R[0]和tail中的样本需要进行特殊处理

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

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