系统仿真（二）：CMRG与combRecursive两随机数生成器性能检验与分析

Sun-wz

已于 2022-03-30 23:28:21 修改

阅读量976

点赞数 1

分类专栏：仿真模拟文章标签： matlab

于 2022-03-29 16:01:33 首次发布

本文链接：https://blog.csdn.net/qq_36318271/article/details/123812426

版权

仿真模拟专栏收录该内容

5 篇文章 4 订阅

订阅专栏

一、题目背景：

将多个并行递归序列生成器结合可以有效实现具有长周期和良好结构特性的随机数生成器。这样的生成器具有很强的鲁棒性，L’ Eculyer(1999）提出的CMRG生成器实现相关方法，即：

$\begin{array}{l} Z_{1, \bar{l}}=\left(1403580 Z_{1, i-2}-810728 Z_{1, i-3}\right)\left[\bmod \left(2^{32}-209\right)\right] \\ Z_{2, i}=\left(527612 Z_{2, i-1}-1370589 Z_{2, i-3}\right)\left[\bmod \left(2^{32}-22853\right)\right] \\ Y_{\bar{i}}=\left(Z_{1, i}-Z_{2 i j}\right)\left[\bmod \left(2^{32}-209\right)\right] \\ U_{i}=\frac{Y_{\dot{i}}}{2^{32}-209} \end{array}$
本文将CMRG和matlab中RNG源为‘combRecursive’生成随机数两种方式在消耗时间、相关性检验、均匀性检验、序列检验和游程检验等五种方法作对比。

二、检验方法对比

消耗时间对比

表1 CMRG和combRecursive消耗时间对比表（单位:s）

数据量	3000	30000	300000	3000000
CMRG	0.0030	0.0110	0.1220	1.1270
combRecursive	0.0110	0.0110	0.0200	0.0130

相关性检验
用于直接计算随机数序列在 j 延迟后的相关性。经测试，在数据量分别为3000和30000时，CMRG和combRecursive的相关性检验均通过。
均匀性检验
对比两种方式生成30000个数据的散点图和直方图可知基本分布均匀。

图2.1 CMRG生成柱状图和combRecursive生成柱状图

图2.2CMRG生成散点图和combRecursive生成散点图通过卡方检验验证观测值与理论值之间的偏离程度

$\chi^{2}=\sum_{j=1}^{k}\left(n_{j}-N / K\right)^{2} /(N / K)$

检验结果：
在这里插入图片描述

④ 序列检验
视为卡方检验的高维版本，用于检验随机数生成器利用（u1,u2）(u1+1,u2+2)生成坐标或数据对下高维均匀性，具有独立意义。CMRG和combRecursive的相关性检验均通过。图3.1和图3.2是两种方法在数据量为30000下的图表。
在这里插入图片描述

图3.1 CMRG法柱状图

在这里插入图片描述

图3.2 combRecursive生成柱状图检验结果：

在这里插入图片描述

⑤ 游程检验
游程检验是根据样本标志表现排列所形成的游程的多少进行判断的检验方法。需统计随机序列的单调递增长度为1，2，3，4，5，和>=6的数目。构建统计量：

$\chi^{2}=\sum_{j=1}^{k}\left(n_{j}-N / K\right)^{2} /(N / K)R=\frac{1}{n} \sum_{i=1}^{6} \sum_{j=1}^{6} a_{i j}\left(r_{j}-n b_{i}\right)\left(r_{j}-n b_{j}\right)$

CMRG和combRecursive的游程检验均通过
在这里插入图片描述

三、总结

针对两种CMRG和matlab中RNG源为‘combRecursive’生成随机序列的方法不同检验方法均通过，不同点在于产生不同数据量大小的时间消耗上不同，combRecursive时间消耗总体上变化不大，而CMRG耗费时间随数据量的增加呈指数型增长。

附录：

clear;clc;close all
randomNumber = 3000; %100000
%函数参数
a1 = 1403580;a2 = 810728;a3 = 2.^32-209;
b1 = 527612;b2 = 1370589;b3 = 2.^32-22853;
z(1,1) = round(rand*1000000);z(1,2) = round(rand*1000000);z1(1,3) = round(rand*1000000);
z(2,1) = round(rand*1000000);z(2,2) = round(rand*1000000);z(2,3) = round(rand*1000000);

n = 1;
switch n
    case 1
        timeBegin = clock;
        for i=4:randomNumber+3
            z(1,i)=mod(a1*z(1,i-2)-a2*z(1,i-3),a3);
            z(2,i)=mod(b1*z(2,i-1)-b2*z(2,i-3),b3);
            y(i)=mod(z(1,i)-z(2,i),a3);%i-3
            U(i-3)=y(i)/a3;
        end
        timeEnd=clock;
        spendTime=etime(timeEnd,timeBegin)
    case 2
        timeBegin=clock;
        rng('shuffle','combRecursive');
        U = rand(1,randomNumber);
        timeend=clock;
        t_rng=etime(timeend,timeBegin)
end
%画图
plot(U,'.')
figure
histogram(U,10)
figure


%*************************相关性检验***************************************
n = randomNumber/3;
for i=1:n
    X(i)=U((i-1)*3+1);
    Y(i)=U((i-1)*3+2);
    Z(i)=U((i-1)*3+3);
end

k=ceil(n.^(1/3));
Hist=zeros(k,k,k);
for i=1:n
%     x1(i)
%     x2(i)
    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
%n=n/3;
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_str1 = '相关性检验通过'
else
    chi2_str1 = '相关性检验不通过'
end

%*********************************均匀性检验*******************************
k=10; %k个等长区间
n=randomNumber;
n1=hist(U,10);%统计计算在每个子区间上的随机数的个数
kf_s = (sum((n1-n/10).^2))*(k/n);  % 第五步
kf = chi2inv(0.95,k-1);%    卡方分布逆累积分布函数
if kf_s<kf
    chi2_str1 ='卡方检验通过'
else%拒绝H0,即不通过
    chi2_str1 ='卡方检验不通过'
end

%********************************序列检验**********************************
for i=1:1:length(U)/2
   x1(i)=U((i-1)*2+1);
   x2(i)=U((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)
n=randomNumber;
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_str1 = '序列检验通过'
else
    chi2_str1 = '序列检验不通过'
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=randomNumber;
x=U;
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
sum(r_ij.*w);

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_str1='游程检验通过'
else
    chi2_str1='游程检验不通过'
end