系统建模与仿真设计报告一

设计一产生十种不同分布的独立的随机数

一、设计内容及要求

任务：产生十种不同分布的独立的随机数，并进行检验。

要求：对随机数进行的统计性检验包括频率检验、参数检验、独立性检验。

二、设计环境及工具

Windows7、MatlabR2010b

三、设计思想及方法

(1) 在对雷达系统进行仿真时，首当其冲的问题就是对电磁环境

的仿、真。其中无用的电磁信号包括三大类，即杂波、噪声和干扰，在模拟仿真时相比于有用的电磁信号也是不可或缺的。其所谓的仿真就是在已知随机变量的统计特性及其参数的情况下，研究如何在计算机上产生服从给定统计特性和参数的随机变量。

(2) 在雷达、导航、声呐、通信和电子对抗等系统中，应用最多

的概率统计模型还是正态分布或高斯分布、指数分布、瑞利分布、莱斯分布或广义瑞利分布、韦尔分布、对数-正态分布、m分布、拉普拉斯分布、复合k分布等。

(3) 在这些随机总体中畸形随机抽样，实际上都是以[0,1]区间上

的均匀分布随机总体为基础的。原则上讲，只要已知[0,1]区间上的均匀分布随机数序列，总可以通过某种方法(数学方法)来获得某已知分布的简单子样。只要给定的均匀分布随机数列满足均匀

且相互独立打的要求，经过严格的数学变换或者严格的数学方法，所产生的任何分布的简单子样都会满足具有相同总体分布和相互独立的要求。 四、设计过程及结果

本次设计的十种随机数包括均匀分布、高斯分布、指数分布、广

义指数分布、瑞利分布、广义瑞利分布、韦尔分布、拉普拉斯分布、柯西分布和2χ分布，使用Matlab 完成设计并给出具体的参数，代码附在最后。 1．均匀分布

已知随机变量ε在[0,1]区间上服从均匀分布，则有概率密度函数

1,01

()0,x f x ≤≤?=??

其他

其分布函数为

0,0F(),01

x x x x x

=≤≤??≥?

1,1

在计算机上利用数学方法产生随机数的方法有平方取中法、移位

寄存器法和各种同余法等，这里我们采用剩同余法产生均匀随机序列，其递推公式为

1(mod )n n x x M λ+=

通常取2k M =，k 为计算机的尾数字长

021b x =+，b 为正整数，x 初始值须为奇数

323a λ=±，a 为正整数

(1) 所谓频率检验就是检验随机数序列的观测频数与理论频数的差异是否显著。详细的论述见教材P24，只要给定一个显著水平α，就可以确定二者之间的差异程度，α一般取值为0.05或0.01。为了有效的进行统计检验，N 值最好大于100，在雷达仿真中，通常1024N ≥。

(2) 随机数的参数检验是随机数分布的各个参数的观测值与理论值的差异是否显著，本设计的检验只包括一二阶矩和方差的检验。详细的论述见教材P25。

(3) 独立性检验就是检验随机数序列中前后各数的统计相关性是否异常，通常包括相关系数检验，联列表检验和连检验等，我们这里只进行相关系数的检验。详细的论述见教材P26，需要注意的是通常取50N j ->，这样就可以根据给定的显著水平α及正态分布表查处临界值了。

(4) 如此根据上述，本次设计的均匀分布的参数分布为：显著水平α取0.05，相应的频率检验中区间分为50k =组，取2χ分布的自由度为39(参考)时的临界值为254.572χ=；参数检验时相应的正态分布显著水平为0.05时的临界值为 1.96z α=；独立性检验时N 取值为100000点，取10050N j -=>，满足条件。

用Matlab 实现上述的100000点的均匀分布随机数的仿真并实现三种统计检验，结果如下图所示：

图1. 100000点均匀分布直方图

相应的检验结果为：

从上述检验结果可以看出，本设计产生的[0,1]区间上的均匀分布的随机数序列通过了各项检验。

