2015-2016第一学期随机过程第二次上机实验报告实验目的：通过随机过程上机实验，熟悉Monte Carlo计算机随机模拟方法，熟悉Matlab的运行环境，了解随机模拟的原理，熟悉随机过程的编码规律即各种随机过程的实现方法，加深对随机过程的理解。

上机内容：(1 )模拟随机游走。

(2)模拟Brown运动的样本轨道。

(3)模拟Markov过程。

实验步骤：(1)给出随机游走的样本轨道模拟结果，并附带模拟程序。

①一维情形%—维简单随机游走% “从0开始，向前跳一步的概率为p，向后跳一步的概率为1-p”n=50;p=0.5;y=[0 cumsum(2.*(rand(1,n-1)v=p)-1)]; % n 步。

plot([0:n-1],y); %画出折线图如下。

w%一维随机步长的随机游动%选取任一零均值的分布为步长，比如，均匀分布。

n=50;x=rand(1,n)-1/2;y=[0 (cumsum(x)-l)];plot([0:n],y);②二维情形%在(u, v)坐标平面上画出点(u(k), v(k)), k=1:n,其中(u(k)) 和(v(k))是一维随机游动。

例%子程序是用四种不同颜色画了同一随机游动的四条轨道。

n=100000;colorstr=['b' 'r' 'g' 'y'];for k=1:4z=2.*(rand(2,n)<0.5)-1;x=[zeros(1,2); cumsum(z')];col=colorstr(k);plot(x(:,1),x(:,2),col);③%三维随机游走 ranwalk3dp=0.5;n=10000; colorstr=['b' 'r' 'g' 'y'];for k=1:4z=2.*(rand(3,n)v=p)-1; x=[zeros(1,3); cumsum(z')];col=colorstr(k);plot3(x(:,1),x(:,2),x(:,3),col);hold on end gridhold onendgrid4：04003?0-200-300-400-2OD20050、-100-200 -20D⑵给出一维，二维Brown运动和Poisson过程的模拟结果，并附带模拟程序，没有结果的也要把程序记录下来。

①一维Brown%这是连续情形的对称随机游动，每个增量W(s+t)-W(s)是高斯分布N(0, t)，不相交区间上的增量是独立的。

典型的模拟它方法是用离散时间的随机游动来逼近。

n=1000;dt=1;y=[0 cumsum(dt A0.5.*randn(1,n))]; % 标准布朗运动。

