機率研究 (3) - 機率分佈生成之 - 應用與分佈的意義

前面第一章，我們在不控制隨機產生器本身的情況下，探討機率相關應用，第二章，我們從機率產生器本身做了一點研究，總結了在第二章中所期望的是「混亂」的機率，談論的是可測與不可測的問題，回到這篇文章，主題是希望按照分佈去生產一個亂數，也就是在「可測」範圍內探討機率分佈的問題。

我們需要搞懂一些東西

對於一般的亂數產生器來說，都是使用均勻分佈的機率 (Uniform) ，我們可以透過第一篇研究的方法，很輕易的決定:

有五個武器 A,B,C,D,E ，這五個武器各自被抽中的機率是多少。

這是一個離散型隨機分布。

回到各種層面的應用上，我們會需要依照各種不同的情境，應用【連續型隨機分布】，

首先，我們得先搞懂要做一個特定分佈的亂數產生器，跟預期的結果是什麼關係。

釐清過程、結果、反推

對於一正一反的硬幣，一次投擲五個，其正反面的個數的機率是多少，這是一個從【結果】反推整個機率系統的過程，也就是統計，但從統計反推回原本的機率，該如何做到，如果硬幣不公平，我們可以怎麼設計。

對於我們需要控制機率，最常見的方式就是從 CDF 機率累積函數反推(機率密度函數的積分就是機率累積函數)，他的演算法叫做 Inverse Transform，是基於把均勻分布的 1x1 正方形，做線性變換映射到 CDF 函數上，這會需要主動拆解 CDF 的反函數，所以要確定函數可逆，也要找到解析解、級數解去把它做映射，得到一個偽隨機數產生器(Random Number Generator: RNG)，且服從你指定的分佈，詳細的數學探討，會在下一篇文章解釋。

import numpy as np
from pylab import *

# Create some test data
dx = 0.01
X  = np.arange(-2, 2, dx)
Y  = exp(-X ** 2)

# Normalize the data to a proper PDF
Y /= (dx * Y).sum()

# Compute the CDF
CY = np.cumsum(Y * dx)

# Plot both
plot(X, Y)
plot(X, CY, 'r--')

常態分佈 RNG Algorithm - Box Muller Transform

在這篇文章中，會稍微詳細解釋常態分佈的常見隨機數產生器做法，用於把常態分佈取線下方到 x 水平軸所有的均勻分佈機率生成，常見有好幾種演算法:

1. Inverse Transform
2. Rejection Sampling
3. Importance Sampling
4. Sampling-Importance-Resampling
5. Markov Chain Monte Carlo

...

在下一篇文章會詳細說明到 Inverse Transform 的推演，在處理【常態分佈】情況下，會看到有兩大演算法:

1. Box Muller
2. Ziggurat [14][15][16]
3. Marsaglia polar method (避免 Box-Muller 做三角函數運算帶來的系統負載的演算法)

在這邊，只打算提出 Box Muller 的做法，上一小節，有提到從分佈中反推機率產生器，需要知道 CDF 可逆函數的解析解或是級數解，常見求可逆以數學分析的方法來看，有很重要的理論在支持有界閉區間上的連續函數可逆的性質，那就是: 只要是初等函數，在其定義的區間上都連續，而且都有可逆的性質。 

離題了，回過頭來談談 Box Muller 之前，先講完我們要如何開始從實作面上進行流程，整體而言，我們希望從 Inverse Transform Method 這個方法回推機率產生器。

Inverse Transform Method

流程是，我們需要找到該分佈函數的機率累積分佈函數 (CDF) 的反函數，從 Wikipedia 搜尋 Normal Distribution，你就可以看到常態分布(高斯分佈) 的 CDF，以下參考 Hulu 文章，使用平移伸縮線性變換後到簡化的形式 [17] (注意，若參考 [17] 的做法，也應該注意 [18] 所列出的參數沒有根號 pi，因為已經做過變換簡化了) :

$$ \text{erf}(x) =  \frac{2}{\sqrt{\pi}} \int_{0}^{x} e^{-t^2} dt $$

