生成特定分布隨機數的方法

生成隨機數是程序設計里常見的需求。一般的編程語言都會自帶一個隨機數生成函數，用於生成服從均勻分布的隨機數。不過有時需要生成服從其它分布的隨機數，例如高斯分布或指數分布等。有些編程語言已經有比較完善的實現，例如Python的NumPy。這篇文章介紹如何通過均勻分布隨機數生成函數生成符合特定概率分布的隨機數，主要介紹Inverse Ttransform和Acceptance-Rejection兩種基礎算法以及一些相關的衍生方法。下文我們均假設已經擁有一個可以生成0到1之間均勻分布的隨機數生成函數，關於如何生成均勻分布等更底層的隨機數生成理論，請參考其它資料，本文不做討論。

• 基礎算法
  • Inverse Transform Method
    • ITM算法描述
    • ITM算法說明
    • ITM實現示例
  • Acceptance-Rejection Method
    • ARM算法描述
    • ARM算法說明
    • ARM實現示例
• 衍生算法
  • 組合算法
  • 生成具有相關性的隨機數
• 更多參考

基礎算法

Inverse Transform Method

最簡單的生成算法是Inverse Transform Method（下文簡稱ITM）。如果我們可以給出概率分布的累積分布函數（下文簡稱CDF）及其逆函數的解析表達式，則可以非常簡單便捷的生成指定分布隨機數。

ITM算法描述

1. 生成一個服從均勻分布的隨機數U∼Uni(0,1)
2. 設F(X)為指定分布的CDF，F−1(Y)是其逆函數。返回X=F−1(U)作為結果

ITM算法說明

這是一個非常簡潔高效的算法，下面說明其原理及正確性。

我們通過圖示可以更直觀的明白算法的原理。下圖是某概率分布的CDF：

我們從橫軸上標注兩點xa和xb，其CDF值分別為F(xa)和F(xb)。

由於U服從0到1之間的均勻分布，因此對於一次U的取樣，U落入F(xa)和F(xb)之間的概率為:

P(U∈(F(xa),F(xb)))=F(xb)−F(xa)

而由於CDF都是單調非減函數，因此這個概率同時等於X落入xa和xb之間的概率，即：

P(U∈(F(xa),F(xb)))=P(F−1(U)∈(F−1(F(xa)),F−1(F(xb))))=P(X∈(xa,xb))=F(xb)−F(xa)

而根據CDF的定義，這剛好說明X服從以F(x)為CDF的分布，因此我們的生成算法是正確的。

ITM實現示例

下面以指數分布為例說明如何應用ITM。

首先我們需要求解CDF的逆函數。我們知道指數分布的CDF為

F(X)=1−e−λX

通過簡單的代數運算，可以得到其逆函數為

F−1(Y)=−1λln(1−Y)

由於U服從從0到1的均勻分布蘊含着1−U服從同樣的分布，因此在實際實現時可以用Y代替1−Y，得到：

F−1(Y)=−1λln(Y)

下面給出一個Python的實現示例程序。

1. import random
2. import math
3. def exponential_rand(lam):
4. if lam <=0:
5. return-1
6. U = random.uniform(0.0,1.0)
7. return(-1.0/ lam)* math.log(U)

Acceptance-Rejection Method

一般來說ITM是一種很好的算法，簡單且高效，如果可以使用的話，是第一選擇。但是ITM有自身的局限性，就是要求必須能給出CDF逆函數的解析表達式，有些時候要做到這點比較困難，這限制了ITM的適用范圍。

當無法給出CDF逆函數的解析表達式時，Acceptance-Rejection Method（下文簡稱ARM）是另外的選擇。ARM的適用范圍比ITM要大，只要給出概率密度函數（下文簡稱PDF）的解析表達式即可，而大多數常用分布的PDF是可以查到的。

ARM算法描述

1. 設PDF為f(x)。首先生成一個均勻分布隨機數X∼Uni(xmin,xmax)
2. 獨立的生成另一個均勻分布隨機數Y∼Uni(ymin,ymax)
3. 如果Y≤f(X)，則返回X，否則回到第1步

ARM算法說明

通過一幅圖可以清楚的看到ARM的工作原理。

ARM本質上是一種模擬方法，而非直接數學方法。它每次生成新的隨機數后，通過另一個隨機數來保證其被接受概率服從指定的PDF。

顯然ARM從效率上不如ITM，但是其適應性更廣，在無法得到CDF的逆函數時，ARM是不錯的選擇。

ARM實現示例

下面使用ARM實現一個能產生標准正態分布的隨機數生成函數。

首先我們要得到標准正態分布的PDF，其數學表示為：

f(x)=1√2πe−x22

為了方便，這里我會直接使用SciPy來計算其PDF。

程序如下。

1. import random
2. import scipy.stats as ss
3.  
4. def standard_normal_rand():
5. whileTrue:
6. X = random.uniform(-3.0,3.0)
7. Y = random.uniform(0.0,0.5)
8. if Y < ss.norm.pdf(X):
9. return X

注意：標准正態分布的x取值范圍從理論上說是(−∞,∞)，但是當離開均值點很遠后，其概率密度可忽略不計。這里只取(−3.0,3.0)，實際使用時可以根據具體需要擴大這個取值范圍。

衍生算法

組合算法

當目標分布可以用其它分布經過四則運算表示時，可以使用組合算法生成對應隨機數。

最常見的就是某分布可以表示成多個獨立同分布（下文簡稱IID）隨機變量之和。例如二項分布可以表示成多個0-1分布之和，Erlang分布可以由多個IID的指數分布得出。

以Erlang分布為例說明如何生成這類隨機數。

設X1,X2,⋯,Xk為服從0到1均勻分布的IID隨機數，則−1λlnX1,−1λlnX2,⋯,−1λlnXk為服從指數分布的IID隨機數，因此

X=−1λlnX1−1λlnX2−⋯−1λlnXk=−1λlnk∏i=1Xi∼Erl(k,λ)

所以生成Erlang分布隨機數的算法如下：

1. 生成X1,X2,⋯,Xk∼Uni(0,1)
2. 返回−1λln∏ki=1Xi

這類分布的隨機數生成算法很直觀，就是先生成相關的n個IID隨機數，然后帶入簡單求和公式或其它四則公式得出最終隨機數。其數學理論基礎是卷積理論，稍微有些復雜，這里不再討論，有興趣的同學可以查閱相關資料。

生成具有相關性的隨機數

現在考慮生成多維隨機數，以最簡單的二維隨機數為例。

如果兩個維度的隨機數是相互獨立的，那么只要分別生成兩個列就可以了。但是如果要求兩列具有一定的相關系數，則需要做一些特殊處理。

下列算法可以生成兩列具有相關系數ρ的隨機數。

1. 生成IID隨機變量X和Y
2. 計算X′=ρX+√1−ρ2Y
3. 返回(X,X′)

可以這樣驗證其正確性：

再分享一下我老師大神的人工智能教程吧。零基礎！通俗易懂！風趣幽默！還帶黃段子！希望你也加入到我們人工智能的隊伍中來！https://blog.csdn.net/jiangjunshow