偽隨機數

01-25

隨機模擬的基礎就是能夠生成隨機數。不過這裡說的隨機數，不是任意一個隨機變數，而是特指［0，1］區間上的服從均勻分布 $U[0,1]$ 的隨機變數。我們這裡解釋隨機數如何用計算機產生的並且如何運行隨機數計算定積分的。

事實上，不可能用計算機生成隨機數的。凡用計算機生成的一列數，肯定是用程序實現的，既然用程序實現，肯定不是隨機的。因此用計算機生成的所謂隨機數只能稱為偽隨機數。只不過用計算機生成的偽隨機數真正的隨機數相比較沒有顯著性差異，可以作為隨機數來用。

下面我們給出常見的偽隨機數生成演算法－線性同餘法。雖然現在偽隨機數的生成演算法各式各樣並且越來越複雜，但思想和該演算法基本上是一樣的。因此通過對該演算法的理解，就能理解偽隨機數的產生原理了。

線性同餘法： $x_n=ax_{n-1}+c mod m$

例子1. $x_0=5,x_n=3x_{n-1}+7 mod 200$ ，用R語言給出前100個迭代結果。

random1<-function(n){ x<-0 x[1]<-5 for(i in 2:n){ x[i]=3*x[i-1]+7 x[i]=x[i]%%200 } x }random1(100)

運行結果如下：

[1] 5 22 73 26 85 62 193 186 165 102 113 146 45 142 33 106 125 182

[19] 153 66 5 22 73 26 85 62 193 186 165 102 113 146 45 142 33 106

[37] 125 182 153 66 5 22 73 26 85 62 193 186 165 102 113 146 45 142

[55] 33 106 125 182 153 66 5 22 73 26 85 62 193 186 165 102 113 146

[73] 45 142 33 106 125 182 153 66 5 22 73 26 85 62 193 186 165 102

[91] 113 146 45 142 33 106 125 182 153 66

從程序來看，生成的100個正整數，因此肯定不是100個[0,1]上的偽隨機數。那麼怎麼辦？把每個數除以200。這樣就得到了100個偽隨機數。

random1(100)/200

運行結果如下：

[1] 0.025 0.110 0.365 0.130 0.425 0.310 0.965 0.930 0.825 0.510 0.565 0.730

[13] 0.225 0.710 0.165 0.530 0.625 0.910 0.765 0.330 0.025 0.110 0.365 0.130

[25] 0.425 0.310 0.965 0.930 0.825 0.510 0.565 0.730 0.225 0.710 0.165 0.530

[37] 0.625 0.910 0.765 0.330 0.025 0.110 0.365 0.130 0.425 0.310 0.965 0.930

[49] 0.825 0.510 0.565 0.730 0.225 0.710 0.165 0.530 0.625 0.910 0.765 0.330

[61] 0.025 0.110 0.365 0.130 0.425 0.310 0.965 0.930 0.825 0.510 0.565 0.730

[73] 0.225 0.710 0.165 0.530 0.625 0.910 0.765 0.330 0.025 0.110 0.365 0.130

[85] 0.425 0.310 0.965 0.930 0.825 0.510 0.565 0.730 0.225 0.710 0.165 0.530

[97] 0.625 0.910 0.765 0.330

從上面迭代的結果來看，運行結果具有周期性。5出現了5次（加粗），周期為20。這樣出現的隨機數肯定也具有周期性。顯然這與隨機數的特徵是不符的。一列隨機數相同的概率是為0的。因此偽隨機數不應該具有一定的周期。事實上，周期性是無法避免的。如果生成的偽隨機數的個數大於m，肯定會有周期發生。如何避免周期性呢？第一m盡量大，第二周期為儘可能做到等於m。是否存在無論取什麼初值，其生成的隨機數的周期為m呢？下面定理解決了這個問題.

定理線性同餘法周期為m當且僅當滿足下面三個條件：（1）m和c互素；（2）a-1可以被m的每個素因子整除；（3）如果m是4的整數倍，那麼a-1也是4的整數倍。

例子2 $m=2^{32},a=1103515245,c=12345$ 。下面驗證滿足定理的3個條件。由於m只有一個素數因子2，且c為奇數，因此m和c互素，這樣條件（1）滿足。顯然條件（2）也滿足。顯然m是4的倍數，而a-1=1103515244=275878811*4，因此條件（3）也滿足。

現在是 R 、Octave 和 Matlab 等軟體（較新版本）的默認隨機數發生器是一個名為 Mersenne Twister（簡稱 MT）的發生器，它的周期長達 $2^{19937}?1$ 。