1. 程式人生 > >正態分佈隨機數生成演算法

正態分佈隨機數生成演算法

最近在學習基於蒙特卡羅的強化學習方法時遇到 生成服從正態分佈的隨機數的演算法,因此做一個回顧和總結。

要程式設計得到服從均勻分佈的偽隨機數是容易的,C、Python、Java語言等都提供了相應的函式。

但是要想生成服從正態分佈的隨機數就沒那麼容易了,生成服從正態分佈的隨機數的基本思想是先得到服從均勻分佈的隨機數,再將服從均勻分佈的隨機數轉變為服從正態分佈。

實現均勻分佈到正態分佈轉變的方法:

  • 利用分佈函式的反函式

使用反函式,先隨機抽出一個服從[0,1]均勻分佈的數字u,然後

那z是正態分佈的

import numpy as np  
from scipy.special import erfinv 
def inverfsampling(mu=0, sigma=1, size=1):  
    z = np.sqrt(2) * erfinv(2 * np.random.uniform(size=size) - 1)  
    return mu + z * sigma
  • Box Muller方法
Box Muller方法的推導過程較為複雜。但得到的結果卻是很令人滿意的。只要用兩個相互獨立的均勻分佈就能得到正態分佈用Box-Muller方法,隨機抽出兩個從均[0,1]勻分佈的數字u和v。然後




均服從正態分佈。

import numpy as np  

def boxmullersampling(mu=0, sigma=1, size=1):  
    u = np.random.uniform(size=size)  
    v = np.random.uniform(size=size)  
    z = np.sqrt(-2 * np.log(u)) * np.cos(2 * np.pi * v)  
    return mu + z * sigma  
  
  • 利用中心極限定理
利用林德伯格-萊維(Lindeberg—Levi)中心極限定理, 獨立同分布的事件,具有相同的期望和方差,則事件服從中心極限定理。他表示了對於抽取樣本,n足夠大的時候,樣本分佈符合x~N(μ,σ^2)。
中心極限定理告訴我們,當樣本量足夠大時,樣本均值的分佈慢慢變成正態分佈。
import numpy as np
import matplotlib.pylab as plb
from scipy import stats
import math

number=[30,100,300,1000,5000,30000]
y=[]
for i in range(len(number)):
    ave_uniform=[]
    for j in range(1000):
        data=np.random.uniform(-1,1,number[i])
        variance=(1.0/3.0)/(1.0*number[i])
        summer=sum(data)
        ave_uniform.append(summer/(1.0*number[i]))
        Range=np.arange(-0.5,0.5,0.001)
    plb.subplot(230+i+1)
    plb.plot(Range,stats.norm.pdf(Range,0,math.sqrt(variance)))
    plb.hist(ave_uniform,bins=100)
plb.show()
關於正態分佈的理論證明,可參考課程:http://open.163.com/movie/2011/6/G/I/M82IC6GQU_M83JBFVGI.html