由于均匀分布是以下9种其他分布的基础，即可以用均匀分布序列通过一定的数学变换得到，若均匀分布通过检验，则基于此产生的其他分布也必定能通过相关的检验，故下述的分布均是由通过检验的均匀分布变换得到，本身就不再进行检验了。

2．高斯分布

高斯分布即正态分布，其概率密度函数为

2

2

() ())

2

x

f x

μ

σ

-

=-

(1) 本次设计中的高斯分布是通过近似抽样法产生，即假设有N 个相互独立的随机变量12,,...,N u u u ，它们有相同的分布，其均值

()i E u m =，方差2()i D u σ=，根据中心极限定理，这N 个随机变量之和服从高斯分布，详细论述见教材P81。

(2) 这里我们取独立的随机变量i u 为均匀分布，通常其数目N 至少要大于8，这里我们将N 取为12，则高斯随机变量为

12

1

(6)j i i y u u σ==-+∑

为简便起见，其中u 取0，σ取1，仿真结果如下图所示

图2. 100000点高斯分布直方图

其中共检验了32组参数相同且相互独立的均匀分布，取其中符合要求的12组仿真生成高斯变量。 3．指数分布

指数分布的概率密度函数为

1,0()0,0x

e x

f x x λ

λ-?≥?=??

(1) 指数分布的仿真采用直接抽样法，详细论述见教材P74，即

知道某一分布的概率密度函数，就可以由其积分形式构成[0,1]区间上的均匀分布从而得到其随机变量。 (2) 由均匀分布得到的指数的随机变量的表达式为

1

ln i i y u λ

=-

其中i u 即为均匀分布随机变量，参数λ取8。仿真结果如下图所

示：

图3. 100000点高斯分布直方图

高斯分布的检验包括之后的7个检验与第一个均匀分布完全相同不再赘述。 4．广义指数分布

广义指数分布的概率密度函数为

0()exp[()]f x x s I =-+

式中s 是信噪比，

由教材P95-96的论述可以得到最终的广义指数分布的的随机变量抽样表达式为

12ln )i i i x u u s π=-++

其中s 是输入信噪比，1i u 、2i u 均为[0,1]区间上的均匀分布，仿真时取信噪比为2，方针的结果如下图所示

图4. 100000点广义指数分布直方图

5．瑞利分布

瑞利分布的概率密度函数为

2

22()exp(),02x

x f x x σσ

=-≥

式中，σ是瑞利分布的参量，但不是其均方根值。

利用直接抽样法得到瑞利分布的随机变量，具体的公式变换详见P77-78，最后得到瑞利分布直接抽样公式为

i ξσ=

其中σ

，i u 为[0,1]区间上的均匀分

布，仿真结果如下图所示

图5. 100000点瑞利分布直方图

6．广义瑞利分布

广义瑞利信号就是所谓的莱斯信号，它是将一个恒值信号叠加在两个相互独立正交的高斯随机变量值上，并取其矢量和而构成的。其概率密度函数为

220222()exp()()2r

r a ar

f x I σσσ

+=-

σ为随机变量r 的分布参数，a 为常数

(1) 因此，只要在正态分布总体中抽取两个相互独立的均值为0的正态分布随机数，再在其中一个加上个常数a ，便可获得广义瑞利分布的随机数。由P94论述得到其分布的表达式为

i r =式中i x 、i y 均为正态分布，且相互独立，a 为常数值 (2) 仿真时，取a 为2，结果如下图所示

图6. 100000点广义瑞利分布直方图

7．韦尔分布

服从三参量的韦尔分布的概率密度函数为

1()()exp(())0,0,0a a

n n a x x x x f x a b x b b b

---=->>≤<∞

式中，n x 、a 、b 是韦尔分布的位置参量、形状参量和标度参量。 利用直接抽样法得到韦尔分布的随机变量，具体的公式变换详见P77-78，最后得到韦尔分布直接抽样公式为

i n x b ξ=+

