随机数生成及其在统计模拟中的应用

1、-作者xxxx-日期xxxx随机数生成及其在统计模拟中的应用【精品文档】随机数生成及其在统计模拟中的应用黄湘云 关键词：随机数; 统计检验; 模拟 审稿：郎大为、边蓓蕾；编辑：吴佳萍揭秘统计软件如 R，Octave，Matlab 等使用的随机数发生器，然后做一些统计检验，再将其应用到独立随机变量和的模拟中，最后与符号计算得到的精确结果比较。除特别说明外，文中涉及到的随机数都是指伪随机数，发生器都是指随机数发生器。背景随机数的产生和检验方法是蒙特卡罗方法的重要部分，另外两个是概率分布抽样方法和降低方差提高效率方法。在 20 世纪 40 年代中期，当时为了原子弹的研制，乌拉姆（）、冯诺依曼（ Ne

2、umann） 和梅特罗波利斯（N. Metropolis） 在美国核武器研究实验室创立蒙特卡罗方法。当时出于保密的需要，与随机模拟相关的技术就代号 “蒙特卡罗”。早期取得的成果有产生随机数的平方取中方法，取舍算法和逆变换法等。这两个算法的内容见统计之都王夜笙的文章。随机数生成讲随机数发生器，不得不提及一个名为 Mersenne Twister（简称 MT）的发生器，它的周期长达 2199371，现在是R 、Octave 和 Matlab等软件（较新版本）的默认随机数发生器。Matlab通过内置的rng函数指定不同的发生器，其中包括1995年Matlab采用 George Marsaglia 在

3、1991年提出的借位减（subtract with borrow，简称SWB）发生器。在Matlab中，设置如下命令可指定发生器及其状态，其中1234是随机数种子，指定发生器的状态，目的是重复实验结果，v5uniform是发生器的名字。rng(1234, 'v5uniform') Octave 通过内置的rand函数指定发生器的不同状态，为获取相同的两组随机数，state 参数得设置一样，如1234（你也可以设置为别的值）。Octave已经放弃了老版本内置的发生器，找不到命令去指定早期的发生器，这个和 Matlab 不一样。rand ('state',1234)

4、rand(1,5)rand ('state',1234)rand(1,5)Python的numpy模块也采用MT发生器，类似地，通过如下方式获得相同的两组随机数import numpy as npa = (1234)(5)array( 0.19151945, 0.62210877, 0.43772774, 0.78535858, 0.77997581)a = (1234)(5)array( 0.19151945, 0.62210877, 0.43772774, 0.78535858, 0.77997581)R的默认发生器也是 MT 发生器，除了设置随机数种子的 seed 参数，

5、还可以通过kind参数指定其他发生器，参数指定产生正态分布随机数的发生器，下面也使用类似的方式产生两组完全一样的随机数。(seed, kind = NULL, = NULL)(1234,kind = "Mersenne-Twister", = "Inversion") # 默认参数设置runif(5)(1234,kind = "Mersenne-Twister", = "Inversion") # 默认参数设置runif(5)SWB发生器中“借位相减”步骤是指序列的第i个随机数zi要依据如下递推关系产生，zi=zi

6、+20zi+5b下标i,i+20,i+5都是对32取模的结果，zi+20与zi+5做减法运算，b是借位，其取值与前一步有关，当zi是正值时，下一步将b置为0，如果计算的zi是负值，在保存zi之前，将其值加1，并在下一步，将b的值设为 253。下面关于随机数生成的效率和后面的统计检验，都以生成224个为基准，是1600多万个，取这么多，一方面为了比较编程语言实现的发生器产生随机数的效率，另一方面是后面的游程检验需要比较大的样本量。Matlab内置的发生器及大部分的函数，底层实现都是 C 或者 Fortran，MathWorks创始人Cleve B. Moler是数值分析领域著名的LINPACK和

7、EISPACK包的作者之一，他当年做了很多优化和封装，进而推出Matlab，只要是调用内置函数，效率不会比C差，自己写的含有循环、判断等语句的代码，性能就因人而异了，对大多数人来说，性能要比C低。这里比较Matlab内置SWB发生器（就当作是C语言实现的好了和用Matlab重写的SWB发生器的效率,代码如下：% matlab codetic % 大约几秒rng(1234, 'v5uniform') % 调用SWB发生器x = rand(1,224);toc% octave codeid = tic % 时间耗费大约一小时randtx('state',0)x =

8、randtx(1,224);toc (id)randtx不是Matlab和Octave内置的函数，而是Cleve 基于Matlab实现的SWB发生器，也是100多行包含嵌套循环等语句的Matlab代码打包的函数，上面的代码运行时间差异之大也就不难理解了，为了能在Octave上跑，我做了少量修改，Octave软件版本为4.2.1，安装Octave时,Blas选择OpenBlas，为了后续检验，在获得随机数后，将其保存到磁盘文件 -mat x R，Octave，Matlab和Python内置的发生器都是MT发生器，与之实现有关的数学库，也是Blas，虽然有开源和进一步优化的商业版本之分，但是矩阵乘

