作者：桂。
時間：2017-03-12 19:31:55
鏈接：http://www.cnblogs.com/xingshansi/p/6539319.html
前言
本文是曲線擬合與分布擬合一文的插曲，進行分布擬合時，碰到一個問題是，如何指定分布的隨機數呢？本文主要包括：
1）連續型隨機數；
2）離散型隨機數；
本文內容為自己的學習筆記，內容多有借鑒他人，在最后一並給出鏈接。
一、連續型隨機數
假設已經擁有U（0，1）的均勻分布數據。
A-逆變換法（Inverse Transform Method）
如果我們可以給出概率分布的累積分布函數（F）及其逆函數的解析表達式，則可以非常簡單便捷的生成指定分布隨機數。
性質：
U是服從[0，1]均勻分布的隨機變量，如果$X = F^{-1}(U)$，則X的分布函數與F相同，即$F_X(x) = F$.
證明：
$F_X(a) = P(X \le a) = P(F^{-1}(U) \le a) = P(U \le F(a)) = F(a)$.即$F_X$與F相同，得證。
算法步驟：
- 生成一個服從均勻分布的隨機數U~Unit(0,1);
- 設F(x)為指定分布的分布函數，則$X = F^{-1}（U）$即為指定分布的隨機數。
示例：生成滿足$\lambda = 2$的指數分布隨機數。
分析：
由$f(x)$得出$F(x)$ —> $F(x) = 1 - {e^{ - \lambda x}}$，進而求得$F(x)$逆函數，得出$X = {F^{ - 1}}(u) = - \frac{1}{\lambda }\ln (1 - u)$.
代碼：
Len = 1000000; u = rand(1,Len); lemda = 2; x = -1/lemda*(log(1-u));
對應結果：
B-舍選法（Acceptance-Rejection Method）
一般來說逆變換法是一種很好的算法，簡單且高效，如果可以使用的話，是第一選擇。但是逆變換法有自身的局限性，就是要求必須能給出分布函數F逆函數的解析表達式，有些時候要做到這點比較困難，這限制了逆變換法的適用范圍。
當無法給出分布函數F逆函數的解析表達式時，舍選法是另外的選擇。舍選法的適用范圍比逆變換法要大，只要給出概率密度函數的解析表達式即可，而大多數常用分布的概率密度函數是可以查到的。
算法思想：
給定輪廓，並在輪廓范圍內生成均勻分布序列；
算法實現：
- 設概率密度函數為$f(x)$，首先生成一個均勻分布隨機數X~Unit($x_{min}$，$x_{max}$);
- 獨立地生成一個均勻分布隨機數Y~Unit($y_{min}$，$y_{max}$);
- 若$Y \le f(X)$，則返回X，否則重復第一步。
給出一張舍選法的原理圖：
代碼（data即為最終的隨機數）：
xmin = 0; xmax = 5; Len = 1000000; x = (xmax-xmin)*rand(1,Len)-xmin; lemda = 2; fx = lemda*exp(-lemda*x); ymax = lemda; ymin = 0; Y = (ymax-ymin)*rand(1,Len)-ymin; data = x(Y<=fx);
結果圖：
舍選法需要對數據重復篩選，性能不如逆變換法，但其適應性更廣，無法得到分布函數的逆函數時，舍選法不失為一個選擇。
其實本質就是取分布以內的隨機數，顛倒過來，給定分布，將分布無限細分，生成各區間對應數量的隨機數，也可以實現近似。
C-組合法
當目標分布可以用其它分布經過四則運算表示時，可以使用組合算法生成對應隨機數。此部分僅以幾個例子簡要介紹。
例一：正態分布（Box Muller方法）
正態分布隨機數的產生，除了上文的方式，經典的就是Box Muller方法，所有正態分布均可由標准正態分布演變得出。
Box Muller理論推導：
設$(X,Y)$是一對相互獨立的服從正態分布$N(0,1)$的隨機變量，則有概率密度函數（指數少一負號）：
轉化為極坐標形式，令
，則R有分布：
Eq.(1)
令
，則分布函數的反函數為：
由逆變換法性質可知：
滿足$F_R$的分布可由$1-Z$~$U(0,1)$得出，亦即：可由$Z$~$U(0,1)$得出；
再分析$P_{\theta}$，同樣對Eq.(1)中表達式關於$r$進行從0到∞的積分，容易得出$\theta$服從均勻分布，因此該隨機數可由均勻分布直接產生；
又易證$P(r;\theta)=P(r)P(\theta)$，即二者獨立，因此二者可由不同的均勻分布分別生成。
以上為Box Muller的原理分析。
Box Muller算法實現：
- 分別生成兩組均勻分布隨機數：$U_1$~$U(0,1)$、$U_2$~$U(0,1)$；
- 生成$R$以及$\theta$的分布：
- 生成兩組獨立的標准正態分布:
- 生成符合給定均值、方差具體要求的正態分布；
代碼(生成標准正態分布)：
Len = 10000; U1 = rand(1,Len); U2 = rand(1,Len); R = sqrt(-2*log(U1)); theta = 2*pi*U2; X = R.*cos(theta); Y = R.*sin(theta);
效果圖：
啰嗦一句：從scatter(X,Y)的結果中，可以直觀看出X、Y不相關，但不能確定不獨立。$X*Y = (X,Y)$與$P_X*P_Y =P (X,Y)$不是一個概念。對於同樣一個圓形區域，X、Y可以按不同的密度分布落入。
再看看均勻分布的scatter(X,Y)：
scatter結果：一個圓形、一個方形，但兩種情況下X、Y都是獨立同分布的。
例二：瑞利分布（Rayleigh Distribution）
其實Box Muller方法的中間過程，包含了瑞利分布，即$P(r)$。給出PDF定義式：
$f(r) = \frac{r}{{{\sigma ^2}}}{e^{ - \frac{{{r^2}}}{{2{\sigma ^2}}}}}$
其中$\sigma > 0$.
理論分析：
設$(X,Y)$是一對相互獨立的服從正態分布$N(0,\sigma ^2)$的隨機變量，則有概率密度函數：
$f\left( {x,y} \right) = f(x)f(y) = \frac{1}{{2\pi {\sigma ^2}}}{e^{ - \frac{{{x^2} + {y^2}}}{{2{\sigma ^2}}}}}$
轉化為極坐標，$dxdy = rdrd\theta $，並對$\theta$求取積分，得：
$f\left( r \right) = \frac{r}{{{\sigma^2}}}{e^{ - \frac{{{r^2}}}{{2{\sigma ^2}}}}}$
可見：設$(X,Y)$是一對相互獨立的服從正態分布$N(0,\sigma ^2)$的隨機變量，按Box Muller的方法構造R，即為對應的瑞利分布。
回顧公式：$R = \sqrt {{X^2} + {Y^2}} $，兩個獨立同正態分布（零均值）信號的包絡是瑞利分布。
算法步驟：
- 生成一組均勻分布隨機數：$U$~$U(0,1)$；
- 根據給定的參數$\sigma$，生成：$R = \sigma \sqrt { - 2\ln U} $;R即為滿足要求的瑞利分布；
對應代碼：
Len = 1000000; U1 = rand(1,Len); U2 = rand(1,Len); sigma = 2; R = sigma*sqrt(-2*log(U1));
結果圖：
例三：泊松分布（Poisson distribution）
- 背景介紹
一段時間內，事件發生的次數，組成計數過程。泊松分布是一種常用的計數過程。
生活中，大量事件是有固定頻率的，即可以通過一段時間內計數找出規律：
- 某醫院平均每小時出生3個嬰兒
- 某公司平均每10分鍾接到1個電話
- 某超市平均每天銷售4包xx牌奶粉
- 某網站平均每分鍾有2次訪問
這里不展開講泊松分布的來龍去脈，直接給出公式：
再來回顧指數分布，指數分布描述的是：計數過程中，相鄰事件發生的時間間隔概率分布。
對應上面的問題，則有
- 嬰兒出生的時間間隔
- 來電的時間間隔
- 奶粉銷售的時間間隔
- 網站訪問的時間間隔
指數分布可由泊松分布推出：
泊松分布 —> 指數分布：
設相鄰兩次事件間隔為$T$，起始時刻為$T_{start}$，則終止時刻為$T_{start}+T$，$P\{ T \ge t \}$表示$[T_{start},T_{start}+t]$時間內沒有事件發生，即：
從而事件發生的概率為：
對應概率密度為：
對於泊松分布，得出結論：
- 兩次事件的時間間隔為負指數分布，且均值為$\frac{1}{\lambda }$;
- 事件1與事件2，事件2與事件3，...事件n與事件n+1，相鄰事件的時間間隔具有獨立同分布的特性。
- 算法實現
此處給出的泊松分布隨機數，產生原理基於指數分布，更多方法可以戳這里。
首先給出算法思想與算法步驟：
算法思想：
根據上文分析可知，Poisson分布對應一段時間內（記為$t_{max}$）事件發生次數，而指數分布exp對應相鄰時間的時間分布;
反過來，已知兩兩相鄰的時間變量，該變量服從指數分布且相互獨立，所有時間變量相加—>並讓其不超過一段時間的總量$t_{max}$，則累加數的分布對應泊松（Poisson）分布。
算法步驟：
- 生成一組均勻分布隨機數：$U$~$U(0,1)$;
- 利用逆變換法生成一系列獨立的指數分布$X_i$ ~ $exp(\lambda )$;
- 記
如果$Y > t_{max}$，則停止，並輸出$k-1$;否則，繼續生成$X_{k+1}$，直到$Y > t_{max}$為止；
- 循環操作過程3；
輸出的一系列整數（值為$k-1$）服從參數為$\mu = \lambda t_{max}$的Poisson分布。
由於上文已經交代如何生成指數分布，代碼中直接調用指數隨機量exprnd函數，不再重復生成；本文主要講究思路，事實上，本文列舉的例子都有現成的函數調用。
給出代碼（data即為生成的數據）：
Iter = 10000; % 迭代次數，即data個數
lambda = 2;
t_max = 11;
%phei = lambda*t_max;%均值
for m=1:Iter,
k=1;
Y(m,1)=exprnd(1/lambda);
while Y(m,k) <= t_max % 時間間隔 (0,t_max]
Y(m,k+1)=Y(m,k)+exprnd(1/lambda);
k=k+1;
end
data(m)=k-1;
end
測試結果(將代碼生成結果與MATLAB自帶函數poissrnd.m的結果對比)：
二、離散型隨機數
假設已經擁有U（0，1）的均勻分布數據，離散型隨機數根據[0,1]划分區間，給出對應變量值即可。
例：離散隨機變量X的分布率如下表，試生成該隨機數。
| X | -1 | 1 | 3 |
| P（x） | 0.7 | 0.2 | 0.1 |
tag = rand(1,1000); data = tag; data(tag>0.9) = 3; data(tag<=0.7) = -1; data(abs(data)<1) = 1;
直方圖結果：
可以看出結果服從指定的分布率。
補充：本文僅簡要介紹一維分布隨機數的生成方式，有兩個主要的點沒有涉及：
- 均勻分布隨機數的產生方式；
- 多維隨機數的產生方式。
參考：