漫談正態分佈的生成

本文作者簡介：王夜笙，就讀於鄭州大學資訊工程學院，感興趣的方向為逆向工程和機器學習，長期從事資料抓取工作（長期與反爬蟲技術作鬥爭~），涉獵較廣（技藝不精……），詳情請見我的個人部落格~

感謝怡軒同學的悉心指導~

之前拜讀了靳志輝（@rickjin）老師寫的《正態分佈的前世今生》，一直對正態分佈懷著一顆敬畏之心，剛好最近偶然看到python標準庫中如何生成服從正態分佈隨機數的原始碼，覺得非常有趣，於是又去查詢其他一些生成正態分佈的方法，與大家分享一下。

設x1,x2,⋯,xn

為獨立同分布的隨機變數序列，均值為μ

，方差為σ2

，則zn=x1+x2+⋯+xn−nμσn√

具有漸近分布n(0,1)

，也就是說當n→∞

時，p→12π−−√∫x−∞e−t22dt

換句話說，n

個相互獨立同分布的隨機變數之和的分布近似於正態分佈，n

越大，近似程度越好。當然也有例外，比如n

個獨立同分布的服從柯西分布隨機變數的算術平均數仍是柯西分布，這裡就不擴充套件講了。

根據中心極限定理，生成正態分佈就非常簡單粗暴了，直接生成n個獨立同分布的均勻分布即可，看**

# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import numpy as np
def getnormal(samplesize,n):
result = 
for i in range(samplesize):
# 利用中心極限定理，[0,1)均勻分布期望為0.5，方差為1/12
iid = (np.mean(np.random.uniform(0,1,n))-0.5)*np.sqrt(12*n)
return result
# 生成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)
normal = getnormal(samplesize, n)
# 繪製直方圖
plt.hist(normal,np.linspace(-4,4,81),facecolor="green",label="n=".format(n))
plt.ylim([0,450])
plt.legend()
plt.show()

得到結果如下圖所示

可以看到，n=1時其實就是均勻分布，隨著n逐漸增大，直方圖輪廓越來越接近正態分佈了~因此利用中心極限定理暴力生成服從正態分佈的隨機數是可行的。但是這樣生成正態分佈速度是非常慢的，因為要生成若干個同分布隨機變數，然後求和、計算，效率是非常低的。

假設u=f(x)

是乙個概率分布函式(cdf)，f−1

是它的反函式，若u

是乙個服從(0,1)

均勻分布的隨機變數，則f−1(u)

服從函式f

給出的分布。例如要生成乙個服從指數分布的隨機變數，我們知道指數分布的概率分布函式(cdf)為f(x)=1–e–λx

，其反函式為f−1(x)=−ln(1−x)λ

，所以只要不斷生成服從(0,1)

均勻分布的隨機變數，代入到反函式中即可生成指數分布。

正態分佈的概率分布函式(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()

結果如下圖所示，

以上兩個方法雖然方便也容易理解，但是效率實在太低，並不實用，那麼在實際中到底是如何生成正態分佈的呢？

說來也巧，某天閒的無聊突然很好奇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

頓時感覺非常神奇，也就是說當x

和y是兩個獨立且服從(0,1)

均勻分布的隨機變數時，cos(2πx)⋅–2ln(1–y)−−−−−−−−−√

和sin(2πx)⋅–2ln(1–y)−−−−−−−−−√

是兩個獨立且服從正態分佈的隨機變數！

後來查了查這個公式，發現這個方法叫做box–muller，其實本質上也是應用了逆變換法，證明方法比較多，這裡我們選取一種比較好理解的

我們把逆變換法推廣到二維的情況，設u1,u2

為(0,1)

上的均勻分布隨機變數，(u1,u2)

的聯合概率密度函式為f(u1,u2)=1(0≤u1,u2≤1)

，若有：

{u1=g1(x,y)u2=g2(x,y)

其中，g1,g2

的逆變換存在，記為

{x=h1(u1,u2)y=h2(u1,u2)

且存在一階偏導數，設j

為jacobian矩陣的行列式

|j|=∣∣∣∣∂x∂u1∂y∂u1∂x∂u2∂y∂u2∣∣∣∣≠0

則隨機變數(x,y)

的二維聯合密度為（回顧直角座標和極座標變換）：

f[h1(u1,u2),h2(u1,u2)]⋅|j|=|j|

根據這個定理我們來證明一下，

{y1=−2lnx1−−−−−−−√cos(2πx2)y2=−2lnx1−−−−−−−√sin(2πx2)

求反函式得

⎧⎩⎨x1=e–y21+y222x2=12πarctany2y1

計算jacobian行列式

|j|=∣∣∣∣∂x1∂y1∂x2∂y1∂x1∂y2∂x2∂y2∣∣∣∣=∣∣∣∣−y1⋅e−12(y21+y22)−y22π(y21+y22)−y2⋅e−12(y21+y22)y12π(y21+y22)∣∣∣∣=e−12(y21+y22)[−y212π(y21+y22)−y222π(y21+y22)]=−

漫談正態分佈的生成

numpy 生成正態分佈

Python numpy生成正態分佈的平均數

Python隨機生成多維正態分佈

漫談正態分佈的生成

numpy 生成正態分佈

Python numpy生成正態分佈的平均數

Python隨機生成多維正態分佈

相關推薦