9、法，向量乘法之类运算的效率也就半斤八两，Julia语言官网给出了一个标准测试2。不同语言的性能表现（C语言在算法中的表现为基准，时间记为1.0）这里再给出用C语言实现的MT发生器，产生同样多的随机数，只需要10秒左右，其中包含编译和写随机数到文件的时间，实际产生随机数的时间应该远小于这个时间。（程序运行环境环境Dev-C+ 5.11，用64位的GCC编译）。统计检验随机数的检验是有一套标准的，如George Marsaglia开发的DieHard检验程序，检验的内容很丰富。这篇文章只能算初窥门径，R内产生真随机数的包是 Dirk Eddelbuettel 开发的random包，它是连接产生真随

10、机数网站的接口。相关性检验先来一个简单的，就用R产生的两个独立同均匀分布样本，调用做相关性检验，看看这两组数是不是足够独立同分布，通过眼球验证，随着样本量增大，相关性趋于0，相关性弱到可以视为独立。如下图所示library(ggplot2)library(viridisLite)library(viridis)(1234)corr <- rep(0, 1000) for(i in seq(from = 1000, to = 1000000, by = 1000) corri/1000 <- (runif(i, min = 0, max = 1), runif(i, min = 0,

11、 max = 1)$estimate ggplot(x = seq(1000), y = corr), aes(x = x, y = y) + geom_hex( = FALSE) + scale_fill_viridis(direction = -1) + xlab("Sample size *103") + ylab("Correlation") 分布检验检验产生的随机数是否服从指定的分布：原假设是样本来自指定的分布，计算的P值比较大，就不能拒绝原假设。(runif(1000), "punif") #分布检验# One-sampl

12、e Kolmogorov-Smirnov test# data: runif(1000)# alternative hypothesis: two-sided检验两样本是否来自同一分布：原假设是两样本来自同一分布，计算的P值比较小，就表示两样本不是来自同一分布。est(runif(1000), runif(1000) #同分布检验# Two-sample Kolmogorov-Smirnov test# data: runif(1000) and runif(1000)# alternative hypothesis: two-sided从结果来看,R内置的随机数发生器通过了检验（嘿嘿，这是肯

13、定的！）。游程检验游程检验对随机数序列的随机性检验，不对序列做任何分布假设，不同于上面的相关性检验和省略没讲的特征统计量（如均值和方差）的检验。先对随机性略作解释，简单起见，我们考虑0-1序列，抛掷均匀的硬币1000次，正面向上记为1，反面向上记为0，这是一个离散的均匀分布，每一次抛掷硬币都无法准确地判断出现的是正面还是反面，若记录的序列中0和1相对集中的出现，不是随机，0和1交替出现，呈现周期性也不是随机，除了这两种情况基本就是随机了。游程检验的原假设是序列满足随机性，序列一旦生成，就是有序的，因为游程检验需要固定位置，再数上升（下降）的游程数，当计算的P值比较大时，不能拒绝原假设，即不能否

14、认这个序列是随机的。对上述0-1序列进行模拟，然后检验，如下所示library(tseries)x <- sample(c(0, 1), 1000, replace = TRUE, prob = c(1/2, 1/2)(factor(x)# Runs Test# data: factor(x)# alternative hypothesis: 现在用游程检验比较SWB发生器（Octave/Matlab版）、MT发生器（R语言版）和MT发生器（C语言版）。对于一般的实数序列，得先指定一个阈值，记为delta，然后比较序列中的值和delta的大小关系，这里类似数上升或下降的游程长度（runs

15、 length），基于这样一个事实：如果序列表现随机，则序列中两个小于delta的值，间隔越远出现的次数越少。可以这样理解，还是以上面抛硬币的例子来说，出现了很多次正面，那么下一次抛掷倾向于出现反面，这是一条人人可接受的经验。为了把这条经验可视化出来，对序列做如下操作：先给定阈值delta为0.01（也可以取别的值），取出序列中的值小于delta的位置，位置序列前面添加0，再差分，然后每个值减1，得到新序列，新序列中的值为0，表明原序列连续两个值小于delta，新序列中的值为1，表明间隔1个数小于delta，新序列中的值为2，表明间隔2个数小于delta，依次类推. 统计所有的情况，用直方图显

16、示，这就获得游程长度与间隔的关系图（间隔数取到100足可示意）。library(gridExtra)library()# 游程频数直方图run_test_fun <- function(x, string, delta) n <- length(x) len <- diff(c(0, which(x < delta), n + 1) - 1 ggplot(x = lenlen < 101), aes(x, fill = .count.) + scale_fill_viridis(direction = -1) + geom_histogram(binwidth =

17、 1, = FALSE) + xlab(string) + ylab("") (1234) # R默认采用Mersenne Twister发生器r_data <- runif(224, 0, 1); # R内生成均匀分布随机数matlabv5_data <- readMat("") # 读取Octave生成的均匀分布随机数temp <- (file = "random_number.txt") # 读取C语言生成的均匀分布随机数 c_data <- c(t(temp)p1 <- run_test_fun(

18、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)(p1, p2, p3, ncol=3)从图中可以看出MT发生器通过了检验，SWT发生器没有通过，在间隔数为27的位置，有一条沟，按规律游程长度不应该减少这么多，这是因为SWB发生器

19、“借位减”步骤，取模32的运算和5一起消耗了间隔为27的数量（读者可以按借位减的递推关系思考是如何消耗的），导致不符合随机性的要求，该算法细节参见Cleve B. Moler的书Numerical Computing with MATLAB第9章第267页 。应用两个均匀分布的统计模拟随机变量X1,X2iid某分布（比如二项分布，泊松分布，正态分布，指数分布，卡方分布，伽马分布），则X1+X2也服从该分布。常见的均匀分布是否具有这样的可加性？具体地说，就是X1,X2iidU(0,1),X1+X2是否服从U(0,2)？如果有一台电脑在旁边，我首先想到的就是敲三五行代码，画个散点图、直方图，看图说

20、话。(1234) x <- runif(10000, min = 0, max = 1) y <- runif(10000, min = 0, max = 1) z <- x + yplot(z) # 散点图hist(z) # 直方图为美观起见，从viridis包调用viridis调色板，颜色越深的地方，相应的数值越大，不管是此处geom_hex绘制的六角形热图，还是geom_histogram绘制的直方图，都遵循这个规律。ggplot(x = seq(10000), y = z), aes(x = x, y = y) + geom_hex( = FALSE) + scale

21、_fill_viridis(direction = -1) + xlab("") + ylab("")显然这不是均匀分布，在 z=1处，散点比较集中，看起来有点像正态分布。如果往中心极限定理上靠，将作如下标准化Y2=X1+X22121122=6(X1+X21) 则Y2的期望为0，方差为1。p4 <- ggplot(x = z), aes(x, fill = .count.) + scale_fill_viridis(direction = -1) + geom_histogram(bins=20, = FALSE) + xlab("&qu

22、ot;) + ylab("") p5 <- ggplot(x = sqrt(6)*(z-1), aes(x, fill = .count.) + scale_fill_viridis(direction = -1) + geom_histogram(bins = 20, = FALSE) + xlab("") + ylab("") (p4, p5, ncol=2)只是变换后的图像和之前基本一致，那么现在看来眼球检验不好使了，那就上P值呗！(sqrt(6)*(z-1), "pnorm") # 分布检验# One

23、-sample Kolmogorov-Smirnov test# data: sqrt(6) * (z - 1)# alternative hypothesis: two-sided也不是正态分布，既然如此，那就在两个随机变量的情况下，把精确分布推导出来。精确分布的推导及计算课本如概率论与数理统计教程 采用卷积的方法求分布函数，这种方法实行起来比较繁琐，也不利于后续编程，下面考虑用特征函数的方法求。我们知道标准均匀分布的特征函数(t)=eit1it考虑X1和X2相互独立，它们的和用S2表示，则随机变量S2的特征函数为2(t)=(t)(t)=(eit1it)2=2(1cos(t)eitt2只要满

24、足条件+|2(t)|dt<S2的密度函数就可以表示为p2(x)=12+eitx2(t)dt经计算+|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要说数值计算一个p(x)近似值，是一点问题没有！且看integrate(function(t,x

25、,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 error<6.6e-05那如果要把上面的积分积出来，获得一个精确的表达式，在n=2的时候还可以手动计算，主要使用分部积分，余弦积化和差公式和一个狄利克雷积分公式+0sin(ax)xdx=2sgn(a)，过程略，最后算得p2(x)=12(2x)sgn(2x)xsgn(x)(1x)sgn(1x)=12(|x|+|x2|)|x1|,0<x<2p2(x

26、)的密度函数图象如下：fun_p2_1 <- function(x) 1 / 2 * (abs(x - 2) - 2 * abs(x - 1) + abs(x) fun_p2_2 <- function(x) x <- (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(

27、x,tempfun) )ggplot(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") 从图中可以看出，两种形式的密度函数在数值计算的结果上很一致，当n=100,1000时，含参量积分的表示形式就很方便啦！任意给定一个n，符号计算上面的含参量积

28、分，这个时候还是用软件计算比较合适，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_

29、lift(a)*2, oo), 0), (Integral(sin(a*t)/t, (t, 0, oo), True)初次见到这样的结果，是不是一脸mb，翻译一下，就是2forperiodicargument(polar_lift2(a),)=001tsin(at)dtotherwise稍为好点，但是还是有一大块看不懂，那个绝对值里是什么 3？还是不要纠结了，路远坑多，慢走不送啊！话说要是计算p2(x)密度函数里的积分，from sympy import * x=symbols('x', real=True)t=symbols('t', real=True,po

30、sitive=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)2<10for4(2x2)2<1G3,14,4(121,