偽隨機亂數生成器(Psedorandom number generator,PRNG)

在模擬真實數據時，我們需要一個穩定、可再現性高的演算法來產生隨機數以透過變數變換得到各種不同分配的數據。

偽隨隨機變數產生器，是給定一個初始值 ... 0 偽隨機亂數生成器(Psedorandomnumbergenerator,PRNG) 統計及資料科學 Mar3 WrittenByoswald 亂數(隨機數)(Randomnumber)上篇文章中講到如何用服從均勻分配的隨機變數造出特定分配的隨機變數，那麼服從均勻分配的隨機變數又要怎麼產生呢?首先我們要瞭解亂數(Randomnumber)，這和我們所熟悉的一個概念隨機變數(Randomvariable)好像有些相似，其實亂數就是一種特定的隨機變數，有著特定的性質。

產生一組好的亂數，為了達到模擬真實情況、安全加密等需求，我們會希望它有以下幾種性質:1.服從均勻分配2.統計獨立3.可再現我們可以透過很多無母數方法檢定驗證第1、2點，如我們所熟悉的卡方檢定(chi-squaretest)、K-Stest可以檢定一組亂數對不同分配的適配程度；runtest或透過列連表做卡方檢定可以檢定亂數的產生是否互相獨立。

那麼如何確保第3點呢?這就要從產生亂數時所用的算法下手了。

最早產生亂數的方法大多是擲骰子、撲克牌之類人為影響小，且較有隨機性的隨機實驗取得。

這種方法的缺點顯而易見，就是速度慢，時間一長，亂數表很容易不敷使用。

此外還有另一個缺點，但同時也是優點，就是可再現性低。

為什麼這麼說呢，如果我們的需求是要模擬真實環境數據尋找特定因子或測試實驗方法的容錯能力，那麼我們一定希望當寫出一個能達到自己目的的演算法能保存起來並且每次使用都能達到同樣的效果，甚至同樣的結果(以節省執行時間和佔用的儲存空間)，這時我們就需要一個可再現性高的演算法；而假如我們的需求是透過亂數提升安全性的加密方式，那麼我們就需要可再現低的演算法，或使用自然亂數，以杜絕透過大量結果逆推原始資訊的情形。

註:自然亂數:透過物理現象取得、存在於自然界的真實亂數，如:無線電波的雜訊、某放射性物質的衰變速度$\cdots$，或是一些特定的無理數，如:$\pi、e\cdots$，一些硬體也可以達到近似的效果，一般來說無法逆推或預測。

偽隨機數生成器(Pseudorandomnumbergenerator,PRNG)其實除了自然亂數外，任何一種演算方法所產生的亂數都是偽隨機數(Pseudorandomnumber)，因為這些數字追根究柢都是透過各種函數得出的固定數，即使是擲骰子、輪盤也可透過個人操作習慣、環境等因子進行分析。

然而隨機數的真偽並不代表任何好壞依據，就如上述所說，是使用時機的差別，偽隨機亂數產生器有固定模式可以修改，且在實務上較有效率；自然(真實)隨機數無法逆推和預測的能大幅提升保密性，很多時候為達成需求需要同時使用。

綜上所說，偽隨機數生成器就是一個透過給定初始值，或稱種子(Randomseed)，產生近似真實隨機數數列的算法，以下會介紹幾種運用廣泛或較為基本的演算法。

中間平方法(Middle-squaremethod) Middle-squaremethod,中譯為中間平方法或平方取中法，由馮·諾伊曼(JohnvonNeumann)開發，方法如下: 如果要生成一$n$位數字組成的隨機數列，則將初始的$n$位數平方，取居中的$n$位數字， 再由得到的數字重複動作得到數列。

如$n=4$，初始值$=7777$，依照中間平方法 $$ 7777^{2}$=60481729\:\:\:4817^{2}=23203489\:\:\:2034^{2}=04137156 $$ 就會得到:$4817\:\:2034\:\:1371,\:\:\cdots$ 如果有某個$n$位數字的平方值不滿$2n$位數的話，可以加入一個leadingzero使其變成$2n$位數，如上面例子中$2034^{2}=04137156$。

由於$n$為奇數時，會找不到唯一居中的$n$位數字，因此此方法只能產生偶數位的隨機數，假如遇到某些情況必須要由奇位數數字生成隨機數時，可以在初始值加入leadingzero，以產生偶數位數的數字，如$540^{2}=291600$，加入leadingzero後，$0540^{2}=00291600$，就可以得到唯一居中的數字。

作為最早期的隨機數生成方法，它的缺點顯而易見。

首先，它的週期普遍很短，且遇到一些特定數字時會產生非常短的重複週期，如:$0540→2916→5030→3009→0000→0000；24→57→24→57，\cdots$ 再者，透過統計方法檢定可以發現，這個方法生成的數列隨機性(分布均勻和線性獨立)並不高。

然而，使用此方法存取資料，可以大幅提升搜尋速度，這也是這個方法當時被開發出來的原因。

註:透過$Weyl$序列進行修正可以改進短周期的缺點，有興趣的讀者們可以詳閱以下文章。

