给定概率分布的随机变量仿真

给定概率分布的随机变量仿真

【摘要】在系统仿真的过程中经常需要产生不同分布的随机变量。利用混合同余法计算机可以产生周期很长的伪随机序列，以此作为均匀分布随机变量序列。在均匀分布随机变量的基础上，基于概率积分变换定理，利用计算机可以很方便地产生其他分布的随机变量，如指数分布、三角分布、正态分布随机变量。高斯白噪声通过线性系统产生ARMA信号序列。

在很多系统仿真的过程中，需要产生不同分布的随机变量。利用计算机可以很方便地产生不同分布的随机变量，各种分布的随机变量的基础是均匀分布的随机变量，有了均匀分布的随机变量，就可以用函数变换等方法得到其他分布的随机变量。

基于概率积分变换定理，如果X是分布函数为F(x)的随机变量，且分布函数F(x)为严格单调升函数，令Y=F(x)，则Y为[0,1]上均匀分布的随机变量。反之，若Y是在[0,1]上均匀分布的随机变量，那么X=F-1(Y)就是分布函数为F(x)的随机变量。这就是随机变量产生的变换法。实际上，当求分布函数的反函数变得复杂时，产生非均匀分布随机变量还有其他方法。

（1）均匀分布随机变量

一个均匀分布随机变量的样本序列由一个一个的随机数组成，由于计算机算术单元只能由有限字长的数字来表示，计算出来的样本数也是有限的。当计算机的算术单元是8位时，它只能表示256个不同的数，当计算机的算术单元是16位时，它只能表示65536个不同的数。与真正均匀分布的连续随机变量相比，这些样本并不是连续地占据某个取值区间。因此常称计算机算出来的随机数为伪随机数。

利用混合同余法产生均匀分布随机变量的递推同公式为：

x(n+1) = a*x(n) + c (MOD M)

y(n+1) = x(n+1)/M

式中，MOD M为取模运算，y(n)是[0,1]上均匀分布随机数。

初值x(0)、乘子a、增量c和模M都取非负的整数，当c=0时得到x(n+1)=a*x(n) (MOD M),此时为乘同余法产生伪随机序列，取M=2的k次方(k>2)，可得到最大的周期为T=2^k-2的伪随机序列，为了得到最大的周期，x(0)及a的选择需满足以下条件：

1)、 乘子a和模M=2k互素，属于同余式：

a= 3(MOD 8) 或者 a=-3(MOD 8)

给定的同余类即a=8*I+3这里I为任意选定正数。

2)、初值x(0)为任一奇数。

在混合同余法产生随机数时，取M=2的k次方可得到T=2^k的整周期序列，为了得到T=2^k 的序列，参数a、c、x(0)的选取应满足下面的条件：

1)、乘子a和模M互素，属于同余式a=1(MOD 4)给定的同余类，即a=4*I+1

这里I为任意选定的正数。

2)、增量c为奇数，初值x(0)为任意非负数。

例如在混合同余法中，我们在PC机上产生随机数，可取T=230、c=227623267，x(0)=1，a有30 种选择，分别对应30 种相互独立序列。

令y(i)=x(i)/T(T=230)，则y(i)均匀分布于 (0,1)的噪声。a的30种选择现列如下：

10753813, 38969437, 10763125, 228181237, 107452989

469773661, 4898493, 590648085, 10784189, 590634277

228206429, 469793853, 10841757, 832230813, 107464661

107408213, 469808301, 590585901, 469770781, 228178005

349045949, 832159853, 10845149, 107441485, 228244373

349062285, 107380685, 10767581, 469795933, 10842045

为了使产生的伪随机数有较大的周期，更接近真正的随机数，M越大越好，一般取M=2的k次方(如M=65536 = 2^16，k=16)，周期大小除了与M有关外，还与其他几个参数有关。

x(n+1) = (4a+1)x(n) + (2b+1) (MOD M), y(n+1) = x(n+1)/M

式中，a和b为任意正整数。

例如：x(n+1) = (2053*x(n) + 13849) MOD 65536, x(0)=1，2053=513*4+1，13849=6924*2+1

y(n+1)=x(n+1)/65536，服从(0,1)上的均匀分布。

对于均匀分布概率密度函数：

均匀分布随机数x(i)的平均值为(b+a)/2，标准差为(b-a)^2/12。

所以x(i)=a+(b-a)y(i), y(i)符合0和1之间的均匀分布。

（2）指数分布随机变量

