迭代取中法，线性同余法及混合同余法产生伪随机数方法的scilab实现

在很多计算机仿真算法及数字信号处理算法中，经常会用到随机数序列，比如进行蒙特卡洛仿真时就会用到。而且产生的随机数序列质量会对仿真结果产生比较大的影响，因此如何在仿真中产生随机数也是一个值得研究的课题。事实上，通过计算机编程不可能实现真正的随机序列，因为产生的序列达到一定的长度后，会出现重复序列，使用的随机序列长度不应大于一个周期。这种随机数只能称作伪随机数。本文介绍常用的伪随机数发生器-迭代取中法，线性同余法和混合同余法在scilab下的实现。

迭代取中法，乘同余法及混合同余法的迭代公式如下：





三种方法中，第二个公式均是为了得到[0,1]上的值，真正取随机迭代作用的是第一个式子。三种不同的方法各有特点。

1）迭代取中法容易退化成0。输入参数主要有种子，随机数长度，其scilab代码为:

function x=rand1(r0,n,s)
r=zeros(1,n+1);
x=zeros(1,n);
r(1)=r0;
for i=2:n+1
r(i)=pmodulo((r(i-1)^2)/10^s,10^(2*s)); //迭代
x(i-1)=r(i)/10^(2*s);                                   //取值范围为[0,1]
end
endfunction

function testofrand1
x=rand1(4532,500,2);               //产生500个数据，种子为4532,，S=2
n=1:size(x,2);
clf();
subplot(211);
plot2d(n,x,style=2);                    //画图
subplot(212);
histplot(10,x,style=3,normalization=%t); //画结果直方图，分为10类并归一化
endfunction

运行后，可以得到如下的结果





可以看出，随着迭代的进行，350个数据之后基本上已经退化为0了，这也可以从直方图中看出来，前面的序列均匀性还是比较好的。

2）乘同余的周期与lambda和M密切相关，根据查到资料，有两组会效果比较好:
Lambda=5^5,M=2^35-31
Lambda=7^5,M=2^31-1

其实现及测试代码如下：

function x=rand2(r0,n,lambda,M)
r=zeros(1,n+1);
x=zeros(1,n);
r(1)=r0;
for i=2:n+1
r(i)=pmodulo(lambda*r(i-1),M);
x(i-1)=r(i)/M;
end
endfunction

function testofrand2
x=rand2(5.0,500,5^5,2^35-31);
clf();
subplot(211);
plot2d(x);
subplot(212);
histplot(10,x,style=2);
endfunction

当Lambda=5^5,M=2^35-31，且序列长度为500时，效果如下





当Lambda=5^5,M=2^35-31，且序列长度为1000时，效果如下





可以看出均匀分布性还是比较好的，只是在个别段上会高一些，而且周期比迭代取中法长多了

3）混合同余法在如下条件下效果会比较好：
M=2^q
Lambda=2^c+1,c取q/2附近的值
Mu=(1/2+sqrt(3))/M
X0为任意非负整数

其scilab实现及测试代码如下：

function x=rand3(r0,n,lambda,Mu,M)
r=zeros(1,n+1);
r(1)=r0;
x=zeros(1,n);
for i=2:n+1
r(i)=pmodulo(lambda*r(i-1)+Mu,M);
x(i-1)=r(i)/M;
end
endfunction
function testofrand3
M=2^10;
lambda=2^4+1;
Mu=(1/2+sqrt(3))/M;
x=rand3(5,10000,lambda,Mu,M);
clf();
subplot(211);
plot2d(x);
subplot(212);
histplot(10,x,style=2);
endfunction

生成长度为10000的随机序列效果如下：





可以看出，整体效果都很好，没有出现退化现象，有资料上说在某个周期内会出现线性增长，可是在这个结果中没有出现。