erf 就是高斯函數的機率累積分佈函數，又是一種誤差函數，不過在這裡，我們如果要使用 Inverse Transform Method，勢必要找到他的反函數，不過，erf 並不是初等函數構成的函數，可能沒有解析解，這會導致很難辦到這件事，在 How to generate Gaussian samples [18] 這篇文章一開始有提到可以用泰勒級數 (c=0 Maclaurin series 特例) 來逼近函數，但顯然我們要面對的問題是要對兩個點做更多的逼近採樣，會造成很大的麻煩，因此， Box-Muller 提出了針對 Inverse Trasnsform 處理的方案。

Box-Muller Transform

把兩個獨立的常態分佈的變數，對分佈做聯合，推演出極座標版的機率累積分佈函數，然後座標進行變換，用極座標來處理樣本，得到服從常態分佈的 (x, y)，假設 x, y 服從常態分佈獨立隨機變數，則聯合機率密度為:

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

這個 (A1) 是怎麼來的? 證明 [27]:

我們得先回歸原點，因為命題需要某些形式帶入，所以我們得要證明出:

$$ \int^{\infty}_{-\infty} e^{\frac{-x^2}{2}} dx = 2\pi $$

這個東西是來自高斯積分的參數變形，如果你對數學分析有一些概念，你可以馬上就計算出答案，不過因為這行式子我們得把 x 變數通用到所有各種情形推廣，做以下證明:

$$ \text{令    }  I =  \int^{\infty}_{-\infty} e^{\frac{-x^2}{2}} dx $$

則

$$ I^2 =  \int^{\infty}_{-\infty} e^{\frac{-x^2}{2}} dx\int^{\infty}_{-\infty} e^{\frac{-y^2}{2}} dy =  \int^{\infty}_{-\infty} e^{-\frac{x^2 + y^2}{2}} dxdy  \tag{式0}$$

然後，把用座標變換，轉換成 Polar Form 的形式處理:

$$ \text{令    } x=r cos \theta \text{,  } y = r sin \theta $$

其中，(式 0) 透過多重積分 Jacobian 變數變換，所以建立 r, θ 跟 x^2 + y^2 的關係式:

$$ x^2 + y^2 = r^2 $$

$$ tan^{-1} \frac{y}{x} = \theta $$

還有 dxdy 用 Jacobian 轉換的 r 關係式:

$$ \frac{ \partial(x,y) }{ \partial(r, \theta) } = \begin{vmatrix} \frac{\partial x}{\partial r} & \frac{\partial x}{\partial \theta } \\ \frac{\partial y}{\partial r} & \frac{\partial y}{\partial \theta} \end{vmatrix}  =  \begin{vmatrix} cos \theta & -r sin \theta \\ sin \theta & r cos \theta \end{vmatrix} = r cos^2 \theta + r sin^2 \theta = r $$

則有:

$$ I^2 = \int^{2\pi}_{0} \int^{0}_{\infty} e^{\frac{-r^2}{2}} r dr d\theta = 2\pi \int^{\infty}_{0} e^{\frac{-r^2}{2}} r dr = 1 \times 2\pi $$

為什麼:

$$ \int^{\infty}_{0} e^{\frac{-r^2}{2}} r dr = 1 $$

證明:

$$ \int^{\infty}_{0} e^{\frac{-r^2}{2}} r dr \\ \\ = \frac{1}{2} \int^{\infty}_{0} e^{\frac{-r^2}{2}} dr^2  \\ \\  = -\int^{\infty}_{0} e^{\frac{-r^2}{2}} d(- \frac{r^2}{2}) \\ \\  =  -e^{\frac{-r^2}{2}} {\Bigm\vert}_0^\infty \\ = -(0-1) = 1 \tag{式 1} $$

以上(式 1)補充證明，感謝彰師大李國瑋教授補充證明高斯積分部分得到 1 的結果。

這是一個很哲學的問題，當你問出: 「為什麼要證明這個?」 回答是: 「因為這樣才做得出來」，不過很多都是在推演過程中碰到問題才回頭解決，因為不能保證所有人都看得懂，所以要做出前情提要，以上的證明就是一種前情提要。

回歸正題，因為 Box-Muller 是把兩個獨立的常態分佈的聯合分佈，所以聯合後才會得到 (A1)，也就是兩個函數的相乘，注意把 2π 移到左式，就會有 1/2π 這個東西，我們把它定義成以下的形式 (f 與 x,y 有關)。

$$ f_{(x, y)} (x,y) = \frac{1}{2\pi} e^{-\frac{x^2 + y ^2}{2}} \tag{A2} $$

