機率研究 (4) - 機率分佈生成之 - 倒推分佈生成以及演算法

機率生成函數，就是這幾篇研究筆記中的重點，我們不是只能用均勻分佈去控制偽公平的骰子，於是找到了幾種方法，可以用來生成特定機率分佈函數的機率。

任意的分佈， 並不只限 Binomial, Exponential, Poisson, Chi-Square ，你也可以自訂任意的機率分佈，如:

$$ x \in [0, \frac{3}{2}] ,\space \space f(x) = \frac{1}{4}x^4 $$

手取隨機數推演起手式: 線性變換

在這裡，我們還沒提及連續以及離散型隨機變數，但可以先以連續的情況假設整個分佈的函數。

"幾乎所有程式語言預設的 random 都是均勻分佈" ，這很好見到，畢竟像 Javascript 常見預設的 random 方式就是 [1]:

Math.floor(Math.random() * 10);

這是一個用來產生 0 ~ 9 的正整數的方法。

為什麼要這樣寫，不難想像，就是從 (0, 1) 開區間的均勻隨機中，把整個空間放大 10 倍，也就是整個空間變成 (0, 10) ，不過 x 不會到 10 ，因為原本的機率分佈也都是 (0, 1) 開區間。

$$ \forall x \in (0, 10), x \ne 10 $$

所以當取整後，只會出現 0 ~ 9 的元素，本身的空間是 [0, 9] 10 個正整數。

對分佈的性質來說，又假設取到整個分佈的情形都是累積而來的，可以看到下面的長條圖 (Histogram)，可以想像成是整個機率分佈的積分。

原本 (0, 1) 均勻分布的積分函數可積，做完線性變換之後也是可積分。

附帶一提，對於離散型的隨機變數的分佈，會使用長條圖來表示，對於連續型的變數分佈，就會使用一般的 Plot 來呈現:

如果把離散的分佈用 Plot 去表示，會呈現折線圖的形式出現，這感覺很不好，會產生一些誤會，因為對離散的情形來講， Plot 出分佈出現的折線圖會令人以為在看趨勢，但事實不是這樣。

手取隨機數: 假設解並推導

在這個環節的研究，我受到了這篇文章的啟發:  《如何生成随机数（上）》 作者: pluskid ，但又覺得推導公式上我想寫詳細一些，對於這篇文章做了推導公式的補充。

*原文把區間限制在 [0, 3]，不過對程式語言來說，不是 [0, 1] 閉區間，所以以下呈現都是用 (0, 1) 替換。

從上一節的說明中，我們現在要做的事大概是像這樣:

我們想造一個介於 (0, 3) 區間的一個連續變數，並且服從某個機率密度函數，在這裡我想直接引用原作者本文中使用的函數做示範，然後我們也做一個自己的試試看。

給定我們想手動干涉產生的機率密度函數是: 

$$  p(x) = \frac{1}{9} x^{2} $$

對 x 的範圍，是我們要產生限制區間的重點，這會隨著你所要的區間而變動:

$$ Fixed(固定) X_0 \in (0, 1) \space, Given(給定) X_1 = 3X_0 $$

這樣 X1 就可以限制在 (0, 3) 的區間，而且也是均勻分佈的 (在實數線上 1-1, onto)。

原文解釋 X1 集合中每一個數字的意義，就是拿 X0 = a 的機率相當 X1 = a/3 的機率，因為已經做 3 倍放大，每一個元素就是原來的 1/3 倍。

基於這個解釋，我們才能確立有以下事情發生:

$$ f(x) = \frac{1}{3}, x \in (0, 3) $$

此 f(x) 就是 X1 的機率密度函數 (以此類推， X0 的機率密度函數不就是常數函數 1 嗎)。

然後，我們要通用性的假設:

$$ Suppose:  Y = \phi (X_1), \newline \phi: Strictly \space Monotonic \space Increasing, \newline Continuous \space on \space (0, 3) $$

基於兩條性質，得反函數存在嚴格遞增、連續；定理說明可見 [3] 反函数连续性 - 第一條。

做了這個假設，本文也有反應這是連續隨機變數的假設，在連續的情形下，等於某個值的機率這個說法是不合理的，但是連續函數可以相容變數為離散的情形。

然後，想要把這個 f 函數處理成一個機率分佈函數，然後去相容帶進來的 x (離散隨機變數)，我們還是可以繼續做計算。

上一節有說過，對於分佈的看法，可以用積分來下手，這就是接 下來的思路，想要求得機率函數的積分值，得到分佈:

$$ F(x) = P(X \leq x) = \int_{-\infty}^{x} f(t) dt $$

對剛才的 f(x) = 1/3 做積分，得到的是就是累積分佈函數 (CDF) ，就是 P( X <= x) ，透過 Wikipedia [4] 中，就可以知道他被定義成 F_x(x) = P(X <= x) ，這個 F 在這邊是單變數，則不需要特別註明是誰的函數。 (相關性質請見 [4] )

直接計算上式，就知道:

$$ F(x) = \frac{1}{3} x  $$

