CMRG 生成器的实现和检验分析
将两个或多个单独的发生器以某种方式结合起来生成最终的随机数，这种组合的发生器较之组成它的任何单个发生器有更长的周期，更好的统计性能。采用组合发生器的缺点在于，得到每个
本文着重分析利用matlab实现L’ Eculyer(1999）提出的CMRG成器实现方法和matlab rng源设为‘combRecursive’生成随机序列的方法的性能。
为方便描述，下文将用方式一代替利用matlab实现L’ Eculyer(1999）提出的CMRG成器实现方法。
用方式二代替matlab rng源设为‘combRecursive’生成随机序列的方法。
N表示生成随机变量个数。
一、CMRG生成器算法原理
简单来说，其思想时用{$Z{1i}$}和{$Z{2i}
例如当
二、检验分析
1、均匀性检验
下方两张图片分别通过方式一和方式二生成30000个随机变量的统计图。从直方图和散点图观察来看，这两种方式均满足分布的均匀性。
通过
chi2_test =
'chi2test pass'2、序列检验
序列检测实际上就是把
两种方式在N=30000，90000，3000000时均通过序列检测，出入如下：
chi2_ser =
'serial test pass'
3、游程检验
当N=30000时
r_ij={4921 6103 2873 777 178 41} %方式一
r_ij={5014 6180 2730 813 177 48} %方式二当N=90000时
r_ij={14933 18599 8329 2369 511 137} %方式一
r_ij={15231 18671 8194 2371 525 118} %方式二当N=3000000时
r_ij={50104 62353 27383 7987 1733 392} %方式一
r_ij={50193 62616 27324 7904 1743 367} %方式二两种方式均通过上升游程检验，运行结果如下：
chi2_run =
'runs up test pass'4、相关性检验
当N=30000时，方式一和方式二都能通过检验
chi2_cor =
'Correlation test pass'当N=90000时，同上。
当N=300000时，结果同上。
两种方式均通过相关性检验。
5、生成时间对比
利用方式一生成30000个随机变量所用时间为
t_code =
0.2520使用方式二生成30000个随机变量所用的时间为
t_rng =
0.0100当N=300000时，
t_code =
40.9300
t_rng =
0.0130如上述所用时间所示，在生成相同数量随机变量时，方式二所用时间更少。
三、结论
方式一和方式二有着近似的分布特征，但随着需要生成的随机变量的数量增加，方式二所需时间变化不大，而方式一生成随机序列所需时间呈指数增长。
四、程序源码
clear ;
close ;
%生成随机数的个数
N = 100000*3;
%两种方式生成
type = 2;
%m1,m2
m1=2.^32-209;m2=2.^32-22853;
%z1系数
a1=1403580;a2=810728;
%z2系数
b1=527612;b2=1370589;
%z初值
z1(1)=2637242579;z1(2)=1747669293;z1(3)=2112977067;
z2(1)=2633436652;z2(2)=3518603782;z2(3)=229738972;
%求随机数out
if type==1
t1=clock;
for i=4:(N+3)
z1(i)=mod(a1*z1(i-2)-a2*z1(i-3),m1);
z2(i)=mod(b1*z2(i-1)-b2*z2(i-3),m2);
y(i-3)=mod(z1(i)-z2(i),m1);
out=y./m1;
end
t2=clock;
t_code=etime(t2,t1)
else
%matlab rng源生成
t1=clock;
rng('shuffle','combRecursive');
out = rand(1,N);
t2=clock;
t_rng=etime(t2,t1)
end
%可视化
plot(out,'.')
figure
histogram(out,10)
figure
%相关性检验
n=N/3;
for i=1:n
X(i)=out((i-1)*3+1);
Y(i)=out((i-1)*3+2);
Z(i)=out((i-1)*3+3);
end
k=ceil(n.^(1/3));
Hist=zeros(k,k,k);
for i=1:n
j1=ceil(X(i)*k);
j2=ceil(Y(i)*k);
j3=ceil(Z(i)*k);
Hist(j1,j2,j3)=Hist(j1,j2,j3)+1;
end
for i=1:10
for j=1:10
for m=1:10
h(i,j,m)=((Hist(i,j,m)-n/(k^3)).^2);
end
end
end
kf_s=sum(sum(sum(h))).*(k^3)/n;
kf = chi2inv(0.95,k^3-1);
if kf_s<kf
chi2_cor='Correlation test pass'
else
chi2_cor='Correlation test fail'
end
%chi2test
k=10;
n=N;
n1=hist(out,k);
kf_s = k/n*(sum((n1-n/k).^2)); % calculate the statistic
kf = chi2inv(0.95,k-1);
if kf_s<kf
chi2_test ='chi2test pass'
else
chi2_test ='chi2test fail'
end
% cdftest
chi2_p=chi2cdf(kf_s,k-1) ;
if chi2_p<0.95
chi2_cdf='chi2cdf pass'
else
chi2_cdf='chi2cdf fail'
end
%序列检验
for i=1:N/2
x1(i)=out((i-1)*2+1);
x2(i)=out((i-1)*2+2);
end
k=100;
Hist=zeros(k);
for i=1:length(x1)
j1=ceil(x1(i)*k);
j2=ceil(x2(i)*k);
Hist(j1,j2)=Hist(j1,j2)+1;
end
bar3(Hist)
for i=1:k
for j=1:k
h(i,j)=((Hist(i,j)-(n/2)/(k^2)).^2);
end
end
kf_s = sum(sum(h)).*(k^2)/(n/2);
kf = chi2inv(0.95,k^2-1);
if kf_s<kf
chi2_ser='serial test pass'
else
chi2_ser='serial test fail'
end
%游程检验
a_ij=[4529.4 9044.9 13568 18091 22615 27892
9044.9 18097 27139 36187 45234 55789
13568 27139 40721 54281 67852 83685
18091 36187 54281 72414 90470 111580
22615 45234 67852 90470 113262 139476
27892 55789 83685 111580 139476 172860];
b_i=[1/6,5/24,11/120,19/720,29/5040,1/840];
w=[1,2,3,4,5,6];
n=N;
x=out;
flag=0;
r_ij=zeros(1,6);
temp=1;
for i=1:length(x)-1
if x(i+1)>x(i)
temp=temp+1;
else
if temp<6
r_ij(temp) = r_ij(temp) + 1;
else
r_ij(6) = r_ij(6) + 1;
end
temp=1;
end
end
for i=1:6
for j=1:6
s(i,j)=a_ij(i,j)*(r_ij(i)-n*b_i(i))*(r_ij(j)-n*b_i(j));
end
end
kf=sum(sum(s))/n;
chi2_p=chi2cdf(6,kf);
if chi2_p<0.95
chi2_run='runs up test pass'
else
chi2_run='runs up test fail'
end