指数分布随机变量的概率分布密度函数f(x)和概率分布函数F(x)分别为：

a>0，指数分布随机数的平均值为1/a，标准差为1/(a*a)。

指数分布随机变量的计算方法：

式中，u是在(0,1)上均匀分布的随机数。

（3）三角分布随机变量

三角分布随机变量概率分布密度函数f(x) 和概率分布函数F(x)分别为：

三角分布随机变量计算：

（4）正态分布随机变量

均值为0方差为1的正态概率分布密度函数f(x)为：

如果取两个均匀分布于0和1之间的随机数u1, u2，利用二元函数变换得到：

x1和x2是两个独立的标准正态分布N(0,1)随机变量。对于均值为m，标准差s为的正态分布随机变量y，可由y=s*x+m变换得到。

（5）ARMA序列产生

设ARMA(n, m)线性系统的传递函数为：

让正态分布白噪声通过该线性系统产生ARMA序列。设w(k)为正态分布随机序列，作为系统的输入，x(k)为系统的输出，则x(k)是ARMA(n,m)序列：

x(k) = b(0)w(k) + b(1)w(k-1) + …… + b(m)w(k-m)

- a(1)x(k-1) - a(2)x(k-2) - … -a(n)x(k-n)

k=0, 1, 2, ……, N-1

这里假设当k小于0时，w(k)和x(k)等于零。

当输入系统的信号为高斯白噪声时，可产生AR、MA、ARMA序列信号：

当m=0, n>0时，产生AR时间序列。

当m>0, n=0时，产生MA时间序列。

当m>0, n>0时，产生ARMA时间序列。

例如：用如下ARMA(4,4)模型产生ARMA随机信号序列：

a(0)=1, a(1)=0.4424, a(2)=-0.2030, a(3)=0.4164, a(4)=0.8853

b(0)=1, b(1)=0.1792, b(2)=0.8202, b(3)=0.2676, b(4)=0.3300

该线性系统的极点为：

0.61977 + 0.74622i, 0.61977 - 0.74622i

-0.84097 + 0.48335i, -0.84097 - 0.48335i

用模和相位角表示为：

0.97exp(±j0.8777) = 0.97exp(±j0.1397*2π) = 0.97exp(±j13.97*2π/100)，对应频率f=14Hz

0.97exp(±j2.6199) = 0.97exp(±j0.4170*2π) = 0.97exp(±j41.7*2π/100)，对应频率f=41.7Hz

该系统的零点为：

0.25826 + 0.81034i, 0.25826 – 0.81034i

-0.34786 + 0.57897i, -0.34786 – 0.57897i

用模和相位角表示为：

0.8505exp(±j1.2623) = 0.8505exp(±j0.2009*2π) = 0.8505exp(±j20.09*2π/100)，对应频率f=20.1Hz

0.6754exp(±j2.1118) = 0.6754exp(±j0.3361*2π) = 0.6754exp(±j33.61*2π/100)，对应频率f=33.6Hz

以下是一段ARMA(4,4)的序列，采样率为100Hz，长度为1分钟(显示6行每行10秒)。

对应的信号功率谱为：

信号的幅度分布直方图为：

参考文献

[1] 赵淑清，郑薇 编著. 随机信号分析[M]. 哈尔滨工业大学出版社，1999年6月第1版，第38-43页。
[2] 徐士良 编. C常用算法程序集[M]. 北京 清华大学出版社，1994年1月第1版，第646-650页。
[3] 马振华 主编. 现代应用数学手册—概率统计与随机过程卷[M]. 北京 清华大学出版社，2000年7月第1版，第723-732页。
[4] 张阿舟 主编. 振动数字信号处理程序库[M]. 北京 科学出版社，1988年8月1版。
[5] 北方交通大学信息所 编，近代数字信号处理通用程序[M]. 北京 科学出版社，1988年9月第1版。

联系作者：chengbowork@163.com

相关文章

多道信号分析软件系统
生物医学信号处理与分析软件系统设计
曲线拟合软件
离散小波变换用于信号滤波
多维特征参数机器学习算法
多维特征参数机器学习软件
Kohonen自组织特征映射神经网络(环形和球面形网络)
矩阵的三维图形显示软件
图片浏览软件工具
主成分分析(K-L变换)与信号的分解与合成(滤波)
信号的样本熵序列计算
信号的双谱分析
信号的经验模态分解(EMD)
希尔伯特(Hilbert)变换信号瞬时频率计算