编写一个产生符合高斯分布的随机数函数.docx

信号检测与估计课程作业作业要求1、利用计算机内部函数产生高斯分布的随机数，分别画出500，10000，100000点的波形，并进行统计分析(分别画出概率密度曲线，计算均值与方差)2、利用计算机自己编写一个产生符合高斯分布的随机数函数，画出100000点的波形，并进行统计分析（同一）提示：这一问分两步做，第一步先产生一个均匀分布的随机数序列（乘同余法、混合同余法等，可以用自己的方法），第二步通过适当变换得到符合高斯分布概率模型的随机数列3、对随机数产生函数和高斯分布进行性能分析，并写出自己对于此次作业和上课的学习体会一、利用内部函数产生高斯分布首先利用matlab自带的内部函数randn()就可以方便的生成所需要的高斯分布随机数，然后画出概率密度曲线并计算出均值与方差即可。程序代码如下：A=randn(500,1);B=randn(10000,1);C=randn(100000,1);subplot(2,3,1);bar(A);subplot(2,3,2);bar(B);subplot(2,3,3);bar(C);f1,x1=ksdensity(A);subplot(2,3,4);plot(x1,f1);title(500点高斯分布概率密度函数);f2,x2=ksdensity(B);subplot(2,3,5);plot(x2,f2);title(10000点高斯分布概率密度函数);f3,x3=ksdensity(C);subplot(2,3,6);plot(x3,f3);title(100000点高斯分布概率密度函数);JZ500=mean(A)JZ1000=mean(B)JZ100000=mean(C)FC500=var(A)FC10000=var(B)FC100000=var(C)运行代码之后，可以得到如下结果：500点的均值为JZ500 =0.0334 10000点的均值为JZ1000 =0.0101 100000点的均值为JZ100000 =0.0016500点的方差为FC500 =0.9649 10000点的方差为FC10000 =1.0105 100000点的方差为FC100000 =0.9999我们可以看到随着实验点数的增加，生成的随机数的各项指标越来越接近于标准正态分布。二、编写随机数函数1、映射原理因为我以前曾经做过一个可以生成指数分布随机数的函数，所以这一次我的想法还是基于反函数的方法生成符合高斯分布的随机数。基于这种思想，首先要对标准高斯分布的概率密度函数和分布函数曲线进行分析。为了得到概率密度曲线和分布函数曲线的图像，可以在matlab中运行如下代码：x=-5:0.02:5;y=exp(-0.5*x.2)/(2*pi)0.5;plot(x,y),axis(-5,5,0,0.5)即可看到标准高斯分布的概率密度曲线图fx=12e-x22再运行如下代码：clearx=-10:0.02:10;y =(1125899906842624*2(1/2)*pi(1/2)*(erf(2(1/2)*x)/2)+1)/5644425081792261;plot(x,y),axis(-10,10,0,1.2),grid on可以得到标准高斯分布的分布函数曲线图像Fx=12-xe-t22dt可以看到F（x）具有如下的优良性质：（1） 单调递增，理论上一定存在反函数；（2） 0F（x）1，且F（x）落在区间（a，b）的概率为b-a，所以令y=F（x），则y可以看成为一个满足（0,1）上的均匀分布的随机变量。正因为分布函数具有如上两个性质，所以在理论上一定可以得到F（x）的反函数，并用这个反函数可以将一个（0,1）上的均匀分布映射成一个标准高斯分布。但是，与指数分布、瑞利分布等其他分布不同，高斯分布的分布函数：Fx=12-xe-t22dt不能直接积分出来，也即不能用有限的解析形式来表示，也就是说虽然F（x）的反函数一定存在但是却写不出来，这给我们的映射带来了困难。但是我们可以考虑利用二维的正态分布来解决这一问题。假定r1与r2是0,1区间的两个独立的均匀分布随机数，现将其作如下的变换，令x1=-2ln(r1)*cos(2*r2)x2=-2ln(r1)*sin(2*r2)此时再将上述两个式子做反变换，可以得到：r1=e-x12+x222r2=12arctanx2x1+c其中c为常数，由此可导出其密度函数为：fx1,x2=12e-x12+x222由这个二维正态分布的概率密度可以知道，x1与x2是两个不相关的标准正态分布:（x1,x2）N(0,0,1,1,0)也即两个相互独立的高斯分布。至此，我们找到了将均匀分布映射成为高斯分布的解析表示式。2、用乘同余法生成均匀分布找到了映射表达式之后，我们的任务就是寻找可以生成均匀分布随机数的方法。这里用到的方法是乘同余法。乘同余法的迭代式如下： Xn+1=Lamda*Xn(mod M) Rn+1=Xn/M当然，这里的参数选取是有一定理论基础的，否则所产生的随机数的周期将较小，相关性会较大。经过前人检验的两组性能较好的素数取模乘同余法迭代式的系数为: Lamda=55,M=235-31 Lamda=75,M=231-13、程序代码因为我们需要产生两个相互独立的均匀分布，为了使相关性更低，我们产生的两个均匀分布随机变量r1、r2分别使用不同的迭代系数。并且为了增大随机型，决定不使用固定种子，而是通过提取系统时间(当前时间的秒数)产生一个时间随机种子。运行代码如下：a=clock; x1(1)=a(6)/100; %设置一个时间种子%for i=1:1:100000 %生成100000个随机数%x1(i+1)=mod(75*x1(i),231-1); %采用第二组系数Lamda=75,M=231-1%r1(i)=x1(i)/(231-1);endb=clock;x2(1)=b(6)/100;for i=1:1:100000x2(i+1)=mod(55*x2(i),235-31); %采用第一组系数Lamda=55,M=235-31%r2(i)=x2(i)/(235-31);endx=sqrt(-2*log(r1).*sin(2*pi*r2); %将均匀分布映射成高斯分布%subplot(1,2,1);bar(x);title(100000点的波形);f,xi=ksdensity(x);subplot(1,2,2);plot(xi,f);title(100000点的高斯分布概率密度);JZ=mean(x) %求均值%FC=var(x) %求方差%可以得到如下结果：平均值为JZ =0.0011方差为FC =0.9941三、性能分析和学习体会1、性能对比分析（1）用matlab内部函数产生的随机数随着生成点数的增加（从500到100000）其统计特性越来越接近于标准高斯分布，当点数达到100000点时，其概率密度曲线已经十分接近标准概率密度曲线，并且均值与方差也打到了非常高的接近度（均值为0.0016，方差为0.9999）（2）利用自己编写的函数所产生的100000点随机数在统计特性上也十分接近标准的概率密度，其概率密度曲线和均值（0.0011）的接近程度与matlab内