如果说最常用的分布是哪个,第一个当属均匀分布,第二个就非高斯分布莫属了,在机器学习中常见的采样方法中介绍了均匀分布的采样,本文来介绍一下高斯分布的采样.

首先,假设随机变量z服从标准正态分布N(0,1),令:

x = \sigma \cdot z + \mu

x服从均值为\mu,方差为\sigma^2的高斯分布N(\mu,\sigma^2).因此,任意高斯分布都可以由标准正态分布通过拉伸和平移得到,所以这里只考虑标准正态分布的采样.

由于高斯分布的累积分布函数并不是一个初等函数,所以使用逆变换法计算比较麻烦,Box-Muller算法给出了简单的解决方案.

Box-Muller算法

假设x,y是两个服从标准正态分布的独立随机变量,他们的联合概率密度为

p ( x , y ) = \frac { 1 } { 2 \pi } e ^ { - \frac { x ^ { 2 } + y ^ { 2 } } { 2 } }

考虑(x,y)在圆盘\{ ( x , y ) | x ^ { 2 } + y ^ { 2 } \leqslant R ^ { 2 } \}上的概率

F ( R ) = \int _ { x ^ { 2 } + y ^ { 2 } \leqslant R ^ { 2 } } \frac { 1 } { 2 \pi } e ^ { - \frac { x ^ { 2 } + y ^ { 2 } } { 2 } } d x d y

通过极坐标变换将(x,y)转化为( r , \theta ),可以很容易求得二重积分:

F ( R ) = 1 - \mathrm { e } ^ { - \frac { R ^ { 2 } } { 2 } } , R \geqslant 0

这里F(R)可以看成是极坐标中r的累积分布函数.由于F(R)的计算公式比较简单,逆函数也容易求得,所以可以利用逆变换法来对r进行采样;对于\theta,在[0,2\pi]上进行均匀采样即可.这样就得到了(r,\theta),经过坐标变换即可得到符合标准正态分布的(x,y).具体采样过程如下:

  1. 产生[0,1]上的两个独立的均匀分布随机数u_1,u_2
  2. \left\{ \begin{array} { l } { x = \sqrt { - 2 \ln \left( u _ { 1 } \right) } \cos 2 \pi u _ { 2 } } \\ { y = \sqrt { - 2 \ln \left( u _ { 1 } \right) } \sin 2 \pi u _ { 2 } } \end{array} \right.,则x,y服从标准正态分布,并且是相互独立的.

Box-muller算法由于需要计算三角函数,相对来说还是比较耗时,而Marsaglia polar method则避开了三角函数的计算,因而更快,其具体采样操作如下:

  1. 在单位圆盘\{ ( x , y ) | x ^ { 2 } + y ^ { 2 } \leq 1 \}上产生均匀分布随机数对(x,y)(在矩阵\{ ( x , y ) | - 1 \leqslant x , y \leqslant 1 \}上利用拒绝采用法即可得到).
  2. s = x ^ { 2 } + y ^ { 2 },则x \sqrt { \frac { - 2 \ln s } { s } } , y \sqrt { \frac { - 2 \ln s } { s } }是两个服从标准正态分布的样本,其中\frac { x } { \sqrt { s } } , \frac { y } { \sqrt { s } }用来代替Box-muller算法中的cos和sin操作.

编程实现

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:int(SampleSize/2-1)])*np.sqrt(-2*np.log(iid[int(SampleSize/2):SampleSize-1]))
    normal2 = np.sin(2*np.pi*iid[0:int(SampleSize/2-1)])*np.sqrt(-2*np.log(iid[int(SampleSize/2):SampleSize-1]))
    return np.hstack((normal1,normal2))

# 生成10000000个数,观察它们的分布情况
SampleSize = 10000000
normal = getNormal(SampleSize)
plt.hist(normal,np.linspace(-4,4,81),facecolor="green")
plt.show()

posted @ 2019-01-22 15:22:46
评论加载中...

发表评论