通过编程实现poisson过程的模拟
应用随机过程实验报告实验序号：1-4日期：2013年5月30日姓名梁光佐学号201005050110实验题目应用随机过程综合实验实验所用软件及版本MATLAB20081、实验目的(1) 通过编程实现poisson过程的模拟，运用matlab画图这样更直观的了解poisson过程，(2) 运用计算机通过编程来辅助解题，这样解决了解题的繁琐，使解题的效力提高了，也节约了时间。2、实验内容实验一实验问题1.编制程序产生并输出100个二项分布的随机数，n10'p0.6.2.进行三次Poisson过程的模拟，“3,n=50,n=100,n=200作图：(在同一直角坐标系下，作出''"()'的关系图实验二一、泊松过程的模拟1. 基本原理根据服务系统接受服务顾客数服从泊松分布这一模型可知，X(n),°是一个计数过程，%几，n三丄'是对应的时间间隔序列，若几)(n=1,2,.)是独立同分布的均值为的指数分布，贝X(n),°是具有参数为久的泊松。2. 具休实现过程实现步骤如下：(1) T.由函数random('exponential',lamda)构造服从指数分布的一：序列。(2) .根据服务系统模型，一一=-o.对任意(：,：-),X(t)=n,由此得到泊松过程的模拟。(1) 3.过程模拟验证设定t=0时刻，计数为0,满足X(0)=0这一条件。(2) TT-：是由random('exponential',lamda)生成，一：间相互独立。二、泊松过程的检验1.检验方法Kolmogorov-Smirnov检验(柯尔莫哥洛夫-斯摩洛夫)，亦称拟合优度检验法，用来检用来检验模拟所得的数据的分布是不是符合一个理论的已知分布。检验步骤及过程：(1) 条件设定：H1:实验产生模拟泊松分布数据的总体分布服从泊松分布。H0:实验产生模拟泊松分布数据的总体分布不服从泊松分布。(2) 检验准备：对于H1,已经假定所产生模拟泊松过程数据X(“)服从泊松分布，而强度九未知，利用函数poissfit(x,alpha)估算出模拟泊松过程的强度几，再利用函数poisscdf(x,lamda)得到泊松分布的累积分布函数P。Kolmogorov-Smirnov检验直接调用Kolmogorov-Smirnov检验函数kstest(x,x,p,alpha)，其中，x为输入模拟泊松序列，P为累积分布函数，1-alpha为置信区间，当结果H1=1时，则输入数据位泊松分布，否则，不是泊松分布。实验三学习以下函数1.Definition2.poissrnd3.poisspdf4.poissinv5.poissfit6.poisscdf实验四：(1) 、PPT例2维正态密度与二维正态密度、赌博问题模拟(带吸收壁的随机游走：结束N次游走所花的平均时间及终止状态分布规律)3、详细设计(包括算法描述和程序)实验一程序代码%该函数用于输出服从所给定的Xb(k,p1)的m个随机数%e,E分别为所输出随机数的和原先的均值%v，V分别为所输出随机数的和原先的方差%p值为给定的Xb(k,p1)的分布列中的概率functionr=bino(m,n,k,p1)g=0:n;y=binocdf(g,k,p1)fori=2:(k+1)p(1)=y(1);p(i)=y(i)-y(i-1);endR=rand(1,m)forj=1:mforl=2:kif0<R(j)&R(j)<p(1)x(j)=0;elseify(l-1)<R(j)&R(j)<y(l)x(j)=l-1;endendendendPxe=mean(x)E=k*p1v=var(x)V=k*p1*(1-p1)end实验二程序代码%该程序用于输出Poisson过程的模拟%m1,m2,m3为Poisson过程的最终人数，n为Poisson过程的强度%x,y,z分别ml，m2，m3个人所到达的的时间间隔%t1,t2,t3分别ml，m2，m3个人所到达的的时刻%k1,k2,k3检验是否与Poisson过程的强度相近functiong=poisson1(m1,m2,m3,n)R1=rand(1,m1);R2=rand(1,m2);R3=rand(1,m3);fori=1:m1x(i)=-(log(1-R1(i)/n;endxfori=1:m2y(i)=-(log(1-R2(i)/n;endyfori=1:m3z(i)=-(log(1-R3(i)/n;endzfori=2:m111(1)=x(1);11(i)=t1(i-1)+x(i);endfori=2:m2t2(1)=y(1);t2(i)=t2(i-1)+y(i);endfori=2:m3t3(1)=z(1);13(i)=t3(i-1)+z(i);endt1t2t3k1=m1/t1(m1)k2=m2/t2(m2)k3=m3/t3(m3)N1=1:m1;N2=1:m2;N3=1:m3;plot(t1,Nl,'g.')holdonplot(t2,N2,'r.')holdonplot(t3,N3,'.')gridonxlabel('t')ylabel('N')title('Poisson')实验三程序代码：1. Definition:2. x=0:15;y=poisspdf(x,5);plot(x,y,'+')x=(0:0.1:10)'p1=ncx2pdf(x,4,2);p=chi2pdf(x,4);plot(x,p,'-',x,p1,'-')x=(0.01:0.1:10.01)'p1=ncfpdf(x,5,20,10);p=fpdf(x,5,20);plot(x,p,'-',x,p1,'-')poissrnd：3. lambda=2;random_sample1=poissrnd(lambda,1,10)random_sample2=poissrnd(lambda,110)random_sample3=poissrnd(lambda(ones(1,10)poisspdf：4. p=poisspdf(0,2)poissinv：5. poissinv(0.95,2)mediandefects=poissinv(0.50,2)poissfit：6. r=poissrnd(5,10,2);l,lci=poissfit(r)poisscdf：probability=l-poisscdf(4,2)probability=poisscdf(4,4)实验四程序代码：()、%ppt例2维正态密度与二维正态密度symsxy;s=3;t=4;mu1=0;mu2=0;sigma1=sqrt(1+sA2);sigma2=sqrt(1+tA2);x=-6:0.1:6;f1=1/sqrt(2*pi*sigma1)*exp(-(x-mu1).A2/(2*sigma1A2);f2=1/sqrt(2*pi*sigma2)*exp(-(x-mu2).A2/(2*sigma2A2);plot(x,f1,'r-',x,f2,'k一.')symsxy;s=3;t=4;mu1=0;mu2=0;sigma1=sqrt(1+sA2);sigma2=sqrt(1+tA2);rho=(1+s*t)/(sigma1*sigma2);f=1/(2*pi*sigma1*sigma2*sqrt(1-rhoA2)*exp(-1/(2*(1-rhoA2)*(x-mu1)A2/sigma1A2-2*rho*(x-mu1)*(y-mu2)/(sigma1*sigma2)+(y-mu2)A2/sigma2A2);ezsurf(f)(2)A赌博问题模拟(带吸收壁的随机游走：结束n次游走所花的平均时间及终止状态分布规律)程序代码：p=1/2;m1=0;m2=0;N=1000;t1=0;t2=0;forn=1:1:Nm=0;a=5;whilea>0&a<10m=m+1;r=2*binornd(1,p)-1;ifr=-1a=a-1;elsea=a+1;endendifa=0t1=t1+m;m1=m1+1;elseB. t2=t2+m;m2=m2+1;endendfprintf('Theaveragetimesofarriving0and10respectivelyare%d,%dn',t1/m1,t2/m2);fprintf('Thefrequenciesofarriving0and10respectivelyare%d,%d.n',m1/N,m2/N).赌博问题模拟(带吸收壁的随机游走：结束N次游走所花的平均时间及终止状态分布规律)程序代码:C. randn('state',200)%setthestateofrandnT=1;N=1000;dt=T/N;dW=sqrt(dt)*randn(1,N);%incrementsW=cumsum(dW);%cumulativesumplot(0:dt:T,0,W,'r-')%plotWagainsttxlabel('t','FontSize',16)ylabel('W(t)','FontSize',16,'Rotation',0)赌博问题模拟(带吸收壁的随机游走：结束N次游走所花的平均时间及终止状态分布规律)程序代码：randn('state',100)%setthestateofrandnT=1;N=1000;dt=T/N;t=dt:dt:1;M=2000;%MpathssimultaneouslydW=sqrt(dt)*randn(M,N);%incrementsW=cumsum(dW,2);%cumulativesumU=exp(repmat(t,M1)+0.5*W);Umean=mean(U);plot(0,t,1,Umean,'b-'),holdon%plotmeanoverMpathsplot(0,t,ones(5,1),U(1:5,:),'r'),holdoff%plot5individualpathsxlabel('t','FontSize',16)ylabel('U(t)','FontSize',16,'Rotation',0,'HorizontalAlignment','right')legend('meanof1000paths','5individualpaths',2)aveerr=norm(Umean一exp(9*t/8),'inf')%sampleerror4、实验结果与分析实验一运行结果：clearbino(100,10,10,0.6)Columns1through11Columns12through220.13820.13930.08370.10460.0006z=Columns1through11Columns12through220.55650.90751.77890.64900.2341t1=Columns1through11Columns12through22实验三学习结果:1.Definitionx=0:15;y=poisspdf(x,5);plot(x,y,'+')p=chi2pdf(x,plot(x,p,'-',x1015x=(0.01:0.1:10.01)'pl=ncfpdf(x,5,20,10);p=fpdf(x,5,20);plot(x,p,'-',x,p1,'-')2.poissrnd：程序代码lambd