https://arxiv.org/abs/1704.00358v5 線性同餘法(Linearcongruentialgenerator) 上世紀最廣泛使用的偽隨機數生成演算法，除了易於理解，且由於其原理是透過一個簡單的線性遞推函數產生序列，對於計算機的計算和存儲資料能力要求不高，因此在當時許多模擬過程都採用此方法。

如上所說線性同餘法是透過一線性函數對某除數作餘數除法得到序列，通式如下: $$ X_{i+1}=(aX_{i}+c)\:\:\:\:mod\:\:\:\:m $$ 其中$a(乘數,Multiplier)、c(增量或稱偏移量,Increment)、m(除數,Modulus)$皆為整數，如果需要介於$(0,1)$的隨機數，可將數列除以除數$m$。

看到這裡不知道大家會不會有一個疑惑，這樣生成方法能生成的數字只在$(1,m)$之間，如果沒有足夠大的$m$，應該很容易出現重複的數字吧? 沒錯，不是任一個線性函數和除數的組合都適合用於生成隨機數，其實要找到好的組合是相當困難的，好的組合可以產生性質良好的隨機數列，除了開頭提到的隨機性質(分布均勻和線性獨立)，同時我們希望初始值(seed)對序列的影響越小越好。

另外，我們希望生成的數列有長週期，長週期沒有明確定義，越長的週期就代表有更多的隨機數列可以選擇，此外，長週期也比短週期更能確保數列的分佈均勻，且透過數列間前後兩項$X_{i+1}、X_{i}$的自相關係數近似公式:$\rho=\frac{1}{a}-(\frac{6c}{am})(1-\frac{c}{m})\pm\frac{a}{m}$可得知最大值發生在$a=\sqrt{m}$，此時$\rho$約等於$\frac{2}{\sqrt{m}}$，當$m$越大數列間就越接近線性獨立，因此在能找到合適的乘數$a$搭配的情況下，我們會希望除數$m$越大越好，以得到長週期和良好的統計性質。

註:週期:經過循環再一次回到初始值所經過數字的個數。

然而我們要如何取得最大的週期呢? 根據Hull-Dobell定理可以得出以下結論: $$ 若一LCG具有全週期性，則\ \Leftrightarrowm、c為互質關係\ \Leftrightarrow若一質數p可整除m，則也可整除a-1\ \Leftrightarrow若4可整除m，則也可整除a-1 $$ 同時滿足這三個條件的$LCG$就有全週期性，如: $$ a=11,c=2,m=5,seed=1\ \ \Rightarrow3,0,2,4,1,3,0,2,4,1,\cdots $$ 從上述的定理不難看出，線性同餘方程對於$m$和$a$非常敏感，相較起來$c$只要不要和$m$有公因數，就能讓它們的公倍數有最大值以求得最大週期。

有些型態的$LCG$甚至會忽略$c$的設置，以減少計算量和內存消耗，如乘法型LCG(MultiplicaitiveLCG)，或稱$Lehmerrandomgenerator$。

註:$c>0$的$LCG$稱為混合型線性同餘產生器。

$MLCG$是$LCG$中$c=0$的特例，由於沒有常數的加減項，所以數列中必不能包含$0$，否則接下來的數列會退化成$0,0,0,\cdots$。

因此$MLCG$必不能達到全週期，僅有最大週期$m-1$，而取得最大週期的條件也和$LCG$有一些差別。

$$ 若一MLCG具有最大週期，則 \ \Rightarrowm為質數 \ \Rightarrowa是m的原根 $$ 假設$a、m$互質，且存在一$v$為滿足$a^{v}\equiv1\:(\:mod\;m)$的最小值，則$v$稱為$a$對$m$的階數$(order)$，而當$v=m-1$時，我們即稱$a$是$m$的原根$(primitiveroot)$。

如以下的例子，通過觀察能找出$v=2=m-1$，由此可知$2$是$3$的一個原根。

$$ a=2、m=3 $$ $$ 2^{1}\equiv2\:(\:mod\;3) \ 2^{2}\equiv1\:(\:mod\;3) $$ 多數的程式語言或數據分析軟體都有產生隨機數的函式，不同的語言、軟體因為用途不同，所選擇的生成器也會有所不同。

線性同餘法概念簡單、易於清楚表達、演算速度快，但對常數的選擇嚴格，且容易因為週期不夠長而無法避免數列間的相關性，因此較適用於快速簡便但統計性質較為粗糙的應用，如:C語言中常用的rand函式、Java的java.util.Random和其相關函式。

梅森旋轉法(MersenneTwister) 上面我們提過長週期的重要性和短週期的缺點，當我們需要產生大量的隨機數列時，不夠長的週期就無法提供良好的隨機性。

大部分NPC問題的計算量過於龐大而難以計算出最佳解，此時只能找出實務上較佳的解法來解決問題。

例如著名的最佳化問題TravellingSalesmanProblem，內容是要找出一個銷售人員從某個城市出發，隨後遍歷多個城市，最後回到原城市的最短路徑，其中每個城市和城市間的距離都不盡相同。

