《随机过程上机实验报告讲解.doc》由会员分享，可在线阅读，更多相关《随机过程上机实验报告讲解.doc（31页珍藏版）》请在三一办公上搜索。
1、2015-2016第一学期随机过程第二次上机实验报告实验目的：通过随机过程上机实验，熟悉Monte Carlo计算机随机模拟方法，熟悉Matlab的运行环境，了解随机模拟的原理，熟悉随机过程的编码规律即各种随机过程的实现方法，加深对随机过程的理解。上机内容：（1）模拟随机游走。（2）模拟Brown运动的样本轨道。（3）模拟Markov过程。实验步骤：（1）给出随机游走的样本轨道模拟结果，并附带模拟程序。 一维情形%一维简单随机游走%“从0开始，向前跳一步的概率为p，向后跳一步的概率为1p” n=50;p=0.5; y=0 cumsum(2.*(rand(1,n-1)=p)-1); % n步。
2、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:4 z=2.*(rand(2,n)0.5)-1; x=zeros(1,2); cumsum(z); col=colors
3、tr(k); plot(x(:,1),x(:,2),col); hold on end grid %三维随机游走 ranwalk3d p=0.5; n=10000; colorstr=b r g y; for k=1:4 z=2.*(rand(3,n)=p)-1; x=zeros(1,3); cumsum(z); col=colorstr(k); plot3(x(:,1),x(:,2),x(:,3),col); hold on end grid(2) 给出一维，二维Brown运动和Poisson过程的模拟结果，并附带模拟程序，没有结果的也要把程序记录下来。一维Brown% 这是连续情形的对称随
4、机游动,每个增量W(s+t)-W(s)是高斯分布N(0, t)，不相交区间上的增量是独立的。典型的模拟它方法是用离散时间的随机游动来逼近。 n=1000; dt=1; y=0 cumsum(dt0.5.*randn(1,n); % 标准布朗运动。 plot(0:n,y);二维Brownnpoints = 5000; dt = 1; bm = cumsum(zeros(1, 3); dt0.5*randn(npoints-1, 3); figure(1); plot(bm(:, 1), bm(:, 2), k); pcol = (bm-repmat(min(bm), npoints, 1)./
5、. repmat(max(bm)-min(bm), npoints, 1); hold on; scatter(bm(:, 1), bm(:, 2), . 10, pcol, filled); grid on; hold off;三维Brown npoints = 5000; dt = 1; bm = cumsum(zeros(1, 3); dt0.5*randn(npoints-1, 3); figure(1); plot3(bm(:, 1), bm(:, 2), bm(:, 3), k); pcol = (bm-repmat(min(bm), npoints, 1)./ . repmat(
6、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:n X(i)=-log(Un(i)/r;endX=subs(X);for i=1:1:n for j=1:1:i temp=temp+X(j); end S(i)=temp
7、; temp=0;endS=subs(S);%检验泊松过程使用第四条for i=1:1:n tem=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)|(c1=0&c2=0) disp(这是一个泊松过程！) %参数估计使用极大似然估计 r_=subs(n/T); if abs(subs(r_-r)=1); end stairs(0:m-1),N);M/
8、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
9、.8; mu=1; end i=0; %initial value, start on level i tjump(1)=0; %start at time 0 systsize(1)=i; %at time 0: level i for k=2:n if i=0 mutemp=0; else mutemp=mu; end time=-log(rand)/(lambda+mutemp); % Inter-step times: % Exp(lambda+mu)-distributed if randlambda/(lambda+mutemp) i=i+1; %jump up: a custom
10、er arrives else i=i-1; %jump down: a customer is departing end %if systsize(k)=i; %system size at time i tjump(k)=time; end %for i tjump=cumsum(tjump); %cumulative jump times stairs(tjump,systsize); ：nargin不为0时nargin=2; if (nargin=0) n=500; lambda=0.8; mu=1; end i=0; %initial value, start on level i
11、 tjump(1)=0; %start at time 0 systsize(1)=i; %at time 0: level i for k=2:n if i=0 mutemp=0; else mutemp=mu; end time=-log(rand)/(lambda+mutemp); % Inter-step times: % Exp(lambda+mu)-distributed if randlambda/(lambda+mutemp) i=i+1; %jump up: a customer arrives else i=i-1; %jump down: a customer is de
12、parting end %if systsize(k)=i; %system size at time i tjump(k)=time; end %for i tjump=cumsum(tjump); %cumulative jump times stairs(tjump,systsize); M/D/1系统% function jumptimes, systsize = simmd1(tmax, lambda) % SIMMD1 simulate a M/D/1 queueing system. Poisson arrivals % of intensity lambda, determin
13、istic 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 def
14、ault parameter values if ommited ：nargin=0nargin=0;if (nargin=0) tmax=1500; lambda=0.95; end arrtime=-log(rand)/lambda; % Poisson arrivals i=1; while (min(arrtime(i,:)=tmax) arrtime = arrtime; arrtime(i, :)-log(rand)/lambda; i=i+1; end n=length(arrtime); % arrival times t_1,.t_n arrsubtr=arrtime-(0:
15、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 process syst
16、time=deptime-arrtime; % system times % plot a histogram of system times hist(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 = sim
17、md1(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=2;if (narg
18、in=0) tmax=1500; lambda=0.95; end arrtime=-log(rand)/lambda; % Poisson arrivals i=1; while (min(arrtime(i,:)0) arrt = sort(rand(npoints, 1)*tmax); else arrt = ; end % uncomment if not available POISSONRND % generate Poisson arrivals % arrt=-log(rand)/lambda; % i=1; % while (min(arrt(i,:)0) arrt = so
19、rt(rand(npoints, 1)*tmax); else arrt = ; 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
20、(n,1); % uniform service times s_1,.,s_k alpha = 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 mean dept = arrt+servt; % departure times % Output is system size process N.
21、B = ones(npoints, 1) arrt; -ones(npoints, 1) dept; Bsort = sortrows(B, 2); % sort jumps in order jumps = 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 实验总结： 通过这次实验，加深了我对随机过程这门课程的理解与认识，对各种随机过程的模拟有了更深刻的了解，同时也认识到模拟编程的重要性。理论与实践相结合才能够更全面地理解和认识一门学科。
