隨機數(shù)生成及其在統(tǒng)計模擬中的應用.doc

隨機數(shù)生成及其在統(tǒng)計模擬中的應用黃湘云 關鍵詞：隨機數(shù); 統(tǒng)計檢驗; 模擬 審稿：郎大為、邊蓓蕾；編輯：吳佳萍揭秘統(tǒng)計軟件如 R，Octave，Matlab 等使用的隨機數(shù)發(fā)生器，然后做一些統(tǒng)計檢驗，再將其應用到獨立隨機變量和的模擬中，最后與符號計算得到的精確結(jié)果比較。除特別說明外，文中涉及到的隨機數(shù)都是指偽隨機數(shù)，發(fā)生器都是指隨機數(shù)發(fā)生器。背景隨機數(shù)的產(chǎn)生和檢驗方法是蒙特卡羅方法的重要部分，另外兩個是概率分布抽樣方法和降低方差提高效率方法。在 20 世紀 40 年代中期，當時為了原子彈的研制，烏拉姆（S.Ulam）、馮諾依曼（J.von Neumann） 和梅特羅波利斯（N. Metropolis） 在美國核武器研究實驗室創(chuàng)立蒙特卡羅方法。當時出于保密的需要，與隨機模擬相關的技術就代號 “蒙特卡羅”。早期取得的成果有產(chǎn)生隨機數(shù)的平方取中方法，取舍算法和逆變換法等。這兩個算法的內(nèi)容見統(tǒng)計之都王夜笙的文章。隨機數(shù)生成講隨機數(shù)發(fā)生器，不得不提及一個名為 Mersenne Twister（簡稱 MT）的發(fā)生器，它的周期長達 2199371，現(xiàn)在是R 、Octave 和 Matlab等軟件（較新版本）的默認隨機數(shù)發(fā)生器。Matlab通過內(nèi)置的rng函數(shù)指定不同的發(fā)生器，其中包括1995年Matlab采用 George Marsaglia 在1991年提出的借位減（subtract with borrow，簡稱SWB）發(fā)生器。在Matlab中，設置如下命令可指定發(fā)生器及其狀態(tài)，其中1234是隨機數(shù)種子，指定發(fā)生器的狀態(tài)，目的是重復實驗結(jié)果，v5uniform是發(fā)生器的名字。rng(1234, v5uniform) Octave 通過內(nèi)置的rand函數(shù)指定發(fā)生器的不同狀態(tài)，為獲取相同的兩組隨機數(shù)，state 參數(shù)得設置一樣，如1234（你也可以設置為別的值）。Octave已經(jīng)放棄了老版本內(nèi)置的發(fā)生器，找不到命令去指定早期的發(fā)生器，這個和 Matlab 不一樣。rand (state,1234)rand(1,5) 0.9664535 0.4407326 0.0074915 0.9109760 0.9392690rand (state,1234)rand(1,5) 0.9664535 0.4407326 0.0074915 0.9109760 0.9392690Python的numpy模塊也采用MT發(fā)生器，類似地，通過如下方式獲得相同的兩組隨機數(shù)import numpy as npa = np.random.RandomState(1234)a.rand(5)array( 0.19151945, 0.62210877, 0.43772774, 0.78535858, 0.77997581)a = np.random.RandomState(1234)a.rand(5)array( 0.19151945, 0.62210877, 0.43772774, 0.78535858, 0.77997581)R的默認發(fā)生器也是 MT 發(fā)生器，除了設置隨機數(shù)種子的 seed 參數(shù)，還可以通過kind參數(shù)指定其他發(fā)生器，normal.kind參數(shù)指定產(chǎn)生正態(tài)分布隨機數(shù)的發(fā)生器，下面也使用類似的方式產(chǎn)生兩組完全一樣的隨機數(shù)。set.seed(seed, kind = NULL, normal.kind = NULL)set.seed(1234,kind = Mersenne-Twister, normal.kind = Inversion) # 默認參數(shù)設置runif(5)1 0.1137034 0.6222994 0.6092747 0.6233794 0.8609154set.seed(1234,kind = Mersenne-Twister, normal.kind = Inversion) # 默認參數(shù)設置runif(5)1 0.1137034 0.6222994 0.6092747 0.6233794 0.8609154SWB發(fā)生器中“借位相減”步驟是指序列的第i個隨機數(shù)zi要依據(jù)如下遞推關系產(chǎn)生，zi=zi+20zi+5b下標i,i+20,i+5都是對32取模的結(jié)果，zi+20與zi+5做減法運算，b是借位，其取值與前一步有關，當zi是正值時，下一步將b置為0，如果計算的zi是負值，在保存zi之前，將其值加1，并在下一步，將b的值設為 253。下面關于隨機數(shù)生成的效率和后面的統(tǒng)計檢驗，都以生成224個為基準，是1600多萬個，取這么多，一方面為了比較編程語言實現(xiàn)的發(fā)生器產(chǎn)生隨機數(shù)的效率，另一方面是后面的游程檢驗需要比較大的樣本量。Matlab內(nèi)置的發(fā)生器及大部分的函數(shù)，底層實現(xiàn)都是 C 或者 Fortran，MathWorks創(chuàng)始人Cleve B. Moler是數(shù)值分析領域著名的LINPACK和EISPACK包的作者之一，他當年做了很多優(yōu)化和封裝，進而推出Matlab，只要是調(diào)用內(nèi)置函數(shù)，效率不會比C差，自己寫的含有循環(huán)、判斷等語句的代碼，性能就因人而異了，對大多數(shù)人來說，性能要比C低。這里比較Matlab內(nèi)置SWB發(fā)生器（就當作是C語言實現(xiàn)的好了和用Matlab重寫的SWB發(fā)生器的效率,代碼如下：% matlab codetic % 大約幾秒rng(1234, v5uniform) % 調(diào)用SWB發(fā)生器x = rand(1,224);toc% octave codeid = tic % 時間耗費大約一小時randtx(state,0)x = randtx(1,224);toc (id)randtx不是Matlab和Octave內(nèi)置的函數(shù)，而是Cleve B.Moler基于Matlab實現(xiàn)的SWB發(fā)生器，也是100多行包含嵌套循環(huán)等語句的Matlab代碼打包的函數(shù)，上面的代碼運行時間差異之大也就不難理解了，為了能在Octave上跑，我做了少量修改，Octave軟件版本為4.2.1，安裝Octave時,Blas選擇OpenBlas，為了后續(xù)檢驗，在獲得隨機數(shù)后，將其保存到磁盤文件random_number.matsave -mat random_number.mat x R，Octave，Matlab和Python內(nèi)置的發(fā)生器都是MT發(fā)生器，與之實現(xiàn)有關的數(shù)學庫，也是Blas，雖然有開源和進一步優(yōu)化的商業(yè)版本之分，但是矩陣乘法，向量乘法之類運算的效率也就半斤八兩，Julia語言官網(wǎng)給出了一個標準測試2。不同語言的性能表現(xiàn)（C語言在算法中的表現(xiàn)為基準，時間記為1.0）這里再給出用C語言實現(xiàn)的MT發(fā)生器，產(chǎn)生同樣多的隨機數(shù)，只需要10秒左右，其中包含編譯和寫隨機數(shù)到文件的時間，實際產(chǎn)生隨機數(shù)的時間應該遠小于這個時間。（程序運行環(huán)境環(huán)境Dev-C+ 5.11，用64位的GCC編譯）。統(tǒng)計檢驗隨機數(shù)的檢驗是有一套標準的，如George Marsaglia開發(fā)的DieHard檢驗程序，檢驗的內(nèi)容很豐富。這篇文章只能算初窺門徑，R內(nèi)產(chǎn)生真隨機數(shù)的包是 Dirk Eddelbuettel 開發(fā)的random包，它是連接產(chǎn)生真隨機數(shù)網(wǎng)站的接口。相關性檢驗先來一個簡單的，就用R產(chǎn)生的兩個獨立同均勻分布樣本，調(diào)用cor.test做相關性檢驗，看看這兩組數(shù)是不是足夠獨立同分布，通過眼球驗證，隨著樣本量增大，相關性趨于0，相關性弱到可以視為獨立。如下圖所示library(ggplot2)library(viridisLite)library(viridis)set.seed(1234)corr - rep(0, 1000) for(i in seq(from = 1000, to = 1000000, by = 1000) corri/1000 - cor.test(runif(i, min = 0, max = 1), runif(i, min = 0, max = 1)$estimate ggplot(data.frame(x = seq(1000), y = corr), aes(x = x, y = y) + geom_hex(show.legend = FALSE) + scale_fill_viridis(direction = -1) + xlab(Sample size *103) + ylab(Correlation) 分布檢驗檢驗產(chǎn)生的隨機數(shù)是否服從指定的分布：原假設是樣本來自指定的分布，計算的P值比較大，就不能拒絕原假設。ks.test(runif(1000), punif) #分布檢驗# One-sample Kolmogorov-Smirnov test# data: runif(1000)# D = 0.022302, p-value = 0.7025# alternative hypothesis: two-sided檢驗兩樣本是否來自同一分布：原假設是兩樣本來自同一分布，計算的P值比較小，就表示兩樣本不是來自同一分布。ks.test(runif(1000), runif(1000) #同分布檢驗# Two-sample Kolmogorov-Smirnov test# data: runif(1000) and runif(1000)# D = 0.04, p-value = 0.4005# alternative hypothesis: two-sided從結(jié)果來看,R內(nèi)置的隨機數(shù)發(fā)生器通過了檢驗（嘿嘿，這是肯定的?。?。游程檢驗游程檢驗對隨機數(shù)序列的隨機性檢驗，不對序列做任何分布假設，不同于上面的相關性檢驗和省略沒講的特征統(tǒng)計量（如均值和方差）的檢驗。先對隨機性略作解釋，簡單起見，我們考慮0-1序列，拋擲均勻的硬幣1000次，正面向上記為1，反面向上記為0，這是一個離散的均勻分布，每一次拋擲硬幣都無法準確地判斷出現(xiàn)的是正面還是反面，若記錄的序列中0和1相對集中的出現(xiàn)，不是隨機，0和1交替出現(xiàn)，呈現(xiàn)周期性也不是隨機，除了這兩種情況基本就是隨機了。游程檢驗的原假設是序列滿足隨機性，序列一旦生成，就是有序的，因為游程檢驗需要固定位置，再數(shù)上升（下降）的游程數(shù)，當計算的P值比較大時，不能拒絕原假設，即不能否認這個序列是隨機的。對上述0-1序列進行模擬，然后檢驗，如下所示library(tseries)x - sample(c(0, 1), 1000, replace = TRUE, prob = c(1/2, 1/2)runs.test(factor(x)# Runs Test# data: factor(x)# Standard Normal = 0.45116, p-value = 0.6519# alternative hypothesis: two.sided現(xiàn)在用游程檢驗比較SWB發(fā)生器（Octave/Matlab版）、MT發(fā)生器（R語言版）和MT發(fā)生器（C語言版）。對于一般的實數(shù)序列，得先指定一個閾值，記為delta，然后比較序列中的值和delta的大小關系，這里類似數(shù)上升或下降的游程長度（runs length），基于這樣一個事實：如果序列表現(xiàn)隨機，則序列中兩個小于delta的值，間隔越遠出現(xiàn)的次數(shù)越少。可以這樣理解，還是以上面拋硬幣的例子來說，出現(xiàn)了很多次正面，那么下一次拋擲傾向于出現(xiàn)反面，這是一條人人可接受的經(jīng)驗。為了把這條經(jīng)驗可視化出來，對序列做如下操作：先給定閾值delta為0.01（也可以取別的值），取出序列中的值小于delta的位置，位置序列前面添加0，再差分，然后每個值減1，得到新序列，新序列中的值為0，表明原序列連續(xù)兩個值小于delta，新序列中的值為1，表明間隔1個數(shù)小于delta，新序列中的值為2，表明間隔2個數(shù)小于delta，依次類推. 統(tǒng)計所有的情況，用直方圖顯示，這就獲得游程長度與間隔的關系圖（間隔數(shù)取到100足可示意）。library(gridExtra)library(R.matlab)# 游程頻數(shù)直方圖run_test_fun - function(x, string, delta) n - length(x) len - diff(c(0, which(x delta), n + 1) - 1 ggplot(data.frame(x = lenlen 101), aes(x, fill = .count.) + scale_fill_viridis(direction = -1) + geom_histogram(binwidth = 1, show.legend = FALSE) + xlab(string) + ylab() set.seed(1234) # R默認采用Mersenne Twister發(fā)生器r_data - runif(224, 0, 1); # R內(nèi)生成均勻分布隨機數(shù)matlabv5_data - readMat(random_number.mat) # 讀取Octave生成的均勻分布隨機數(shù)temp - read.table(file = random_number.txt) # 讀取C語言生成的均勻分布隨機數(shù) c_data - c(as.matrix(t(temp)p1 - run_test_fun(x = r_data, string = R, delta = 0.01)p2 - run_test_fun(x = matlabv5_data$x, string = Matlab v5, delta = 0.01)p3 - run_test_fun(x = c_data, string = C, delta = 0.01)grid.arrange(p1, p2, p3, ncol=3)從圖中可以看出MT發(fā)生器通過了檢驗，SWT發(fā)生器沒有通過，在間隔數(shù)為27的位置，有一條溝，按規(guī)律游程長度不應該減少這么多，這是因為SWB發(fā)生器 “借位減”步驟，取模32的運算和5一起消耗了間隔為27的數(shù)量（讀者可以按借位減的遞推關系思考是如何消耗的），導致不符合隨機性的要求，該算法細節(jié)參見Cleve B. Moler的書Numerical Computing with MATLAB第9章第267頁 。應用兩個均勻分布的統(tǒng)計模擬隨機變量X1,X2iid某分布（比如二項分布，泊松分布，正態(tài)分布，指數(shù)分布，卡方分布，伽馬分布），則X1+X2也服從該分布。常見的均勻分布是否具有這樣的可加性？具體地說，就是X1,X2iidU(0,1),X1+X2是否服從U(0,2)？如果有一臺電腦在旁邊，我首先想到的就是敲三五行代碼，畫個散點圖、直方圖，看圖說話。set.seed(1234) x - runif(10000, min = 0, max = 1) y - runif(10000, min = 0, max = 1) z - x + yplot(z) # 散點圖hist(z) # 直方圖為美觀起見，從viridis包調(diào)用viridis調(diào)色板，顏色越深的地方，相應的數(shù)值越大，不管是此處geom_hex繪制的六角形熱圖，還是geom_histogram繪制的直方圖，都遵循這個規(guī)律。ggplot(data.frame(x = seq(10000), y = z), aes(x = x, y = y) + geom_hex(show.legend = FALSE) + scale_fill_viridis(direction = -1) + xlab() + ylab()顯然這不是均勻分布，在 z=1處，散點比較集中，看起來有點像正態(tài)分布。如果往中心極限定理上靠，將作如下標準化Y2=X1+X22121122=6(X1+X21) 則Y2的期望為0，方差為1。p4 - ggplot(data.frame(x = z), aes(x, fill = .count.) + scale_fill_viridis(direction = -1) + geom_histogram(bins=20, show.legend = FALSE) + xlab() + ylab() p5 - ggplot(data.frame(x = sqrt(6)*(z-1), aes(x, fill = .count.) + scale_fill_viridis(direction = -1) + geom_histogram(bins = 20, show.legend = FALSE) + xlab() + ylab() grid.arrange(p4, p5, ncol=2)只是變換后的圖像和之前基本一致，那么現(xiàn)在看來眼球檢驗不好使了，那就上P值唄！ks.test(sqrt(6)*(z-1), pnorm) # 分布檢驗# One-sample Kolmogorov-Smirnov test# data: sqrt(6) * (z - 1)# D = 0.025778, p-value = 3.381e-06# alternative hypothesis: two-sided也不是正態(tài)分布，既然如此，那就在兩個隨機變量的情況下，把精確分布推導出來。精確分布的推導及計算課本如概率論與數(shù)理統(tǒng)計教程 采用卷積的方法求分布函數(shù)，這種方法實行起來比較繁瑣，也不利于后續(xù)編程，下面考慮用特征函數(shù)的方法求。我們知道標準均勻分布的特征函數(shù)(t)=eit1it考慮X1和X2相互獨立，它們的和用S2表示，則隨機變量S2的特征函數(shù)為2(t)=(t)(t)=(eit1it)2=2(1cos(t)eitt2只要滿足條件+|2(t)|dtS2的密度函數(shù)就可以表示為p2(x)=12+eitx2(t)dt經(jīng)計算+|2(t)|dt=4+01cos(t)t2dt=4+0(sin(x)x)2dx=2那么p2(x)=12+eitx2(t)dt=2+0(1cos(t)cos(t(1x)t2dt=2+0cos(2(1x)t)(sin(t)t)2dt一般地，n個獨立隨機變量的和n(t)=(eit1it)n=(sin(t/2)eit2t/2)n那么，同理pn(x)=2+0cos(2(n/2x)t)(sin(t)t)ndt要說數(shù)值計算一個p(x)近似值，是一點問題沒有！且看integrate(function(t,x,n) 2/pi*cos(n-2*x)*t)*(sin(t)/t)n ,x = 1,n = 2,lower = 0,upper = Inf,subdivisions = 1000) # 0.9999846 with absolute error6.6e-05那如果要把上面的積分積出來，獲得一個精確的表達式，在n=2的時候還可以手動計算，主要使用分部積分，余弦積化和差公式和一個狄利克雷積分公式+0sin(ax)xdx=2sgn(a)，過程略，最后算得p2(x)=12(2x)sgn(2x)xsgn(x)(1x)sgn(1x)=12(|x|+|x2|)|x1|,0x2p2(x)的密度函數(shù)圖象如下：fun_p2_1 - function(x) 1 / 2 * (abs(x - 2) - 2 * abs(x - 1) + abs(x) fun_p2_2 - function(x) x - as.matrix(x) tempfun - function(x) integrate(function(t, x, n) 2 / pi * cos(n - 2 * x) * t) * (sin(t) / t) n, x=x, n= 2,lower = 0, upper = Inf, subdivisions = 1000)$value return( sapply(x,tempfun) )ggplot(data.frame(x = c(0, 2), aes(x = x) + stat_function(fun = fun_p2_2, geom = point, colour = #2A768EFF) + stat_function(fun = fun_p2_1, geom = line, colour = #78D152FF) 從圖中可以看出，兩種形式的密度函數(shù)在數(shù)值計算的結(jié)果上很一致，當n=100,1000時，含參量積分的表示形式就很方便啦！任意給定一個n，符號計算上面的含參量積分，這個時候還是用軟件計算比較合適，R 的符號計算僅限于求導，積分運算需要借助Ryacas，rSymPy，可惜的是，這些包更新緩慢，即使+0sin(at)tdt也算不出來，果斷直接使用Python的sympy模塊from sympy import * a=symbols(a, real = True)t=symbols(t, real = True, positive = True)print(integrate(sin(a*t)/t, (t, 0, oo)# Piecewise(pi/2, Eq(Abs(periodic_argument(polar_lift(a)*2, oo), 0), (Integral(sin(a*t)/t, (t, 0, oo), True)初次見到這樣的結(jié)果，是不是一臉mb，翻譯一下，就是2forperiodicargument(polar_lift2(a),)=001tsin(at)dtotherwise稍為好點，但是還是有一大塊看不懂，那個絕對值里是什么 3？還是不要糾結(jié)了，路遠坑多，慢走不送??！話說要是計算p2(x)密度函數(shù)里的積分，from sympy import * x=symbols(x, real=True)t=symbols(t, real=True,positive=True)print(integrate(2/pi*cos(2*t*(1-x)*(sin(t)/t)*2,(t,0,oo)# Piecewise(Piecewise(2*x, (2*x - 2)*2/4 1), (0, 4/(2*x - 2)*2 1), (meijerg(1/2,), (1, 1, 3/2), (1/2, 1, 0), (1/2,), polar_lift(-2*x + 2)*2/4), True)/2, Eq(Abs(periodic_argument(polar_lift(-2*x + 2)*2, oo), 0), (Integral(2*sin(t)*2*cos(2*t*(-x + 1)/(pi*t*2), (t, 0, oo), True)那就更長了。2xfor14(2x2)210for4(2x2)21G3,14