建模与仿真作业-机数生成.docVIP

建模与仿真作业 随机数生成 本次作业使用MATLAB编写程序。 线性同余法 基本公式： 其中a为乘子，为种子，c为常数，M为模 迭代步骤： （1）设定a, c, M值，并给定初始种子 （2）令 （3） （4） （5），转入（3） 本题采用著名的 Coveyou与 Macpherson 混合同余发生器 及Kobayashi 混合同余发生器 程序代码： 线性同余法函数： function r = MixMOD(x0,n,type) %随机数种子：x0 %产生的随机数个数：n %采用的混合同余法的公式类型：type %产生的随机数序列：r format long; M1 = power(2,31); M2 = power(2,35); a1 = 314159269; a2 = power(5,15); c1 = 453806245; c2 = 1; r = zeros(n,1); x = zeros(n+1,1); x(1) = x0; if type == 1 for i=2:n+1 y = a1*x(i-1)+c1; x(i) = mod(y, M1); r(i-1) = x(i)/M1; end else for i=2:n+1 y = a2*x(i-1)+c2; x(i) = mod(y, M2); r(i-1) = x(i)/M2; end end format short; 对（0，1)均匀分布随机数进行矩检验函数： function [s1,s2,s3,s4]=moment_test(R) %对（0，1)均匀分布随机数进行矩检验 n=length(R); R_mean=mean(R);R_var=var(R);R_std=std(R); u1=sqrt(12*n)*(R_mean-0.5); if abs(u1)1.96 s1=pass; else s1=*; end % 对方差进行检验 var(R) u2=sqrt(180*n)*(R_var-1/12) if abs(u2)1.96 s2=pass; else s2=*; end % 对偏度进行检验 u3=mean(((R-R_mean)/R_std).^3)*0.408248*sqrt(n); if abs(u3)1.96 s3=pass; else s3=*; end % 对峰度进行检验 uu=mean(((R-0.5)/sqrt(1/12)).^4)-1.75 u4=uu*0.204124*sqrt(n); if abs(u4)1.96 s4=pass; else s4=*; end 相关系数检验函数： function [sacf1,sacf2,sacf3,sacf4,sacf5,sacf6,sacf7]=acf_test(R) % 独立性的自相关AFC检验 R_mean=mean(R);R_var=var(R);n=length(R); for i=1:7 rou(i)=sum(((R(1:n-i).*R(i+1:n)-R_mean^2))/R_var)*sqrt(1/(n-i)); end rou if abs(rou(1))1.96 sacf1=pass; else sacf1=*; end if abs(rou(2))1.96 sacf2=pass; else sacf2=*; end 线性同余法和移位寄存器法联合使用： function mt=initializeGenerator(seed) %初始化随机序列数组的函数 N=624; mt=zeros(1,N); mt(1)=seed; for i=2:N y=1812433253*bitxor(mt(i-1),bitshift(mt(i-1),-30))+i-1; mt(i)=uint32(bitand(uint64(y),hex2dec(ffffffff))); end; function [y2,mt1]=extractNumber(mt1,index) %根据index的值提取数组中的某个数来生成随机数 if index==1 mt1=generateNumbers(mt1); y1=mt1(index); y1=bitxor(y1,bitshift(y1,-11)); y1=bitxor(y1,bitand(bitsh