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

信号检测与估计课程作业

作业要求

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.0016

500点的方差为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])

即可看到标准高斯分布的概率密度曲线图

再运行如下代码：

clear

x=-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

可以得到标准高斯分布的分布函数曲线图像

可以看到F（x）具有如下的优良性质：

（1）单调递增，理论上一定存在反函数；

（2）0≤F（x）≤1，且F（x）落在区间（a，b）的概率为b-a，所以令y=F（x），则y 可以看成为一个满足（0,1）上的均匀分布的随机变量。

正因为分布函数具有如上两个性质，所以在理论上一定可以得到F（x）的反函数，并用这个反函数可以将一个（0,1）上的均匀分布映射成一个标准高斯分布。

但是，与指数分布、瑞利分布等其他分布不同，高斯分布的分布函数：

不能直接积分出来，也即不能用有限的解析形式来表示，也就是说虽然F（x）的反函数一定存在但是却写不出来，这给我们的映射带来了困难。但是我们可以考虑利用二维的正态分布来解决这一问题。

假定r1与r2是［0,1］区间的两个独立的均匀分布随机数，现将其作如下的变换，令

此时再将上述两个式子做反变换，可以得到：

其中c为常数，由此可导出其密度函数为：

由这个二维正态分布的概率密度可以知道，x1与x2是两个不相关的标准正态分布:

（）

也即两个相互独立的高斯分布。至此，我们找到了将均匀分布映射成为高斯分布的解析表示式。

2、用乘同余法生成均匀分布

找到了映射表达式之后，我们的任务就是寻找可以生成均匀分布随机数的方法。这里用到的方法是乘同余法。

乘同余法的迭代式如下：

Xn+1=Lamda*Xn(mod M)

Rn+1=Xn/M

当然，这里的参数选取是有一定理论基础的，否则所产生的随机数的周期将较小，相关性会较大。经过前人检验的两组性能较好的素数取模乘同余法迭代式的系数为:

Lamda=5^5,M=2^35-31

Lamda=7^5,M=2^31-1

3、程序代码

因为我们需要产生两个相互独立的均匀分布，为了使相关性更低，我们产生的两个均匀分布随机变量r1、r2分别使用不同的迭代系数。并且为了增大随机型，决定不使用固定种子，而是通过提取系统时间(当前时间的秒数)产生一个时间随机种子。

运行代码如下：

a=clock;

x1(1)=a(6)/100; %设置一个时间种子%

for i=1:1:100000 %生成100000个随机数%