泊松過程的生成及其統(tǒng)計分析.docVIP

泊松過程的生成及其統(tǒng)計分析實驗報告 班級：6041姓名：韓麗媛 學號：3116036015 一、實驗題目假設一個交換系統(tǒng)有M部電話，每個用戶在很短的時間（單位時間內(nèi)）呼叫一次的概率為P；用戶間呼入的時刻相互獨立，當M很大，P很小時，時間t內(nèi)到達交換機的呼叫次數(shù)構成泊松過程N(t)。1、 確定此泊松過程的參數(shù)。2、 利用計算機仿真N(t)的生成過程。注意合理選擇M和P，時間分辨率為一個單位時間。3、 為了比較生成的N(t)與理論模型的吻合程度。取N(t)的多個樣本并選取3個典型時間,，得到,三個隨機變量的樣本，在一張圖上畫出其直方圖及理論分布曲線，并將兩者對照。比較M選取不同時的效果。注意：樣本個數(shù)足夠多。4、 驗證N(t)的增量平穩(wěn)性。5、 畫出任意相鄰兩次呼叫間隔的直方圖，和理論值進行對照。驗證其與其它相鄰兩次呼叫間隔隨機變量的獨立性。二、實驗過程1、確定此泊松過程的參數(shù)由題目容易知道，在很短的時間內(nèi)M個用戶的呼叫一次的概率為MP，而由定義知道，時間內(nèi)到達交換機的呼叫一次的概率為，故有 （1）從而有。2、利用計算機仿真N(t)的生成過程對每個用戶，在時間內(nèi)呼叫一次的概率P很小，可以用rand函數(shù)生成一組0,1的隨機數(shù)，當隨機數(shù)小于P時，則認為有呼叫，將其置為1，否則認為沒有呼叫，置為0；有M部電話，則生成M組0,1的隨機數(shù)，對每組隨機數(shù)用上訴方法得到一個只有0和1的邏輯矩陣，用來表示某一時刻是否有呼叫。下面是,M=3000,總時間為T=5的實驗結果：圖1 N(t)的生成結果可以看到呼叫的計數(shù)過程，是遞增的，并且可以計算，時間T=5內(nèi)呼叫總次數(shù)平均為,多次時間結果最后的呼叫次數(shù)都在15次左右。程序：clcclearclose allp=10(-6);M=3000;dt=0.001;T=5;x=rand(M,T/dt);y=;for i=1:M for j=1:T/dt if x(i,j)p x(i,j)=1; else x(i,j)=0; end endendy=(sum(x)=0);m=;m(1)=0;for i=1:T/dt m(i+1)=m(i)+y(i);endt=1:T/dt+1;t=t*dt;plot(t,m)此外，matlab中還有二項分布生成函數(shù)binornd，可以用x=binornd(1,p,M,T/dt)代替中間的兩個for循環(huán)，這個函數(shù)的功能是對一個發(fā)生概率為P的事件隨機試驗一次，若發(fā)生置為1，不發(fā)生置為0，此實驗要對M個電話實驗T/dt次，故生成的是M行，T/dt的矩陣，運行結果是一樣的。3、比較生成的N(t)與理論模型的吻合程度（1）,的統(tǒng)計直方圖和理論分布曲線下面是,M=3000,總時間為T=1.2，選取時間t1=0.3，t2=0.6，t3=0.9作2000次試驗統(tǒng)計的實驗結果：圖2 ,的統(tǒng)計直方圖和理論分布曲線在圖2中，圓圈代表的統(tǒng)計直方圖，正方形代表的統(tǒng)計直方圖，五角星代表的直方圖。從圖中可以看出，雖然有較小的誤差，但是生成的N(t)和理論模型還是基本吻合的。程序中主要用到了直方圖統(tǒng)計函數(shù)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次呼叫的概率，因而有了程序中的相關修正。程序：clcclearclose allp=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;for k=1:loop %作loop次試驗 x=rand(M,T/dt); for i=1:M for j=1:T/dt if x(i,j)p x(i,j)=1; else x(i,j)=0; end end end tt=dt*find(sum(x)=0)=1); %每次試驗各個呼叫發(fā)生的時刻 Nt1(k)=sum(ttt1); %每次試驗在時間（0，t1）內(nèi)呼叫的次數(shù) Nt2(k)=sum(ttt2); Nt3(k)=sum(ttt3);endN1,index1=hist(Nt1,max(Nt1)-min(Nt1); %（0，t1）內(nèi)呼叫次數(shù)的統(tǒng)計直方圖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=;for k=1:length(index1) p1=p1,(a*t1)index1(k)*exp(-a*t1)/factorial(index1(k); %理論值endfor k=1:length(index2) p2=p2,(a*t2)index2(k)*exp(-a*t2)/factorial(index2(k);endfor k=1:length(index3) p3=p3,(a*t3)index3(k)*exp(-a*t3)/factorial(index3(k);endstem(index1,N1/loop,r);hold onplot(index1,p1,r)hold on stem(index2,N2/loop,bs);hold onplot(index2,p2,b)hold on stem(index3,N3/loop,gp);hold onplot(index3,p3,g)hold on（2）比較M不同時的實驗效果 對于上面的參數(shù)，我們選擇t2時刻，M分別取1000、2000、3000得到的統(tǒng)計直方圖如圖3所示，圓形對應的是M=1000，正方形對應的是M=2000，五角星對應的是M=3000，從圖3中可以看到，當M值增大時，直方圖和；理論曲線都往右移動，從理論上分析，在P和不變時，M值越大，強度常數(shù)越大，相同時間內(nèi)呼叫的次數(shù)更多，所以在呼叫次數(shù)多的地方概率更大，曲線往右移動。圖3 M不同時的實驗效果對比4、驗證N(t)的增量平穩(wěn)性 增量平穩(wěn)性數(shù)學表示為，對任何s和t，PN(s+t)-N(s)=n=PN(t)=n，即在相同時間內(nèi)呼叫n次的概率相等。下圖是取了三個相等的時間間隔進行的呼叫次數(shù)的直方圖統(tǒng)計結果：圖4 增量平穩(wěn)性驗證曲線由于只需要相同時間內(nèi)呼叫相同次數(shù)的概率相同，為了簡化程序和計算量，在直方圖統(tǒng)計中沒有對第一個直方條進行修正，并不影響實驗的結論，從圖4中可以看到，三個相等的時間間隔呼叫次數(shù)的概率分布曲線基本重合，說明相同時間間隔內(nèi)呼叫次數(shù)相同的概率基本相同，從而驗證了增量平穩(wěn)性。程序：clcclearclose allp=5*10(-6);M=3000;dt=0.003;T=1.8;loop=2000;for k=1:loop x=rand(M,T/dt); for i=1:M for j=1:T/dt if x(i,j)p x(i,j)=1; else x(i,j)=0; end end end y=(sum(x)=0); m=; m(1)=0; for i=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);endN1,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);hold on;plot(index2,N2/loop,b);hold on;plot(index3,N3/loop,g);hold on;5、（1）畫出任意相鄰兩次呼叫間隔的直方圖，和理論值進行對照。由理論可知，任意兩次的呼叫間隔的概率分布函數(shù)為負指數(shù)分布：， 下面是,M=3000,總時間為T=3，選取第二次和第一次呼叫的時間間隔得到的統(tǒng)計實驗結果：圖5 呼叫時間間隔分布直方圖從圖5中可以看出，相鄰兩次呼叫間隔滿足負指數(shù)分布，與理論相符。編程時，將時間間隔平均分在50個直方條中，在求理論值時，需要對負指數(shù)型概率密度函數(shù)在每個直方條中求積分，需要注意的是積分的區(qū)間。程序：clcclearclose allp=5*10(-6);M=3000;dt=0.003;a=M*p/dt;T=3;loop=3000; for k=1:loop x=rand(M,T/dt); for i=1:M for j=1:T/dt if x(i,j)p x(i,j)=1; else x(i,j)=0; end end end y=sum(x)=0; tt=dt*find(y=1);% for i=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);hold onfor i=1:50 t=(index1(i)-dh):0.001:(index1(i)+dh); l=a*exp(-a*t); q(i)=trapz(t,l);endplot(index1,q)Eb=sum(b)/length(b);Ec=sum(c)/length(c);Ebc=sum(b.*c)/length(b);Db=sum(b.*b)/length(b)-Eb2;Dc=sum(c.*c)/length(c)-Ec2;Covbc=Ebc-Eb*Ec;