随机数产生原理及实现

电子信息与通信工程学院
实验报告
实验名称随机数的产生
课程名称随机信号分析
姓名顾康学号U201413323 日期6月6日地点南一楼东204 成绩教师董燕
以上为6种分布的实验结果
1.均匀分布
随机变量X~U(0，1)的一组样本值的模拟值一般采用某种数值计算方法产生随机数序列，在计算机上运算来得到，通常是利用递推公式:
Xn=f(Xn-1,.....,Xn-k)
1.1 同余法
Xn+1 = λXn(mod M)
Rn=Xn/M
R1 R2...Rn即为（0,1）上均匀分布的随机数列。

而上述方法是伪随机的，{Rn}本质上是递推公式给定的周期序列，周期T可看做logλ（M）。

解决方法是：选择模拟参数并对序列进行统计检验。

1.2选择模拟参数
1）周期长度取决于Xo,λ， M的选择
2）通过选取适当的参数可以改善随机数的性质
几组参考的取值
Xo =1 , λ=7 , M=10^10
Xo =1 , λ=5^13 , M=2 *10^10
Xo =1 , λ=5^17 , M=10^12
1.3对数列进行统计检验
对应序列能否看作X的独立同分布样本，须检验其独立性和均匀性
for i=2:1:size %同余法均匀分布
x(i)= mod ( v*x(i-1), M);
y(i)=x(i)/M;
end
subplot(2,3,1);
hist(y,100)
[ahat,bhat,ACI,BCI]=unifit(y)% 以0.95的置信度估计样本的参数
首先我们的标准是U ~（0，1），而实验值，ACI表示ahat的范围[-0.0030,0], BCI表示bhat的范围[1.0000,1.0030]。

同时样本的均值和方差分别为0.4932
和0.0830，结论与理论值很接近。

该样本以0.95的可信度服从（0,1）均匀分布。

2.伯努利分布
2.1算法原理
若随机变量R服从（0,1），
P(X=Xi)=Pi
P(0)=0, P(n)=∑Pi
P{P(n-1)<R<=P(n)}=P(n)-P(n-1)=Pn
令{P(n-1)<X<=P(n)}={X=Xn} 有P(X=Xn)=Pn
从理论上讲，已经解决了产生具有任何离散型随机分布的问题。

具体执行仍有困难，如X取值为无穷时。

2.2算法
对于伯努利分布只需用到上述算法最简单的情形，即取n为2。

%伯努利分布
k1=0,k2=1
p1=0.2,p2=0.8;
r=zeros(1,size);
for j=1:1:size
if y(j)>p1
r(j)=k2;
else
r(j)=k1;
end
end
subplot(2,3,2);
hist(r)
title('伯努利分布');
[PHAT,PCI]=binofit(r,1000)%以0.95的置信度检验样本参数
PHAT=0.198，而PCI= [0.195,0.212 ] ，而我设置的P=0.2，与实验结果十分接近，可见该样本的性质较好。

该样本以0.95的可信度服从0.2的伯努利分布。

3.正态分布
3.1算法
设有两个在（0,1）上独立均匀分布的随机数R1 ,R1;作如下变换
Y1= （-2㏑R1）½COS(2πR2）
Y2= （-2㏑R1）½SIN(2πR2）
其逆变换为
R1=exp(- (Y1²+Y2²)/2)
R2=1/2πarctag(Y2/Y1)
可导出Y1，Y2的联合分布函数
f(Y1,Y2)= 1/2πexp(- (Y1²+Y2²)/2)故知Y1，Y2相互独立且服从N（0,1）再作变换Xi = σYi+μ,可得到服从N(μ，σ)的X
3.2代码
%正态随机分布
y=rand(1,size);
z=rand(1,size);
m=sqrt(-2*log(y)).*cos(2*pi.*z);
n=sqrt(-2*log(y)).*sin(2*pi.*z);
t=[m,n];
subplot(2,3,3);
hist(t,100)
title('正态分布');
[muhat,sigmahat,muci,sigmaci]=normfit(t)%以0.95置信度检验样本参数
实验结果：muhat=0.0060, sigmahat=1.0080;
Muci=[-0.0382,0.0502] ,sigmaci=[0.9777,1.0402];
理论上假设的是标准正态分布，则该样本以0.95的可信度服从标准正态分布。

4.指数分布
4.1算法
已知指数分布的分布函数为F=1-exp(-λx),F∈[0,1]。

而利用服从U(0，1)的变量Y, 带入方程左边反解出x：
X= -log(1-y)/λ
即我们得到了用均匀分布产生指数分布的方法
4.2代码
%指数分布
y=rand(1,size);
lamda=3;
x=-log(1-y)/lamda;
subplot(2,3,4);
hist(x,30)
title('指数分布')
[muhat,muci]=expfit(x)%以0.95的置信度检验样本参数
理论上λ取为3，而实验值muhat=0.3204, muci=[0.3015,0.3402] ,可得到该样本以0.95的可信度服从E(3).
5.泊松分布
5.1算法
利用一组在（0,1）独立且同均匀分布的变量X，首先理解泊松分布的到达间隔服从指数分布：
在一次独立实验中∏X k≤exp(-λ)作为判断成功条件，输出为k，即X的次数。

5.2代码
%泊松分布
lamda=20;
p=exp(-lamda);
y=zeros(1,100);
for cnt=1:1:100
i=0;
q=1;
while q>=p
q=q*rand(1);
i=i+1;
end
y(cnt)=i-1;
end
subplot(2,3,5);
hist(y,25)
title('泊松分布')
[lamdahat,lamdaci]=possifit(y)%以0.95的置信度检验样本参数
注意到lamda初值被赋为20，实验值lamdahat=20.1900,
Lamdaci=[19.3093,21.0707].所以样本以0.95的可信度服从P(20).
6.几何分布
6.1算法原理
先由几何分布的定义出发P(X=k)= q^(k-1) * p;则可利用与产生泊松分布随机数相似的方法，从均匀分布入手。

产生一组（0,1）独立地均匀分布随机数，分别与P比较。

如果Xi<P,则继续X（i+1）；但如果Xk≥P,则输出K（代表比较的次数）
6.2代码
%几何分布
y=zeros(1,100);
p=0.6;
q=0.4;
test=0;
cnt=1;
while cnt<=100
if test==1
x=rand(1,size);
end
for i=1:1:size
if x(i)>p
y(cnt)=i;
cnt=cnt+1;
test=1;
break
end
end
end
subplot(2,3,6);
hist(y,20)
title('几何分布')
[phat,pci]=mle(y,'distribution','geometric')%极大似然估计函数，以0.95置信度预测p
我给P的初值为0.6，但是值得注意的现象是，多次运算得到的phat都在0.3左右。

当改变P值后，phat仍然会分布在另一个值附近。

通过pci的范围我们可以判断样本以0.95的可信度属于几何分布，但概率P并不与初值一致。

我认为这可能是算法导致的，并无法保证初值P仍然是样本的概率。

7.总结
在这次实验中，最大的收获是体会到均匀分布随机数是很多其它类型随机数的基础，所以均匀分布随机数的质量十分重要。

不仅是产生随机数，判断随机数质量的过程加深了我对所谓伪随机算法的理解，我们需要不断完善伪随机的算法以达到近似随机的效果。

8.参考文献
1... 《MATLAB统计分析与应用》北京航空航天大学出版社。