對於有界區間的嚴格遞增的連續函數 ϕ 來說，可以對 Y <= y* 套上合成函數，也能對不等式成立。

$$ Y \leq y^{*} $$

等價於:

$$ \boxed{X_1 = \phi^{-1}({Y})} \leq \phi^{-1}({y^{*}}) $$

對於機率累積函數而言，也有相同的等價性質:

$$ P(Y \leq y^{*}) = P(X_1 \leq \phi^{-1}({Y})) $$

基於這個討論下，就有 (最左 G(y) 和最右式 F 的延伸):

$$ \boxed{G(y) = P(Y \leq y^{*})} = \boxed{P(X_1 \leq \phi^{-1}({y})) = F(\phi^{-1}({y}))} $$

然後，在最前面討論到的 CDF 是 PDF(機率密度函數) 的積分，所以把 CDF 微分回去就變成 PDF，此時， G(y) 是我們期望的手動干涉產生的機率累積函數 (CDF)，不過這並不是我們原始表示的樣子，我們原始表達的樣子是 PDF 型式，所以，對上式做微分，得到:

$$ \boxed{g(y) = \frac{dG(y)}{dy}} = f(\phi^{-1}({y})) \frac{d}{dy} \phi^{-1}({y}) \tag{1} $$

從這個式子當作  Baseline，對於最初我們期望的手動干涉產生的機率密度函數，就符合 g(y) 設定的必要條件 (在所設區間連續、嚴格遞增)，就可以令它是 g:

$$ g(y) =  \frac{1}{9} x^{2} $$

$$  \boxed{\frac{1}{9} y^{2} = g(y)} = \boxed{f(\phi^{-1}({y})) \frac{d}{dy} \phi^{-1}({y}) =  \frac{1}{3} \frac{d}{dy} \phi^{-1}({y})} $$

則對兩邊積分消除微分算子 d/dy (分離係數法到左右兩邊，積分)，我們就有這個結果:

$$ \phi^{-1}({y}) = \frac{1}{9}y^3 $$

然後，對 ϕ 求反函數，得到:

$$ \phi({y}) = (9y)^{\frac{1}{3}} $$

將 X1 帶入 ϕ(X1) ，整個區間就都是我們期望的結果。

import numpy as np
import matplotlib.pyplot as plot
 
N = 100000
X0 = np.random.rand(N)
X1 = 3*X0
Y  = np.power(9*X1, 1.0/3)
 
t = np.arange(0.0, 3.0, 0.01)
y = t*t/9


plot.plot(t, y, 'r-', linewidth=1)
plot.hist(Y, bins=50, density = True, facecolor='green', alpha=0.75, histtype='bar', ec='black')
plot.show()

自製實驗

最後，我們可以拿任意一個函數出來實驗:

$$ x \in [0, 5] ,\space \space f(x) = \frac{1}{4}x^2 $$

對每個 X1 的線性變換而言，都是 X0 的 1/4 倍:

$$ X_1 = 5X_0 $$

對 X1 的機率密度函數就是 f(x) = 1/5，求 X1 的機率分佈函數就是積分: F(x) = (1/5)x。

令

$$ g(x) = \frac{1}{4}x^2 $$

透過 (1) 的 Formula ，我們有:

$$  \boxed{\frac{1}{4}y^2 = g(y)} = \boxed{f(\phi^{-1}({y})) \frac{d}{dy} \phi^{-1}({y}) =  \frac{1}{5} \frac{d}{dy} \phi^{-1}({y})} $$

解得:

$$ \phi^{-1}({y}) = \frac{5y^3}{12} $$

求回反函數，得到:

$$ \phi({y}) = \frac{12}{5} x^{\frac{1}{3}}$$

import numpy as np
import matplotlib.pyplot as plot
 
N = 1000
X0 = np.random.rand(N)
X1 = 5*X0
Y  =  np.power((12/5)*X1, 1.0/3)
 
t = np.arange(0.0, 5.0, 0.01)
y = t*t/4


plot.plot(t, y, 'r-', linewidth=1)
plot.hist(Y, bins=50, density = True, facecolor='green', alpha=0.75, histtype='bar', ec='black')
plot.show()

從這張圖可以觀測到，這個積分點只到 y = 1 的時候就沒有繼續增長了，實際機率沒有到 3 的區間這麼多，從任意一個很誇張的分佈函數:

$$ \frac{1}{4}x^4, \forall x \in (0, \frac{3}{2}) $$

所得到的函數也是如此，因為斜率太高，導致太快衝超過 1 ，無法服從分佈，如果只取 [0, 3] 區間中的隨機數，就得不到你要的結果，因為這個分佈的隨機數全部落在大於 4 的區域上頭。

對二次情形來講，這個算法的變參數只有你的 g(x) 的積分，再除右式的 f(x) 機率密度常數函數，然後求反函數，馬上就可以求解。

例如剛才的範例，我就直接先求積分:

$$ \int \frac{1}{4} x^2 = \frac{x^3}{12} $$

然後，剛才的右式的 f 其實是跟區間有關，區間上的機率密度函數是 1/3 就是直接乘上 3 倍給 g(x) 的積分值:

