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

隨機(jī)數(shù)生成及其在統(tǒng)計(jì)模擬中的應(yīng)用黃湘云 關(guān)鍵詞：隨機(jī)數(shù); 統(tǒng)計(jì)檢驗(yàn); 模擬 審稿：郎大為、邊蓓蕾；編輯：吳佳萍揭秘統(tǒng)計(jì)軟件如 R，Octave，Matlab 等使用的隨機(jī)數(shù)發(fā)生器，然后做一些統(tǒng)計(jì)檢驗(yàn)，再將其應(yīng)用到獨(dú)立隨機(jī)變量和的模擬中，最后與符號計(jì)算得到的精確結(jié)果比較。除特別說明外，文中涉及到的隨機(jī)數(shù)都是指偽隨機(jī)數(shù)，發(fā)生器都是指隨機(jī)數(shù)發(fā)生器。背景隨機(jī)數(shù)的產(chǎn)生和檢驗(yàn)方法是蒙特卡羅方法的重要部分，另外兩個(gè)是概率分布抽樣方法和降低方差提高效率方法。在 20 世紀(jì) 40 年代中期，當(dāng)時(shí)為了原子彈的研制，烏拉姆（S.Ulam）、馮諾依曼（J.von Neumann） 和梅特羅波利斯（N. Metropolis） 在美國核武器研究實(shí)驗(yàn)室創(chuàng)立蒙特卡羅方法。當(dāng)時(shí)出于保密的需要，與隨機(jī)模擬相關(guān)的技術(shù)就代號 “蒙特卡羅”。早期取得的成果有產(chǎn)生隨機(jī)數(shù)的平方取中方法，取舍算法和逆變換法等。這兩個(gè)算法的內(nèi)容見統(tǒng)計(jì)之都王夜笙的文章。隨機(jī)數(shù)生成講隨機(jī)數(shù)發(fā)生器，不得不提及一個(gè)名為 Mersenne Twister（簡稱 MT）的發(fā)生器，它的周期長達(dá) 2199371，現(xiàn)在是R 、Octave 和 Matlab等軟件（較新版本）的默認(rèn)隨機(jī)數(shù)發(fā)生器。Matlab通過內(nèi)置的rng函數(shù)指定不同的發(fā)生器，其中包括1995年Matlab采用 George Marsaglia 在1991年提出的借位減（subtract with borrow，簡稱SWB）發(fā)生器。在Matlab中，設(shè)置如下命令可指定發(fā)生器及其狀態(tài)，其中1234是隨機(jī)數(shù)種子，指定發(fā)生器的狀態(tài)，目的是重復(fù)實(shí)驗(yàn)結(jié)果，v5uniform是發(fā)生器的名字。rng(1234, v5uniform) Octave 通過內(nèi)置的rand函數(shù)指定發(fā)生器的不同狀態(tài)，為獲取相同的兩組隨機(jī)數(shù)，state 參數(shù)得設(shè)置一樣，如1234（你也可以設(shè)置為別的值）。Octave已經(jīng)放棄了老版本內(nèi)置的發(fā)生器，找不到命令去指定早期的發(fā)生器，這個(gè)和 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ā)生器，類似地，通過如下方式獲得相同的兩組隨機(jī)數(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的默認(rèn)發(fā)生器也是 MT 發(fā)生器，除了設(shè)置隨機(jī)數(shù)種子的 seed 參數(shù)，還可以通過kind參數(shù)指定其他發(fā)生器，normal.kind參數(shù)指定產(chǎn)生正態(tài)分布隨機(jī)數(shù)的發(fā)生器，下面也使用類似的方式產(chǎn)生兩組完全一樣的隨機(jī)數(shù)。set.seed(seed, kind = NULL, normal.kind = NULL)set.seed(1234,kind = Mersenne-Twister, normal.kind = Inversion) # 默認(rèn)參數(shù)設(shè)置runif(5)1 0.1137034 0.6222994 0.6092747 0.6233794 0.8609154set.seed(1234,kind = Mersenne-Twister, normal.kind = Inversion) # 默認(rèn)參數(shù)設(shè)置runif(5)1 0.1137034 0.6222994 0.6092747 0.6233794 0.8609154SWB發(fā)生器中“借位相減”步驟是指序列的第i個(gè)隨機(jī)數(shù)zi要依據(jù)如下遞推關(guān)系產(chǎn)生，zi=zi+20zi+5b下標(biāo)i,i+20,i+5都是對32取模的結(jié)果，zi+20與zi+5做減法運(yùn)算，b是借位，其取值與前一步有關(guān)，當(dāng)zi是正值時(shí)，下一步將b置為0，如果計(jì)算的zi是負(fù)值，在保存zi之前，將其值加1，并在下一步，將b的值設(shè)為 253。下面關(guān)于隨機(jī)數(shù)生成的效率和后面的統(tǒng)計(jì)檢驗(yàn)，都以生成224個(gè)為基準(zhǔn)，是1600多萬個(gè)，取這么多，一方面為了比較編程語言實(shí)現(xiàn)的發(fā)生器產(chǎn)生隨機(jī)數(shù)的效率，另一方面是后面的游程檢驗(yàn)需要比較大的樣本量。Matlab內(nèi)置的發(fā)生器及大部分的函數(shù)，底層實(shí)現(xiàn)都是 C 或者 Fortran，MathWorks創(chuàng)始人Cleve B. Moler是數(shù)值分析領(lǐng)域著名的LINPACK和EISPACK包的作者之一，他當(dāng)年做了很多優(yōu)化和封裝，進(jìn)而推出Matlab，只要是調(diào)用內(nèi)置函數(shù)，效率不會(huì)比C差，自己寫的含有循環(huán)、判斷等語句的代碼，性能就因人而異了，對大多數(shù)人來說，性能要比C低。這里比較Matlab內(nèi)置SWB發(fā)生器（就當(dāng)作是C語言實(shí)現(xiàn)的好了和用Matlab重寫的SWB發(fā)生器的效率,代碼如下：% matlab codetic % 大約幾秒rng(1234, v5uniform) % 調(diào)用SWB發(fā)生器x = rand(1,224);toc% octave codeid = tic % 時(shí)間耗費(fèi)大約一小時(shí)randtx(state,0)x = randtx(1,224);toc (id)randtx不是Matlab和Octave內(nèi)置的函數(shù)，而是Cleve B.Moler基于Matlab實(shí)現(xiàn)的SWB發(fā)生器，也是100多行包含嵌套循環(huán)等語句的Matlab代碼打包的函數(shù)，上面的代碼運(yùn)行時(shí)間差異之大也就不難理解了，為了能在Octave上跑，我做了少量修改，Octave軟件版本為4.2.1，安裝Octave時(shí),Blas選擇OpenBlas，為了后續(xù)檢驗(yàn)，在獲得隨機(jī)數(shù)后，將其保存到磁盤文件random_number.matsave -mat random_number.mat x R，Octave，Matlab和Python內(nèi)置的發(fā)生器都是MT發(fā)生器，與之實(shí)現(xiàn)有關(guān)的數(shù)學(xué)庫，也是Blas，雖然有開源和進(jìn)一步優(yōu)化的商業(yè)版本之分，但是矩陣乘法，向量乘法之類運(yùn)算的效率也就半斤八兩，Julia語言官網(wǎng)給出了一個(gè)標(biāo)準(zhǔn)測試2。不同語言的性能表現(xiàn)（C語言在算法中的表現(xiàn)為基準(zhǔn)，時(shí)間記為1.0）這里再給出用C語言實(shí)現(xiàn)的MT發(fā)生器，產(chǎn)生同樣多的隨機(jī)數(shù)，只需要10秒左右，其中包含編譯和寫隨機(jī)數(shù)到文件的時(shí)間，實(shí)際產(chǎn)生隨機(jī)數(shù)的時(shí)間應(yīng)該遠(yuǎn)小于這個(gè)時(shí)間。（程序運(yùn)行環(huán)境環(huán)境Dev-C+ 5.11，用64位的GCC編譯）。統(tǒng)計(jì)檢驗(yàn)隨機(jī)數(shù)的檢驗(yàn)是有一套標(biāo)準(zhǔn)的，如George Marsaglia開發(fā)的DieHard檢驗(yàn)程序，檢驗(yàn)的內(nèi)容很豐富。這篇文章只能算初窺門徑，R內(nèi)產(chǎn)生真隨機(jī)數(shù)的包是 Dirk Eddelbuettel 開發(fā)的random包，它是連接產(chǎn)生真隨機(jī)數(shù)網(wǎng)站的接口。相關(guān)性檢驗(yàn)先來一個(gè)簡單的，就用R產(chǎn)生的兩個(gè)獨(dú)立同均勻分布樣本，調(diào)用cor.test做相關(guān)性檢驗(yàn)，看看這兩組數(shù)是不是足夠獨(dú)立同分布，通過眼球驗(yàn)證，隨著樣本量增大，相關(guān)性趨于0，相關(guān)性弱到可以視為獨(dú)立。如下圖所示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) 分布檢驗(yàn)檢驗(yàn)產(chǎn)生的隨機(jī)數(shù)是否服從指定的分布：原假設(shè)是樣本來自指定的分布，計(jì)算的P值比較大，就不能拒絕原假設(shè)。ks.test(runif(1000), punif) #分布檢驗(yàn)# One-sample Kolmogorov-Smirnov test# data: runif(1000)# D = 0.022302, p-value = 0.7025# alternative hypothesis: two-sided檢驗(yàn)兩樣本是否來自同一分布：原假設(shè)是兩樣本來自同一分布，計(jì)算的P值比較小，就表示兩樣本不是來自同一分布。ks.test(runif(1000), runif(1000) #同分布檢驗(yàn)# 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)置的隨機(jī)數(shù)發(fā)生器通過了檢驗(yàn)（嘿嘿，這是肯定的?。?。游程檢驗(yàn)游程檢驗(yàn)對隨機(jī)數(shù)序列的隨機(jī)性檢驗(yàn)，不對序列做任何分布假設(shè)，不同于上面的相關(guān)性檢驗(yàn)和省略沒講的特征統(tǒng)計(jì)量（如均值和方差）的檢驗(yàn)。先對隨機(jī)性略作解釋，簡單起見，我們考慮0-1序列，拋擲均勻的硬幣1000次，正面向上記為1，反面向上記為0，這是一個(gè)離散的均勻分布，每一次拋擲硬幣都無法準(zhǔn)確地判斷出現(xiàn)的是正面還是反面，若記錄的序列中0和1相對集中的出現(xiàn)，不是隨機(jī)，0和1交替出現(xiàn)，呈現(xiàn)周期性也不是隨機(jī)，除了這兩種情況基本就是隨機(jī)了。游程檢驗(yàn)的原假設(shè)是序列滿足隨機(jī)性，序列一旦生成，就是有序的，因?yàn)橛纬虣z驗(yàn)需要固定位置，再數(shù)上升（下降）的游程數(shù)，當(dāng)計(jì)算的P值比較大時(shí)，不能拒絕原假設(shè)，即不能否認(rèn)這個(gè)序列是隨機(jī)的。對上述0-1序列進(jìn)行模擬，然后檢驗(yàn)，如下所示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)在用游程檢驗(yàn)比較SWB發(fā)生器（Octave/Matlab版）、MT發(fā)生器（R語言版）和MT發(fā)生器（C語言版）。對于一般的實(shí)數(shù)序列，得先指定一個(gè)閾值，記為delta，然后比較序列中的值和delta的大小關(guān)系，這里類似數(shù)上升或下降的游程長度（runs length），基于這樣一個(gè)事實(shí)：如果序列表現(xiàn)隨機(jī)，則序列中兩個(gè)小于delta的值，間隔越遠(yuǎn)出現(xiàn)的次數(shù)越少?？梢赃@樣理解，還是以上面拋硬幣的例子來說，出現(xiàn)了很多次正面，那么下一次拋擲傾向于出現(xiàn)反面，這是一條人人可接受的經(jīng)驗(yàn)。為了把這條經(jīng)驗(yàn)可視化出來，對序列做如下操作：先給定閾值delta為0.01（也可以取別的值），取出序列中的值小于delta的位置，位置序列前面添加0，再差分，然后每個(gè)值減1，得到新序列，新序列中的值為0，表明原序列連續(xù)兩個(gè)值小于delta，新序列中的值為1，表明間隔1個(gè)數(shù)小于delta，新序列中的值為2，表明間隔2個(gè)數(shù)小于delta，依次類推. 統(tǒng)計(jì)所有的情況，用直方圖顯示，這就獲得游程長度與間隔的關(guān)系圖（間隔數(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默認(rèn)采用Mersenne Twister發(fā)生器r_data - runif(224, 0, 1); # R內(nèi)生成均勻分布隨機(jī)數(shù)matlabv5_data - readMat(random_number.mat) # 讀取Octave生成的均勻分布隨機(jī)數(shù)temp - read.table(file = random_number.txt) # 讀取C語言生成的均勻分布隨機(jī)數(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ā)生器通過了檢驗(yàn)，SWT發(fā)生器沒有通過，在間隔數(shù)為27的位置，有一條溝，按規(guī)律游程長度不應(yīng)該減少這么多，這是因?yàn)镾WB發(fā)生器 “借位減”步驟，取模32的運(yùn)算和5一起消耗了間隔為27的數(shù)量（讀者可以按借位減的遞推關(guān)系思考是如何消耗的），導(dǎo)致不符合隨機(jī)性的要求，該算法細(xì)節(jié)參見Cleve B. Moler的書Numerical Computing with MATLAB第9章第267頁 。應(yīng)用兩個(gè)均勻分布的統(tǒng)計(jì)模擬隨機(jī)變量X1,X2iid某分布（比如二項(xiàng)分布，泊松分布，正態(tài)分布，指數(shù)分布，卡方分布，伽馬分布），則X1+X2也服從該分布。常見的均勻分布是否具有這樣的可加性？具體地說，就是X1,X2iidU(0,1),X1+X2是否服從U(0,2)？如果有一臺(tái)電腦在旁邊，我首先想到的就是敲三五行代碼，畫個(gè)散點(diǎn)圖、直方圖，看圖說話。set.seed(1234) x - runif(10000, min = 0, max = 1) y - runif(10000, min = 0, max = 1) z - x + yplot(z) # 散點(diǎn)圖hist(z) # 直方圖為美觀起見，從viridis包調(diào)用viridis調(diào)色板，顏色越深的地方，相應(yīng)的數(shù)值越大，不管是此處geom_hex繪制的六角形熱圖，還是geom_histogram繪制的直方圖，都遵循這個(gè)規(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處，散點(diǎn)比較集中，看起來有點(diǎn)像正態(tài)分布。如果往中心極限定理上靠，將作如下標(biāo)準(zhǔn)化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)在看來眼球檢驗(yàn)不好使了，那就上P值唄！ks.test(sqrt(6)*(z-1), pnorm) # 分布檢驗(yàn)# One-sample Kolmogorov-Smirnov test# data: sqrt(6) * (z - 1)# D = 0.025778, p-value = 3.381e-06# alternative hypothesis: two-sided也不是正態(tài)分布，既然如此，那就在兩個(gè)隨機(jī)變量的情況下，把精確分布推導(dǎo)出來。精確分布的推導(dǎo)及計(jì)算課本如概率論與數(shù)理統(tǒng)計(jì)教程 采用卷積的方法求分布函數(shù)，這種方法實(shí)行起來比較繁瑣，也不利于后續(xù)編程，下面考慮用特征函數(shù)的方法求。我們知道標(biāo)準(zhǔn)均勻分布的特征函數(shù)(t)=eit1it考慮X1和X2相互獨(dú)立，它們的和用S2表示，則隨機(jī)變量S2的特征函數(shù)為2(t)=(t)(t)=(eit1it)2=2(1cos(t)eitt2只要滿足條件+|2(t)|dtS2的密度函數(shù)就可以表示為p2(x)=12+eitx2(t)dt經(jīng)計(jì)算+|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個(gè)獨(dú)立隨機(jī)變量的和n(t)=(eit1it)n=(sin(t/2)eit2t/2)n那么，同理pn(x)=2+0cos(2(n/2x)t)(sin(t)t)ndt要說數(shù)值計(jì)算一個(gè)p(x)近似值，是一點(diǎn)問題沒有！且看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那如果要把上面的積分積出來，獲得一個(gè)精確的表達(dá)式，在n=2的時(shí)候還可以手動(dòng)計(jì)算，主要使用分部積分，余弦積化和差公式和一個(gè)狄利克雷積分公式+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ù)值計(jì)算的結(jié)果上很一致，當(dāng)n=100,1000時(shí)，含參量積分的表示形式就很方便啦！任意給定一個(gè)n，符號計(jì)算上面的含參量積分，這個(gè)時(shí)候還是用軟件計(jì)算比較合適，R 的符號計(jì)算僅限于求導(dǎo)，積分運(yùn)算需要借助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稍為好點(diǎn)，但是還是有一大塊看不懂，那個(gè)絕對值里是什么 3？還是不要糾結(jié)了，路遠(yuǎn)坑多，慢走不送?。≡捳f要是計(jì)算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