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)<=p)-1)]; % n步。
plot([0:n-1],y); %画出折线图如下。
%一维随机步长的随机游动%选取任一零均值的分布为步长, 比如，均匀分布。
n=50;x=rand(1,n)-1/2;y=[0 (cumsum(x)-1)];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);hold onendgrid③%三维随机游走ranwalk3d p=0.5;n=10000;colorstr=['b' 'r' 'g' 'y'];for k=1:4z=2.*(rand(3,n)<=p)-1;x=[zeros(1,3); cumsum(z')]; col=colorstr(k);plot3(x(:,1),x(:,2),x(:,3),col);hold onendgrid(2) 给出一维，二维Brown运动和Poisson过程的模拟结果，并附带模拟程序，没有结果的也要把程序记录下来。
①一维Brown% 这是连续情形的对称随机游动,每个增量W(s+t)-W(s)是高斯分布N(0, t)，不相交区间上的增量是独立的。
典型的模拟它方法是用离散时间的随机游动来逼近。
n=1000;dt=1;y=[0 cumsum(dt^0.5.*randn(1,n))]; % 标准布朗运动。
plot(0:n,y);②二维Brownnpoints = 5000;dt = 1;bm = cumsum([zeros(1, 3); dt^0.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;③三维Brownnpoints = 5000;dt = 1;bm = cumsum([zeros(1, 3); dt^0.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;n=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/3)^(1/2)); p2=(1/2)*(n+1.96*(n/3)^(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('这不是一个泊松过程！') end⑤%二维poisson2d lambda=100;nmb=poissrnd(lambda)x=rand(1,nmb);y=rand(1,nmb);gridscatter(x,y,5,5.*rand(1,nmb));⑥%三维poisson3d%单位体积的泊松点数强度为lambda 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));（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);endstairs((0:m-1),N);②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 ommited Ⅰ：nargin=0nargin=0;if (nargin==0)n=500;lambda=0.8;mu=1;endi=0; %initial value, start on level itjump(1)=0; %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 itjump(k)=time;end %for itjump=cumsum(tjump); %cumulative jump times stairs(tjump,systsize);Ⅱ：nargin不为0时nargin=2;if (nargin==0)n=500;lambda=0.8;mu=1;endi=0; %initial value, start on level itjump(1)=0; %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 itjump(k)=time;end %for itjump=cumsum(tjump); %cumulative jump timesstairs(tjump,systsize);③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)% 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-Oct-02% set default parameter values if ommitedⅠ：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=[0;cumsum(jumps)]; % M/D/1 processsysttime=deptime-arrtime'; % system times% plot a histogram of system timeshist(systtime,30);Ⅱ：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 queue% systtime - system times% Authors:% v1.2 07-Oct-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=[0;cumsum(jumps)]; % M/D/1 processsysttime=deptime-arrtime'; % system times% plot a histogram of system timeshist(systtime,30);④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 ommitedⅠ: 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 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,:))<=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^(-1/(alpha-1))-1; % stationary renewal processservt = [servt; rand(npoints-1,1).^(-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Ⅱ: 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, SIMGEOD1, SIMMM1, SIMMD1, SIMMG1.% set default parameter values if ommitednargin=2;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 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,:))<=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^(-1/(alpha-1))-1; % stationary renewal process servt = [servt; rand(npoints-1,1).^(-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实验总结：通过这次实验，加深了我对随机过程这门课程的理解与认识，对各种随机过程的模拟有了更深刻的了解，同时也认识到模拟编程的重要性。