然後，開始做極座標轉換，先對參數做限制:

$$ \text{令    } x = R cos \theta \text{,    } y = R sin \theta \text{, 其中 } \theta \in [0,2\pi]  $$

則，R 有以下這個分佈函數 (有 u 這個變數變換):

$$ P(R \leq r) =  \int^{2\pi}_{0}\int^{0}_{x} \frac{1}{2\pi} e^{-\frac{u^2}{2}}u du d\theta =  \int^{0}_{x} e^{-\frac{u^2}{2}}u du = 1 - e^{-\frac{x^2}{2}}\tag{B1} $$

以上的做法不只有一種，也有使用 Jacobian 處理的證明，可以見 [28][29][17][18][26][21][22]。

現在，在這裡可以把 (B1) 看成是極座標 r (radius) 的累積分佈函數， (B1) 的反函數也可以求了，所以繼續做以下推導:

$$ \text{令    } F_R (r) = 1 - e^{-\frac{x^2}{2}} $$

$$ \text{令    } R = {F_R}^{-1} (Z) = \sqrt{-2 ln(1-Z)} $$

所以，如果 x 本身是均勻分佈 (Uniform) 隨機產生器出來的均勻變數，則 R 的分佈函數就是 FR(r)，也就是 R 這個值去計算 CDF 的反函數生成的分佈函數。

其中 1-Z 是常見的形式 1-p ，詳情也可以回顧 (odds ratio: p/(1-p))，可以把 uniform random 產生的變數等於 1-Z 去處理，即:

$$ \text{令    } u_1 = (1-Z), \text{令    }\theta = 2\pi u_2 $$

所以:

$$ x  = R cos \theta = \sqrt{-2 ln(u_1)} cos(2\pi u_2) \text{,    } y = R sin \theta = \sqrt{-2 ln(u_1)} sin(2\pi u_2) $$

總而言之，就是你的 u_1, u_2 是從 math.random 這個函數給定的，而 x, y 則是兩種解，你可以用 x 或 y ，任一種當作產生常態分佈的隨機變數。

程式碼，兩行解決:

u1 = random()
print(sqrt(-2*log(u1))*cos(2*pi*u1)) # pi 是預設就有的

破除迷思

上一節提到了眾多可以從分佈函數中逆推機率產生器的方式，顯然都不簡單，要求一個【簡單的廣義分佈隨機數生成作法】並不存在，每個分佈都有它的挑戰。

在以上的例子及推演，展示手作分佈的難點，就是在生成分佈的演算法這件事，要操作機率，你需要了解的事實上比操作 uniform 來得更複雜一些。

常態分佈，以及其他的分佈

常態分佈

from numpy import random
import matplotlib.pyplot as plt
import seaborn as sns

sns.distplot(random.normal(size=1000), hist=False)

plt.show()

Bernoulli / Binomial 機率分佈生成 (Matlab)

第一章，我們見過 Bernoulli 的 PMF 樣子，現在，希望從生成器中產生一個相似的機率分佈:

N = 10:10:100

r1 = binornd( N, 1./N )

從 10 ~ 100，每次增長 10 的陣列做為事件個數 N。

輸出後:

N =

    10    20    30    40    50    60    70    80    90   100

r1 =

     0     3     0     0     1     1     0     1     0     3

要觀察結果，可以從 sort 後看到圖形:

plot(sort( r1 ))

Bernoulli / Binomial 分佈生成 (Python)

from numpy import random
import matplotlib.pyplot as plt
import seaborn as sns

sns.distplot(random.normal(loc=50, scale=5, size=1000), hist=False, label='normal')
sns.distplot(random.binomial(n=100, p=0.5, size=1000), hist=False, label='binomial')

plt.show()

Poisson 分佈生成 (Python)

from numpy import random
import matplotlib.pyplot as plt
import seaborn as sns

sns.distplot(random.normal(loc=50, scale=7, size=1000), hist=False, label='normal')
sns.distplot(random.poisson(lam=50, size=1000), hist=False, label='poisson')

plt.show()

均勻分佈生成 (Python)

from numpy import random
import matplotlib.pyplot as plt
import seaborn as sns

sns.distplot(random.uniform(size=1000), hist=False)

plt.show()

Logistic 分佈生成 (Python)