plot(0:n,y);②二维Brownnp oints = 5000;dt = 1;bm = cumsum([zeros(1,3); dt A0.5*randn(npoints-1,3)]); figure(1);plot(bm(:, 1), bm(:, 2), 'k');pcol = (bm-repmat(min(bm), npoints, 1))./ ...repmat(max(bm)-min(bm), npoints, 1);hold on;scatter(bm(:, 1), bm(:, 2),...10, pcol, 'filled');grid on;hold off;200③三维Brown叩oints = 5000;dt = 1;bm = cumsum([zeros(1,3); dt A0.5*randn(npoints-1,3)]); figure(1); plot3(bm(:, 1), bm(:, 2), bm(:, 3), 'k');pcol = (bm-repmat(min(bm), npoints, 1))./ ...repmat(max(bm)-min(bm), npoints, 1);hold on;scatter3(bm(:, 1), bm(:, 2), bm(:, 3),...10, pcol, 'filled');grid on;hold off;syms Un X S;15010050 o-5050④%泊松过程的模拟、检验及参数估计100n=10;%生成n*n个随机数r=1;%参数temp=0;tem=0;Un=rand(n,1);%共产生n*n个随机数for i=1:1:nX(i)=-log(Un(i))/r;endX=subs(X);for i=1:1:nfor j=1:1:itemp=temp+X(j);endS(i)=temp;temp=0;endS=subs(S);%检验泊松过程使用第四条for i=1:1:ntem=tem+S(i);endsigmaN=tem;T=S(n);alpha=0.05;%置信水平p=sigmaN/T;p1=(1/2)*(n-1.96*(n/3F(1/2));p2=(1/2)*(n+1.96*(n/3F(1/2));c1=subs(p-p1)c2=subs(p-p2)if (c1<=0&c2>=0)|(c1>=0&c2<=0) disp('这是一个泊松过程！') %参数估计使用极大似然估计r_=subs(n/T);if abs(subs(r_-r))<0.1disp('参数估计正确！')disp('参数估计值为：')r_end%绘制轨迹y=0;x=0:0.001:subs(S(1));Plot(x,y)for k=1:1:ny=k;x=subs(S(k)):0.001:subs(S(k+1));hold onplot(x,y)endelsedisp（'这不是一个泊松过程！'）end96 : ---------------------------------------------------------------------------------5 ” ---------4「3 -2 _ _____________________〔. _______________ ―0 ---------------------------------------- > -------------------------------------150 5 10⑤%二维poisson2d lambda=100; nmb=poissrnd(lambda)x=rand(1,nmb); y=rand(1,nmb);gridscatter(x,y,5,5.*rand(1,nmb));1「0.9 --0.8 - *07 - . , *+0.6 - ，05 -0 4 -0.3 - *■0.2 ；*0.1 - ” ,0 ------------- L -------- L --------- L ----------L_0 0 1 0.2 0.3 0 4⑥%三维poisson3d%单位体积的泊松点数强度为lambda=100;nmb=poissrnd(lambda)x=rand(1,nmb);y=rand(1,nmb);z=rand(1,nmb);gridscatter3(x,y,z,5,5.*rand(1,nmb));0.5 0.£07 0.8 0.9 1 lambda(3) Markov过程的模拟结果。

①离散服务系统中的缓冲动力学%离散服务系统中的缓冲动力学m=200;P=0.2；N=zeros(1,m); %初始化缓冲区A=geornd(1-p,1,m); %生成到达序列模型，比如，几何分布for n=2:mN(n)=N(n-1)+A(n)-(N(n-1)+A(n)>=1);end stairs((0:m-1),N);21.8 -1.6 -1.4 -1.2 -② M/M/1模型% [tjump, systsize] = simmm1(n, lambda, mu)% Inputs: n - number of jumps% lambda - arrival intensity% mu - intensity of the service times% Outputs: tjump - cumulative jump times% systsize - system size% set default parameter values if ommitedI: nargin=0nargin=0;if (nargin==0) n=500;£0 80 100 120 140160180 200 1 -0 0 -0.6 -0.4 -0.2 --lambda=0.8;mu=1;endi=0; %initial value, start on level itjump(1)=O; %start at time 0systsize(1)=i; %at time 0: level ifor k=2:nif i==0mutemp=0;elsemutemp=mu;endtime=-log(rand)/(lambda+mutemp); % Inter-step times:%Exp(lambda+mu)-distributedif rand<lambda/(lambda+mutemp)i=i+1; %jump up: a customer arriveselsei=i-1; %jump down: a customer is departingend %ifsystsize(k)=i;%system size at time i tjump(k)=time;end %for in ： nargin 不为0时 nargin=2;if (nargin==0)n=500;lambda=0.8;mu=1;endtjump=cumsum(tjump);stairs(tjump,systsize);%cumulative jump timesi=0; %initial value, start on level itjump(1)=O; %start at time 0systsize(1)=i; %at time 0: level ifor k=2:nif i==0mutemp=0;elsemutemp=mu;endtime=-log(rand)/(lambda+mutemp); % Inter-step times:%Exp(lambda+mu)-distributedif rand<lambda/(lambda+mutemp)i=i+1; %jump up: a customer arriveselsei=i-1; %jump down: a customer is departingend systsize(k)=i;%system size at time i tjump(k)=time;end %for i③ M/D/1系统% function [jumptimes, systsize] = simmd1(tmax, lambda)% SIMMD1 simulate a M/D/1 queueing system. Poisson arrivals % of intensity lambda, deterministic service times S=1.% [jumptimes, systsize] = simmd1(tmax, lambda)%iftjump=cumsum(tjump);stairs(tjump,systsize);%cumulative jump times% Inputs: tmax - simulation interval% lambda - arrival intensity% Outputs: jumptimes - time points of arrivals or departures % systsize - system size in M/D/1 queue% systtime - system times% Authors:% v1.2 07-0ct-02% set default parameter values if ommitedI: nargin=0nargin=0;if (nargin==0)tmax=1500;lambda=0.95;endarrtime=-log(rand)/lambda; % Poisson arrivalsi=1;while (min(arrtime(i,:))<=tmax)arrtime = [arrtime; arrtime(i, :)-log(rand)/lambda];i=i+1;endn=length(arrtime); % arrival times t_1,…t_narrsubtr=arrtime-(0:n-1)'; % t_k-(k-1)arrmatrix=arrsubtr*ones(1,n);deptime=(1:n)+max(triu(arrmatrix)); % departure times%u_k=k+max(t_1,..,t_k- k+1)B=[ones(n,1) arrtime ; -ones(n,1) deptime'];Bsort=sortrows(B,2); % sort jumps in order jumps=Bsort(:,1);jumptimes=[0;Bsort(:,2)];systsize=[O;cumsum(jumps)]; % M/D/1 processsysttime=deptime-arrtime'; % system times% plot a histogram of system timeshist(systtime,30);140n：nargin 不为0% function [jumptimes, systsize] = simmd1(tmax, lambda)% SIMMD1 simulate a M/D/1 queueing system. Poisson arrivals % of intensity lambda, deterministic service times S=1.%% [jumptimes, systsize] = simmd1(tmax, lambda)%% Inputs: tmax - simulation interval% lambda - arrival intensity% Outputs: jumptimes - time points of arrivals or departures% systsize - system size in M/D/1 queuesysttime - system times% Authors:% v1.2 07-0ct-02% set default parameter values if ommitednargin=2;if (nargin==0)tmax=1500;lambda=0.95;endarrtime=-log(rand)/lambda; % Poisson arrivalsi=1;while (min(arrtime(i,:))<=tmax)arrtime = [arrtime; arrtime(i, :)-log(rand)/lambda];i=i+1;endn=length(arrtime); % arrival times t_1,…t_n arrsubtr=arrtime-(0:n-1)'; % t_k-(k-1)arrmatrix=arrsubtr*ones(1,n);deptime=(1:n)+max(triu(arrmatrix)); % departure times%u_k=k+max(t_1,..,t_k-k +1)B=[ones(n,1) arrtime ; -ones(n,1) deptime'];Bsort=sortrows(B,2); % sort jumps in order jumps=Bsort(:,1);jumptimes=[0;Bsort(:,2)];systsize=[O;cumsum(jumps)]; % M/D/1 processsysttime=deptime-arrtime'; % system times% plot a histogram of system timeshist(systtime,30);900 5 10 15 20 25 308070605040和2010④ M/G/infinity 系统%function [jumptimes, systsize] = simmginfty(tmax, lambda)% SIMMGINFTY simulate a M/G/infinity queueing system. Arrivals are% a homogeneous Poisson process of intensity lambda. Service times% Pareto distributed (can be modified).% [jumptimes, systsize] = simmginfty(tmax, lambda) %% Inputs: tmax - simulation interval% lambda - arrival intensity% Outputs: jumptimes - times of state changes in the system% systsize - number of customers in system% See SIMSTMGINFTY , SIMGEOD1, SIMMM1, SIMMD1, SIMMG1.% set default parameter values if ommitedI : nargin=0nargin=0;if (nargin==0)tmax=1500;lambda=1;end% generate Poisson arrivals% the number of points is Poisson-distributednpoints = poissrnd(lambda*tmax);% conditioned that number of points is N,% the points are uniformly distributed if(npoints>0)arrt = sort(rand(npoints, 1)*tmax);elsearrt =[];end% uncomment if not available POISSONRND % generate Poisson arrivals% arrt=-log(rand)/lambda;% i=1;% while (min(arrt(i,:))<=tmax)% arrt = [arrt; arrt(i, :)-log(rand)/lambda]; % i=i+1;% end% npoints=length(arrt); % servt=50.*rand(n,1); s_1,...,s_kalpha = 1.5; % arrival times t_1,...,t_n% uniform service times% Pareto service timesservt = rand A(-1/(alpha-1))-1; % stationary renewal processservt = [servt; rand(npoints-1,1)4(-1/alpha)-1]; % arbitrary choice of mean% Output is system size process N.B = [ones(npoints, 1) arrt; -ones(npoints, 1) dept]; Bsort = sortrows(B, 2);% sort jumps inorderjumps = Bsort(:, 1); jumptimes = [0; Bsort(:, 2)]; systsize = [0; cumsum(jumps)]; % M/G/infinitysystem size% processstairs(jumptimes, systsize); xmax = max(systsize)+5; axis([0 tmax 0 xmax]); gridservt =10.*servt;dept = arrt+servt; % departure timesn : nargin不为0时%function [jumptimes, systsize] = simmginfty(tmax, lambda)% SIMMGINFTY simulate a M/G/infinity queueing system. Arrivals are% a homogeneous Poisson process of intensity lambda. Service times% Pareto distributed (can be modified).% [jumptimes, systsize] = simmginfty(tmax, lambda) %% Inputs: tmax - simulation interval% lambda - arrival intensity% Outputs: jumptimes - times of state changes in the system% systsize - number of customers in system% See SIMSTMGINFTY , SIMGE0D1, SIMMM1, SIMMD1, SIMMG1.% set default parameter values if ommited nargin=2;if (nargin==O)tmax=1500;lambda=1;end% generate Poisson arrivals% the number of points is Poisson-distributednpoints = poissrnd(lambda*tmax);% conditioned that number of points is N,% the points are uniformly distributedif (npoints>0)arrt = sort(rand(npoints, 1)*tmax);elsearrt =[];end% uncomment if not available POISSONRND% generate Poisson arrivals% arrt=-log(rand)/lambda;% i=1;% while (min(arrt(i,:))v=tmax)% arrt = [arrt; arrt(i, :)-log(rand)/lambda];% i=i+1;% end% npoints=length(arrt); % arrival times t_1,...,t_n% servt=50.*rand(n,1); % uniform service times s_1,...,s_kalpha = 1.5; % Pareto service times servt = rand A(-1/(alpha-1))-1; % stationary renewal process servt = [servt; rand(npoints-1,1).A(-1/alpha)-1]; servt =10.*servt; % arbitrary choice of meandept = arrt+servt; % departure times% Output is system size process N.B = [ones(npoints, 1) arrt; -ones(npoints, 1) dept];Bsort = sortrows(B, 2); % sort jumps in orderjumps = Bsort(:, 1);jumptimes = [0; Bsort(:, 2)];systsize = [0; cumsum(jumps)]; % M/G/infinity system size% process stairs(jumptimes, systsize);xmax = max(systsize)+5;axis([0 tmax 0 xmax]);grid实验总结:通过这次实验，加深了我对随机过程这门课程的理解与认识，对各种随机过程的模拟有了更深刻的了解，同时也认识到模拟编程的重要性。