其中公式中位置参数n x 取3，形状参量a 取3，标度参量b 取2，

i u 为[0,1]区间上的均匀分布，仿真结果如下图所示

图7. 100000点韦尔分布直方图

8．拉普拉斯分布

拉普拉斯分布的概率密度函数为

|)|exp(2

)(m x a a

x f la --=

式中m 为均值，a 为形状参数。

为简便起见这里只考虑m 为0，a 为1的情况，即

|)|exp(2

1

)(x x f -=

可以看出，该分布为双指数分布，则得出有两个相同指数分布随机变量之差服从拉普拉斯分布的结论，最后有

)ln(21i

i

i u u =ξ

式中1,2i i u u 为 [0,1]区间的均匀分布随机数。

仿真时，直接调用产生两个均匀分布随机序列即可，结果如下图所示

图8. 100000点拉普拉斯分布直方图

9．柯西分布

柯西分布有概率密度函数

21

()(1)cau f x x π=

+

其中均值()0E x =，方差不存在。 随机数产生公式

a u

b i i +-=)]2

1

(tan[πξ

式中：i u 为[0,1]区间的均匀分布随机数。 仿真时参数a 、b 分别取,0、1，结果如下图所示

图9. 100000点柯西分布直方图

10．2

χ分布

2χ分布的概率密度函数为

212

21exp(),02

()2()

20,

0n n

x x x n f x p x χ+?->??=???≤? 根据正态分布之平方和为2

χ分布这一基本原理，可以直接获得

2χ分布随机变量抽样

21

()n

j j n x ξ==∑

式中j x 是采样间相互独立地均值为0，方差为1的正态分布随机变量：n 是2

χ分布的自由度，则可以看出自由度n 是决定该分布的唯一分布参量。

这里为方便起见，取自由度为4，则其为两个指数分布之和的分布，仿真结果为

图10. 自由度为4的2 分布直方图

五、附录

0.均匀分布

function u=evenlydis(a,b,m,N) %3?í?óà·¨

lam=8*a-3;

M=pow2(m);

x=zeros(1,N);

x(1)=pow2(b)+1;

for i=2:N;

y=lam*x(i-1);

x(i)=mod(y,M);

end

u=x/M;

end

1.分布检验

function result=check(x)

%?ì?é????°üà¨:?μ?ê?ì?é,ò??×??,?t?×??,·?2?,?àá￠D?; N=length(x);

%%%%%%%%%%%%%?μ?ê?ì?é,??μ?í3??á?A%%%%%%%%%%%%%

L=50; %?±·?í?é?μ?????′óD?

A=0;

n=zeros(1,L);

for i=1:N

for j=1:L

if

(x(i)>=(j-1)*1/L)&&(x(i)<=j*1/L) %?μ?èí3??

n(j)=n(j)+1;

end

end

end

for i=1:L

A=A+((n(i)-N/L)^2)/(N/L); %???±·?í?é?μ?????????è?μ?′ó·?2???ò?

end

%%%%%%%%%%%%%2?êy?ì?é,??μ?í3??á?Z1,Z2,Z%%%%%%%%%%%%% M1=0;

M2=0;

for i=1:N

M1=M1+x(i);

M2=M2+x(i)^2;

end

M1=M1/N;

M2=M2/N;

S=M2-M1+1/4;

Z1=sqrt(12*N)*(M1-1/2);%ò??×??í3??á?

Z2=1/2*sqrt(45*N)*(M2-1/3);%?t?×??í3??á?

Z=(sqrt(180*N))*(S-1/12);%·?2?í3??á?%%%%%%%%%%%%%?àá￠D??ì?é,??μ??à1??μêyí3??á?p%%%%%%%%%% %%%?

j=N-100; %í¨3￡è?N-j>50?′?é￡?′?′|è?100

sum=0;

for i=1:(N-j)

sum=sum+x(i)*x(j+i);

end

p=(sum/(N-j)-M1*M1)/S*sqrt(N-j); %%%%%%%%%%%%%?ù?Yò?é?í3??á?￡??ì?é???úDòáDx?ü·?í¨1y?ì?é%%%%%%%%%%%%%

ALPHA=0.05; %????????

Guass_Value=1.96; %±ê×??yì?·?2??????????a5%μ?áù???μLamenda_Value=54.572; %×?óé?è?a39μ?x^2·?2??????????a5 %μ?áù???μ

%%%%%??ê??á1?%%%%%%%

if

(abs(A)

disp('?ù?è·?2?í¨1y?ì?é');

result=1;

else

disp('?ù?è·?2??′í¨1y?ì?é'); %???°×￠êíμ?μ??-òòê?ê2?′

result=0;

end

disp('í¨1y????￡??á1?è???￡o');

disp(' ?μ?ê?ì?é ò??×?? ?t?×?? ·?2?

?à1??μêy ×?óé?è ????D????? áù???μ1 áù???μ2'); disp([A,Z1,Z2,Z,p,39,ALPHA,Guass_Value,Lamenda_Value] );

end

2.高斯分布

function z=guassdis(a,b,m,u,sigma,N)

%??oˉêyó?à′2úéú???1DòáD

%oˉêyμ÷ó?D?ê??az=GuassDist(s,u,sigma,N)

%s?a?ùó?·?2?μ?a

%b?a?ùó?·?2?μ?b

%m?a?ùó?·?2?μ?m

%u?a???1·?2?μ??ù?μ

%sigma?a???1·?2?μ?·?2?

undis=zeros(12,N);

i=1;

x0=b;

z=zeros(1,N);

seed=zeros(1,12);

%%%%%%%2úéú12???àá￠μ?í¨1y?ì?éμ?[0,1]?????ù?è·?2?μ????úDòáD￡?3¤?è?aN%%%%%%%%%

while(i<=12)

undis(i,:)=evenlydis(a,x0,m,N);

s=check(undis(i,:));

while(s==0)

x0=x0+1;

undis(i,:)=evenlydis(a,x0,m,N);

s=check(undis(i,:));

end

seed(i)=x0;

x0=x0+1;

i=i+1;

end

disp('?ùó?·?2?μ?bμ??μ?a￡o');

disp(seed);

for i=1:N

for j=1:12

z(i)=z(i)+undis(j,i);

end

z(i)=z(i)-6;

z(i)=z(i)*sigma+u;

end

end

3.指数分布

function z=expdis(a,b,lamda,m,N)

%??oˉêyó?à′2úéú2?êy?alamdaμ???êy·?2????úDòáD,3¤?è?aN? %s?a?ù?è·?2?μ?a?

%...

%lamda?a??êy·?2?μ?2?êy

z=zeros(1,N);

x0=b;

%%%%%%%%%%%%ê×?è2úéúí¨1y?ì?éμ?[0,1]????μ??ù?è·?2????úDòáD,3¤?è?aN%%%%%%%%%%%%%%?

undis=evenlydis(a,x0,m,N);

s=check(undis);

while(s==0)

x0=x0+1;

undis=evenly(a,x0,m,N);

s=check(undis);

end

for i=1:N

z(i)=-log(undis(i))/lamda;

end

end

4.广义指数分布

function z=eexpdis(a,b,sn,m,N)

%??oˉêyó?à′2úéú2?êy?asnμ?1?ò???êy·?2????úDòáD,3¤?è?aN ?

%sn?aê?è?D???±è?

%a?a?ùó?·?2?μ?3?ê?2?êy

%...

undis=zeros(2,N);

i=1;

x0=b;

z=zeros(1,N);

%%%%%%%2úéú2???àá￠μ?í¨1y?ì?éμ?[0,1]?????ù?è·?2?μ????úDòáD￡?3¤?è?aN%%%%%%%%%?

while(i<=2)

undis(i,:)=evenlydis(a,x0,m,N);

s=check(undis(i,:));

while(s==0)

x0=x0+1;

undis(i,:)=evenlydis(a,x0,m,N);

s=check(undis(i,:));

end

x0=x0+1;

i=i+1;

end

for i=1:N

z(i)=(-log(undis(1,i)))+2*sqrt(-sn*log(undis(1,i)))*c os(2*pi*undis(2,i))+sn;

end

end

5.瑞利分布

function z=rayleighdis(a,b,m,sigma,N)

%??oˉêyó?à′2úéú2?êy?asigmaμ?èeà?·?2????úDòáD,3¤?è?aN? %a?a?ù?è·?2?3?ê?2?êy

z=zeros(1,N);

x0=b;

%%%%%ê×?è2úéúí¨1y?ì?éμ?[0,1]????μ??ù?è·?2????úDòáD,3¤?è?aN%%%%%

undis=evenlydis(a,x0,m,N);

s=check(undis);

while(s==0)

x0=x0+1;

undis=evenlydis(a,x0,m,N);

s=check(undis);

end

for i=1:N

z(i)=sigma*sqrt((-2)*log(undis(i)));

end

end

6.广义瑞利分布

function z=erayleighdis(a,b,m,n,N)

%GYRelayDist?2úéú1?ò?èeà?·?2????úêyDòáD

undis=zeros(2,N);

i=1;

x0=b;

M1=zeros(1,N);

M2=zeros(1,N);

z=zeros(1,N);

%%%%%%%2úéú2???àá￠μ?í¨1y?ì?éμ?[0,1]?????ù?è·?2?μ????úDòáD￡?3¤?è?aN%%%%%%%

while(i<=2)

undis(i,:)=evenlydis(a,x0,m,N);

s=check(undis(i,:));

while(s==0)

x0=x0+1;

undis(i,:)=evenlydis(a,x0,m,N);

s=check(undis(i,:));

end

x0=x0+1;

i=i+1;

end

for i=1:N

M1(i)=sqrt(-2*log(undis(1,i))).*cos(2*pi*(undis(2,i)) );

M2(i)=sqrt(-2*log(undis(1,i))).*sin(2*pi*(undis(2,i)) );

z(i)=sqrt((M1(i)+n)^2+M2(i)^2);

end

7.韦尔分布

function z=vaildis(a,b,m,xn,beta,lamda,N)

%??oˉêyó?à′2úéú?¤??·?2????úDòáD,3¤?è?aN

%s?a?ù?è·?2?μ?a?

%...

%lamd?a?¤??·?2?μ?2?êy

z=zeros(1,N);

x0=b;

%%%%%%%%%%%%ê×?è2úéúí¨1y?ì?éμ?[0,1]????μ??ù?è·?2????ú

DòáD,3¤?è?aN%%%%%%%%%%%%%%?

undis=evenlydis(a,x0,m,N);

s=check(undis);

while(s==0)

x0=x0+1;

undist=evenly(a,x0,m,N);

s=check(undist);

end

for i=1:N

t=1/lamda;

y=-log(undis(i));

z(i)=xn+beta*power(y,t);

end

end

8.拉普拉斯分布

function z=laplacedis(a,b,m,N)

%??oˉêyó?à′2úéú3¤?è?aNμ?Laplacian·?2????úDòáD?

%b?a3?ê??ù?è·?2?μ?2?êy

undis=zeros(2,N);

i=1;

x0=b;

z=zeros(1,N);

%%%%%%%2úéú2???àá￠μ?í¨1y?ì?éμ?[0,1]?????ù?è·?2?μ????úDòáD￡?3¤?è?aN%%%%%%%%%?

while(i<=2)

undis(i,:)=evenlydis(a,x0,m,N);

s=check(undis(i,:));

while(s==0)

x0=x0+1;

undis(i,:)=evenlydis(a,x0,m,N);

s=check(undis(i,:));

end

x0=x0+1;

i=i+1;

end

for i=1:N

z(i)=log(undis(1,i)/undis(2,i));

end

end