$$  3 \times \frac{x^3}{12} $$ 

乘完後直接求反函數 (懶得算可以問 Wolframe Alpha): inverse of 3 * (x^3) / 12 ，不跟 12 約分其實是保留積分的倍數， 3 其實可以隨著區間去變動，計算也比較方便，解得:

$$ (-2)^{2/3} x^{1/3} $$

wolframe alpha 給出來的解析解有點暴力，手動把它直接變成統一的形式就好了:

$$  (4x)^{1/3} $$

連續型分佈適用的隨機生成演算法

連續型的隨機變數，特性是能取到不相同的值是有限個，或可數無限個。

1. 逆變換法 (Inverse Transform Method, 簡稱: ITM)

這就是上一節手動推演的作法，在《信号处理——生成给定分布随机数》這一篇文章中，也提到「如果我们可以给出概率分布的累积分布函数（F）及其逆函数的解析表达式，则可以非常简单便捷的生成指定分布随机数。」，這不難想像最早期的文章研究中 [5]，看到了從 CDF 變換成機率產生器的一種用法。

2. 取捨法 (Acceptance-Rejection Method, 簡稱: ARM)

ITM 算法其實是需要找出 CDF 的反函數之解析解，要對所有函數做到這點其實有點難，所以 ARM 算法中，只要給出機率密度函數 PDF 的解析表達式就可以產生。

3. 組合法

如果期望的干涉分佈，是可以用各種分佈，經過四則運算(線性變換)以及線性組合來表示時，就可以用組合法去實現，可用的分佈例子有 Normal Distri(Box Muller Method)、 Rayleigh Distri、Poisson Distri。

2, 3 的方法中，可以參考 《Coding Labs - 生成特定分布随机数的方法》這篇文章也給了一個很好的解釋。 [6][7]

離散型分佈適用的隨機生成演算法

連續型的隨機變數，特性是能取到的值可以是連續的，在某個區間可以取任意的一個實數元素。

對於離散分佈適用的隨機數生成演算法，就是第二篇研究所用到的內容，它包含在均勻分布隨機數生成中，也就是說連續型分佈適用的隨機生成，就是非均勻隨機數生成。

回顧之前提到的方法以及新的方法: [14]

1. 線性同餘 LCG

2. FSR 產生器

3. 組合產生器

4. 蒙地卡羅法 [13]

5. Fibonacci

6. RC4

7. 雙橢圓曲線法 Dual Elliptic Curve Deterministic Random Bit Generator

8. 梅森旋轉算法[10]

9. Cryptographically secure pseudo-random number generator

10. C++ Random Device [14]

[11][12]

統計、分佈、機率系統性的文章，可以參考 [8][9] 的文獻繼續研究。

Reference:

[1] https://www.w3schools.com/js/js_random.asp

[2] https://blog.pluskid.org/?p=430

[3] https://baike.baidu.com/item/%E8%BF%9E%E7%BB%AD%E5%87%BD%E6%95%B0

[4] https://zh.wikipedia.org/wiki/%E7%B4%AF%E7%A7%AF%E5%88%86%E5%B8%83%E5%87%BD%E6%95%B0

[5] https://alpynepyano.github.io/healthyNumerics/posts/sampling_arbitrary_distributions_with_python.html

[6] https://blog.codinglabs.org/articles/methods-for-generating-random-number-distributions.html

[7] https://www.cnblogs.com/xingshansi/p/6539319.html

[8 - 統計概觀] http://support.ptc.com/help/mathcad/zh_CN/index.html#page/PTC_Mathcad_Help%2Fabout_statistics.html%23

[9] https://reference.wolframcloud.com/language/tutorial/RandomNumberGeneration.html.zh?source=footer

[10] https://zh.wikipedia.org/wiki/%E6%A2%85%E6%A3%AE%E6%97%8B%E8%BD%AC%E7%AE%97%E6%B3%95

[11] https://zh.wikipedia.org/wiki/%E4%BC%AA%E9%9A%8F%E6%9C%BA%E6%95%B0%E7%94%9F%E6%88%90%E5%99%A8

[12] https://zh.wikipedia.org/wiki/Category:%E4%BC%AA%E9%9A%8F%E6%9C%BA%E6%95%B0%E7%94%9F%E6%88%90%E5%99%A8

[13] https://medium.com/%E6%95%B8%E5%AD%B8-%E4%BA%BA%E5%B7%A5%E6%99%BA%E6%85%A7%E8%88%87%E8%9F%92%E8%9B%87/%E6%8A%BD%E6%A8%A3%E8%88%87%E8%92%99%E5%9C%B0%E5%8D%A1%E7%BE%85-%E4%B8%80-%E9%9A%A8%E6%A9%9F%E6%95%B8%E7%94%9F%E6%88%90-4c4951fc5523

[14] https://www.math.pku.edu.cn/teachers/lidf/docs/statcomp/html/_statcompbook/rng-nonuni.html

[15] https://blog.gtwang.org/programming/cpp-random-number-generator-and-probability-distribution-tutorial/