假設銷售員需要經過$4$個城市，只有$3!=6$種走法，但當需要經過$20$個城市時，就變成了$19!=121645100408832000$種，這類複雜度以指數或階乘函數增加的問題往往需要模擬。

作為生成隨機序列的來源，簡單的線性遞推式無法產生如此長的週期，模擬的品質肯定就不會高。

梅森旋轉產生器基於二進位制的線性旋轉反饋移位暫存產生器(TwistedGeneralizedFeedbackShiftRegister,TGFSR)產生隨機數列，$TGFSR$和廣義的反饋移位暫存器(GeneralizedFeedbackShiftRegister,GFSR)一樣都是由初始值透過一個遞推式產生一個長度為$n$的數列，將取得的一串數列(旋轉鍊)再進行線性反饋遞推，而$TGFSR$和$GFSR$不同的地方就在於取得序列後再次遞推的遞推式，$GFSR$的遞推規則是遍歷旋轉鍊將各個數以二進制表示並依照$x[l]\leftarrowx[(l+m)modn]\bigoplusx[l]$轉換，也就是把一個位元的狀態和某些位元進行異或運算後的結果回傳以產生新的隨機數列，其中$m、n$都是常數，$m$是控制週期的參數；$n$是遞迴的長度。

而$TGFSR$是在$GFSR$的基礎上將遞推式改寫成: $$ x[l]\leftarrowx[(l+m)\;mod\;\;n]\;\bigoplus\;shiftright(x[l])\;\bigoplus \begin{cases} 0&if\;\;LSB\;\;of\;\;x[l]=0\\ a&if\;\;LSB\;\;of\;\;x[l]=1 \end{cases} $$ 和$GFSR$的遞推式概念大概相同，不過在進行異或運算前會先對該二進位制數進行右移$1$個位元的運算，如$76$這個數字在二進制時是$1001100$，進行右移運算後就變成$0100110$也就是$38$，其實某數字在二進位制時全部右移一位，也就相當於在十進位時除以$2$，這非常直觀。

不過發明者並沒有這個意思，只是希望透過讓$2$個二進制數字間位元狀態差異更大，再做異或運算時能讓結果更加沒有規律。

註:異或運算是數位邏輯中的一種邏輯運算，記作XOR或$\bigoplus$，規則大致上是兩運算元不同時取值為真；相同時取值為假，在二進制中位元狀態相同即產生$0$；不同則產生$1$例:$100010\bigoplus010010$就會產生$110000$。

除此之外這個算法還多加了一個異或運算，在式子的最右邊可以看到，假如$x[l]$這個二進位制數的最低有效位(LeastSignificantBit)是$0$就和$0$這個數字進行異或運算；如果是$1$就和$a$這個常數轉換成二進制後進行異或運算，其中最低有效位其就是一個數字轉成二進制後最右邊的位元，也就是代表$2$的$0$次方的位元，因此$x[l]$的最低有效位等於$1$還是等於$0$這個問題也等價於$x[l]$這個數字是奇數還是偶數。

註:這種能接收函式做為引數或能回傳函式做為結果回傳的函式稱為高階函式(Higherorderfunction)。

整體來說$TGFSR$是$GFSR$一種變體，它們都具有非常長的週期性質，能產生高品質的隨機數列，有良好的統計性質，能通過大部分的統計檢定，如:k-stest、runtest、weightdistributiontest，其中weightdistributiontest在這裡是對於二進制序列獨立性的檢定，能用於生成器的$TGFSR$都能通過如用於梅森旋轉產生器的$TGFSR$，以MT19937-32為例，它的部分引數為$(w,n,m,r)=(32,624,397,31)$，其中$w$是二進制數字最大長度，也就是說這個算法不會產生比$2^{32}$還大的數字；$r$是低位遮罩，可以把它想成一個區分點，$r=31$就代表第$0$到第$31$位元都是低位。

除此之外，$TGFSR$在存取效率、初始化的便利程度也都有很大的提升。

現在大多數的程式語言和軟體中都有用梅森旋轉產生器產生隨機數的函式庫，在$R$、$Python$、$Matlab$、$PHP$、$Julia$中都是預設的隨機數產生器，目前有兩個版本$MT-19937-32$和$MT-19937-64$，其中$32$和$64$代表的是二進制數字最大長度這在上面有提過，而$19937$代表它們能產生週期為$2^{19937}-1$的數列，其中$2^{19937}-1$是一個梅森質數，這也是梅森旋轉產生器的名稱的由來。

參考資料:https://en.wikipedia.org/wiki/Pseudorandom_number_generatorhttps://en.wikipedia.org/wiki/Mersenne_Twisterhttps://en.wikipedia.org/wiki/Mersenne_primehttps://en.wikipedia.org/wiki/Linear_congruential_generatorhttp://commons.apache.org/proper/commons-rng/https://www.youtube.com/watch?v=1UCaZjdRC_chttps://en.wikipedia.org/wiki/Travelling_salesman_problemHull,T.E.andA.R.Dobell“RandomNumberGenerators.”MakotoMatsumoto“TwistedGFSRGeneratorsApril27,1992(RevisedNov.15).” oswald Previous Previous 機率積分轉換(UniversalityofUniform)