from numpy import random
import matplotlib.pyplot as plt
import seaborn as sns

sns.distplot(random.normal(scale=2, size=1000), hist=False, label='normal')
sns.distplot(random.logistic(size=1000), hist=False, label='logistic')

plt.show()

多項式分佈生成 (Python)

from numpy import random
import matplotlib.pyplot as plt
import seaborn as sns

sns.distplot(random.multinomial(n=6, pvals=[1/6, 1/6, 1/6, 1/6, 1/6, 1/6]), hist=False, label='multinomial')

plt.show()

指數分佈生成 (Python)

from numpy import random
import matplotlib.pyplot as plt
import seaborn as sns

sns.distplot(random.exponential(size=1000), hist=False)

plt.show()

卡方分佈生成 (Python)

from numpy import random
import matplotlib.pyplot as plt
import seaborn as sns

sns.distplot(random.chisquare(df=1, size=1000), hist=False)

plt.show()

Rayleigh 分佈生成 (Python)

from numpy import random import matplotlib.pyplot as plt import seaborn as sns sns.distplot(random.rayleigh(size=1000), hist=False) plt.show()

Pareto 分佈生成 (Python)

from numpy import random import matplotlib.pyplot as plt import seaborn as sns sns.distplot(random.pareto(a=2, size=1000), kde=False) plt.show()

Zipf 分佈生成 (Python)

from numpy import random import matplotlib.pyplot as plt import seaborn as sns x = random.zipf(a=2, size=1000) sns.distplot(x[x<10], kde=False) plt.show()

幾何分佈生成 (Python)

from numpy import random import matplotlib.pyplot as plt import seaborn as sns sns.distplot(random.geometric(p=0.35, size=10000), hist=False) plt.show()

Gamma 分佈生成 (Python)

from numpy import random import matplotlib.pyplot as plt import seaborn as sns sns.distplot(random.gamma(2., 2., 1000), hist=False) plt.show()

分佈產生關係圖:

Reference:

[1] https://blog.csdn.net/GregoryHanson/article/details/77823081

[2] https://blog.csdn.net/ljyljyok/article/details/86657942

[3] https://cn.comsol.com/blogs/sampling-random-numbers-from-probability-distribution-functions/

[4] https://www.itl.nist.gov/div898/handbook/eda/section3/eda3667.htm

[5] https://ivanky.wordpress.com/2014/05/23/random-numbers-in-matlab/

[6] https://stackoverflow.com/questions/26099170/bernuli-geometric-simulation-on-matlab

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

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

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

[10] https://www.thinbug.com/q/25471457

[11] https://blog.csdn.net/fengying2016/article/details/80593266

[12] https://blog.csdn.net/wangyangzhizhou/article/details/80088490

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

[14] https://en.wikipedia.org/wiki/Ziggurat_algorithm

[15] https://blog.csdn.net/Real_Myth/article/details/51013680

[16] 隨機測驗資料的產生模式 - 許玟弒、施慶麟 

[17] https://zhuanlan.zhihu.com/p/32578539

[18] https://medium.com/mti-technology/how-to-generate-gaussian-samples-3951f2203ab0

[19] http://4rdp.blogspot.com/2008/06/random-variable-of-normal-distribution.html

[20] 中大: http://www.math.ncu.edu.tw/~yu/sc97/boards/lec11_sc1_97.pdf

[21] http://financelab.nctu.edu.tw/www/course/QF/NormalGen.pdf

[22] https://ichi1234567.github.io/2012/07/21/box-muller/

[23] https://glowingpython.blogspot.com/2013/01/box-muller-transformation.html

[24] https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform

[25] http://www.math.ncu.edu.tw/~yu/sc97/boards/lec11_sc1_97.pdf

[26] https://zhuanlan.zhihu.com/p/38638710

[27] Box-Muller方法的证明

[28] Random variable of normal distribution (作者：Bridan) - 刊載於程式人雜誌 -- 2014 年 9 月號

[29] https://blog.csdn.net/weixin_39910711/article/details/101773776

[30] https://stackoverflow.com/questions/57335218/how-to-generate-random-x-and-y-coordinates-using-box-muller-transformation

[31] https://www.w3schools.com/python/numpy_random.asp

[32] http://www.hmwu.idv.tw/web/R_AI_M/AI-M1-hmwu_R_Stat&Prob_v2.pdf