泊松过程的生成及其统计分析.docx
《泊松过程的生成及其统计分析.docx》由会员分享，可在线阅读，更多相关《泊松过程的生成及其统计分析.docx（15页珍藏版）》请在冰豆网上搜索。
泊松过程的生成及其统计分析
泊松过程的生成及其统计分析
实验报告
班级：
6041
姓名：
韩丽媛
学号：
3116036015
一、实验题目
假设一个交换系统有M部电话，每个用户在很短的时间
（单位时间内）呼叫一次的概率为P；用户间呼入的时刻相互独立，当M很大，P很小时，时间t内到达交换机的呼叫次数构成泊松过程N（t）。
1、确定此泊松过程的参数
。
2、利用计算机仿真N（t）的生成过程。
注意合理选择M和P，时间分辨率为一个单位时间
。
3、为了比较生成的N（t）与理论模型的吻合程度。
取N（t）的多个样本并选取3个典型时间
，得到
三个随机变量的样本，在一张图上画出其直方图及理论分布曲线，并将两者对照。
比较M选取不同时的效果。
注意：
样本个数足够多。
4、验证N（t）的增量平稳性。
5、画出任意相邻两次呼叫间隔的直方图，和理论值进行对照。
验证其与其它相邻两次呼叫间隔随机变量的独立性。
二、实验过程
1、确定此泊松过程的参数
由题目容易知道，在很短的时间
内M个用户的呼叫一次的概率为MP，而由定义知道，
时间内到达交换机的呼叫一次的概率为
，故有
（1）
从而有
。
2、利用计算机仿真N（t）的生成过程
对每个用户，在
时间内呼叫一次的概率P很小，可以用rand函数生成一组[0,1]的随机数，当随机数小于P时，则认为有呼叫，将其置为1，否则认为没有呼叫，置为0；有M部电话，则生成M组[0,1]的随机数，对每组随机数用上诉方法得到一个只有0和1的逻辑矩阵，用来表示某一时刻是否有呼叫。
下面是
M=3000,总时间为T=5的实验结果：
图1N（t）的生成结果
可以看到呼叫的计数过程，是递增的，并且可以计算，时间T=5内呼叫总次数平均为
多次时间结果最后的呼叫次数都在15次左右。
程序：
clc
clear
closeall
p=10^（-6）;
M=3000;
dt=0.001;
T=5;
x=rand（M,T/dt）;
y=[];
fori=1:
M
forj=1:
T/dt
ifx（i,j）x（i,j）=1;
else
x（i,j）=0;
end
end
end
y=（sum（x）~=0）;
m=[];
m
（1）=0;
fori=1:
T/dt
m（i+1）=m（i）+y（i）;
end
t=1:
T/dt+1;
t=t*dt;
plot（t,m）
此外，matlab中还有二项分布生成函数binornd，可以用x=binornd（1,p,M,T/dt）代替中间的两个for循环，这个函数的功能是对一个发生概率为P的事件随机试验一次，若发生置为1，不发生置为0，此实验要对M个电话实验T/dt次，故生成的是M行，T/dt的矩阵，运行结果是一样的。
3、比较生成的N（t）与理论模型的吻合程度
（1）
的统计直方图和理论分布曲线
下面是
M=3000,总时间为T=1.2，选取时间t1=0.3，t2=0.6，t3=0.9作2000次试验统计的实验结果：
图2
的统计直方图和理论分布曲线
在图2中，圆圈代表
的统计直方图，正方形代表
的统计直方图，五角星代表
的直方图。
从图中可以看出，虽然有较小的误差，但是生成的N（t）和理论模型还是基本吻合的。
程序中主要用到了直方图统计函数hist，生成max（Nt1）-min（Nt1）个直方条间的间隔刚好是1，此时的坐标分别为0.5、1.5、2.5……，并且0.5的直方条包括了0次呼叫和1次呼叫的的概率，1.5、2.5、3.5等等依次代表的是2次、3次、4次呼叫的概率，因而有了程序中的相关修正。
程序：
clc
clear
closeall
p=5*10^（-6）;
M=3000;
dt=0.003;
a=M*p/dt;
T=1.2;
loop=2000;
t1=0.3;
t2=0.6;
t3=0.9;
fork=1:
loop%作loop次试验
x=rand（M,T/dt）;
fori=1:
M
forj=1:
T/dt
ifx（i,j）x（i,j）=1;
else
x（i,j）=0;
end
end
end
tt=dt*find（（sum（x）~=0）==1）;%每次试验各个呼叫发生的时刻
Nt1（k）=sum（ttNt2（k）=sum（ttNt3（k）=sum（ttend
[N1,index1]=hist（Nt1,max（Nt1）-min（Nt1））;%（0，t1）内呼叫次数的统计直方图
[N2,index2]=hist（Nt2,max（Nt2）-min（Nt2））;
[N3,index3]=hist（Nt3,max（Nt3）-min（Nt3））;
index1=[min（Nt1）,index1+0.5];%作相关修正
index2=[min（Nt2）,index2+0.5];
index3=[min（Nt3）,index3+0.5];
N1=[sum（Nt1==min（Nt1））,N1
（1）-sum（Nt1==min（Nt1））,N1（2:
end）];
N2=[sum（Nt2==min（Nt2））,N2
（1）-sum（Nt2==min（Nt2））,N2（2:
end）];
N3=[sum（Nt3==min（Nt3））,N3
（1）-sum（Nt3==min（Nt3））,N3（2:
end）];
p1=[];p2=[];p3=[];
fork=1:
length（index1）
p1=[p1,（a*t1）^index1（k）*exp（-a*t1）/factorial（index1（k））];%理论值
end
fork=1:
length（index2）
p2=[p2,（a*t2）^index2（k）*exp（-a*t2）/factorial（index2（k））];
end
fork=1:
length（index3）
p3=[p3,（a*t3）^index3（k）*exp（-a*t3）/factorial（index3（k））];
end
stem（index1,N1/loop,'r'）;
holdon
plot（index1,p1,'r'）
holdon
stem（index2,N2/loop,'bs'）;
holdon
plot（index2,p2,'b'）
holdon
stem（index3,N3/loop,'gp'）;
holdon
plot（index3,p3,'g'）
holdon
（2）比较M不同时的实验效果
对于上面的参数，我们选择t2时刻，M分别取1000、2000、3000得到
的统计直方图如图3所示，圆形对应的是M=1000，正方形对应的是M=2000，五角星对应的是M=3000，从图3中可以看到，当M值增大时，直方图和；理论曲线都往右移动，从理论上分析，在P和
不变时，M值越大，强度常数
越大，相同时间内呼叫的次数更多，所以在呼叫次数多的地方概率更大，曲线往右移动。
图3M不同时的实验效果对比
4、验证N（t）的增量平稳性
增量平稳性数学表示为，对任何s和t，P{N（s+t）-N（s）=n}=P{N（t）=n}，即在相同时间内呼叫n次的概率相等。
下图是取了三个相等的时间间隔
进行的呼叫次数的直方图统计结果：
图4增量平稳性验证曲线
由于只需要相同时间内呼叫相同次数的概率相同，为了简化程序和计算量，在直方图统计中没有对第一个直方条进行修正，并不影响实验的结论，从图4中可以看到，三个相等的时间间隔呼叫次数的概率分布曲线基本重合，说明相同时间间隔内呼叫次数相同的概率基本相同，从而验证了增量平稳性。
程序：
clc
clear
closeall
p=5*10^（-6）;
M=3000;
dt=0.003;
T=1.8;
loop=2000;
fork=1:
loop
x=rand（M,T/dt）;
fori=1:
M
forj=1:
T/dt
ifx（i,j）x（i,j）=1;
else
x（i,j）=0;
end
end
end
y=（sum（x）~=0）;
m=[];
m
（1）=0;
fori=1:
T/dt
m（i+1）=m（i）+y（i）;
end
Nt1（k）=m（201）-m
（1）;
Nt2（k）=m（401）-m（201）;
Nt3（k）=m（601）-m（401）;
end
[N1,index1]=hist（Nt1,max（Nt1）-min（Nt1））;
[N2,index2]=hist（Nt2,max（Nt2）-min（Nt2））;
[N3,index3]=hist（Nt3,max（Nt3）-min（Nt3））;
plot（index1,N1/loop,'r'）;holdon;
plot（index2,N2/loop,'b'）;holdon;
plot（index3,N3/loop,'g'）;holdon;
5、
（1）画出任意相邻两次呼叫间隔的直方图，和理论值进行对照。
由理论可知，任意两次的呼叫间隔的概率分布函数为负指数分布：
，下面是
M=3000,总时间为T=3，选取第二次和第一次呼叫的时间间隔得到的统计实验结果：
图5呼叫时间间隔分布直方图
从图5中可以看出，相邻两次呼叫间隔满足负指数分布，与理论相符。
编程时，将时间间隔平均分在50个直方条中，在求理论值时，需要对负指数型概率密度函数在每个直方条中求积分，需要注意的是积分的区间。
程序：
clc
clear
closeall
p=5*10^（-6）;
M=3000;
dt=0.003;
a=M*p/dt;
T=3;
loop=3000;
fork=1:
loop
x=rand（M,T/dt）;
fori=1:
M
forj=1:
T/dt
ifx（i,j）x（i,j）=1;
else
x（i,j）=0;
end
end
end
y=sum（x）~=0;
tt=dt*find（y==1）;
%fori=1:
length（tt）-1
%b（i）=tt（i+1）-tt（i）;
%end
b（k）=tt
（2）-tt
（1）;
c（k）=tt（4）-tt（3）;
end
[N1,index1]=hist（b,50）;
dh=（max（b）-min（b））/100;
stem（index1,N1/loop,'r'）;
holdon
fori=1:
50
t=（index1（i）-dh）:
0.001:
（index1（i）+dh）;
l=a*exp（-a*t）;
q（i）=trapz（t,l）;
end
plot（index1,q）
Eb=sum（b）/length（b）;
Ec=sum（c）/length（c）;
Ebc=sum（b.*c）/length（b）;
Db=sum（b.*b）/length（b）-Eb^2;
Dc=sum（c.*c）/length（c）-Ec^2;
Covbc=Ebc-Eb*Ec;
rov=Covbc/sqrt（Db）/sqrt（Dc）
（2）验证其与其它相邻两次呼叫间隔随机变量的独立性
可以再选取第四次和第三次呼叫的时间间隔来验证，要验证其独立性，我们可以求它们的相关系数，程序见上面，计算可以得到相关系数rov=0.0107，M=4000时，rov=0.0106，,当M增大时，rov减小，当M足够大时，rov将接近0，独立性得证，但是这是不严格的，如果要整个证明两个随机变量x、y的独立性，我们需要证明
，从数值上来证明比较困难。
（范文素材和资料部分来自网络，供参考。
可复制、编制，期待你的好评与关注）
升级会员 