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

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

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

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

【版权声明】本文为华为云社区用户原创内容,转载时必须标注文章的来源(华为云社区)、文章链接、文章作者等基本信息, 否则作者和本社区有权追究责任。如果您发现本社区中有涉嫌抄袭的内容,欢迎发送邮件进行举报,并提供相关证据,一经查实,本社区将立刻删除涉嫌侵权内容,举报邮箱: cloudbbs@huaweicloud.com
  • 点赞
  • 收藏
  • 关注作者

评论(0

0/1000
抱歉,系统识别当前为高风险访问,暂不支持该操作

全部回复

上滑加载中

设置昵称

在此一键设置昵称,即可参与社区互动!

*长度不超过10个汉字或20个英文字符,设置后3个月内不可修改。

*长度不超过10个汉字或20个英文字符,设置后3个月内不可修改。

举报
请填写举报理由
0/200