随机数的产生

lovehaboy

normal.m
%产生n个N(a,b)正态分布随机数
%其中a为均值，b为方差
%function x=normal(a,b,n)
function x=normal(a,b,n)
m=48;%
for i=1:n
    r=rand(1,m);
    x(i)=a+sqrt(b)*(sum(r)-m/2)/sqrt(m/12);
end

uinform.m
%产生一组[0，1]分布的随机数，采用余数法
%从文献中得知下列参数组较为适用：
%y的初值为1，w=7,模M=10^10(伪随机数周期为5*10^7)
%y的初值为1，w=5^13,模M=2^36(伪随机数周期为2^34,约2*10^10)
%y的初值为1，w=5^17,模M=2^42(伪随机数周期为2^40,约10^12)
%--------------------------------------------------------
%function x=uniform(y,w,M,n)
%n为要产生的随机数个数
function x=uniform(y,w,M,n)
for i=1:n
    x(i)=y./M;
    y=rem(w*y,M);%rem()为求余函数
end

riddling
%产生一个N(0,1)正态分布随机数
%采用筛选法，精度较高
%其他方法可参阅《现代应用数学手册--概率统计与随机过程卷》清华大学出版社马振华主编
%function y=riddling()
function y=riddling()
sign=0;
while 1
    x=rand(1,2);
    v1=2*x(1)-1;v2=2*x(2)-1;
    s=v1^2+v2^2;
    if s<=1
        a=sqrt(-2*log(s)/s);
        y=v1*a;
        sign=1;
    end
    if sign==1
        break;
    end
end