随机过程上机实验报告讲解.docx
- 文档编号:13731650
- 上传时间:2023-06-16
- 格式:DOCX
- 页数:27
- 大小:358.09KB
随机过程上机实验报告讲解.docx
《随机过程上机实验报告讲解.docx》由会员分享,可在线阅读,更多相关《随机过程上机实验报告讲解.docx(27页珍藏版)》请在冰点文库上搜索。
随机过程上机实验报告讲解
2015-2016第一学期随机过程第二次上机实验报告
实验目的:
通过随机过程上机实验,熟悉MonteCarlo计算机随机模拟方法,熟悉Matlab的运行环境,了解随机模拟的原理,熟悉随机过程的编码规律即各种随机过程的实现方法,加深对随机过程的理解。
上机内容:
(1)模拟随机游走。
(2)模拟Brown运动的样本轨道。
(3)模拟Markov过程。
实验步骤:
(1)给出随机游走的样本轨道模拟结果,并附带模拟程序。
①一维情形
%一维简单随机游走
%“从0开始,向前跳一步的概率为p,向后跳一步的概率为1-p”
n=50;
p=0.5;
y=[0cumsum(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'];
fork=1:
4
z=2.*(rand(2,n)<0.5)-1;
x=[zeros(1,2);cumsum(z')];
col=colorstr(k);
plot(x(:
1),x(:
2),col);
holdon
end
grid
③%三维随机游走ranwalk3d
p=0.5;
n=10000;
colorstr=['b''r''g''y'];
fork=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);
holdon
end
grid
(2)给出一维,二维Brown运动和Poisson过程的模拟结果,并附带模拟程序,没有结果的也要把程序记录下来。
①一维Brown
%这是连续情形的对称随机游动,每个增量W(s+t)-W(s)是高斯分布N(0,t),不相交区间上的增量是独立的。
典型的模拟它方法是用离散时间的随机游动来逼近。
n=1000;
dt=1;
y=[0cumsum(dt^0.5.*randn(1,n))];%标准布朗运动。
plot(0:
n,y);
②二维Brown
npoints=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);
holdon;
scatter(bm(:
1),bm(:
2),...
10,pcol,'filled');
gridon;
holdoff;
③三维Brown
npoints=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);
holdon;
scatter3(bm(:
1),bm(:
2),bm(:
3),...
10,pcol,'filled');
gridon;
holdoff;
④%泊松过程的模拟、检验及参数估计
symsUnXS;
n=10;%生成n*n个随机数
r=1;%参数
temp=0;
tem=0;
Un=rand(n,1);%共产生n*n个随机数
fori=1:
1:
n
X(i)=-log(Un(i))/r;
end
X=subs(X);
fori=1:
1:
n
forj=1:
1:
i
temp=temp+X(j);
end
S(i)=temp;
temp=0;
end
S=subs(S);
%检验泊松过程使用第四条
fori=1:
1:
n
tem=tem+S(i);
end
sigmaN=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);
ifabs(subs(r_-r))<0.1
disp('参数估计正确!
')
disp('参数估计值为:
')
r_
end
%绘制轨迹
y=0;
x=0:
0.001:
subs(S
(1));
plot(x,y)
fork=1:
1:
n
y=k;
x=subs(S(k)):
0.001:
subs(S(k+1));
holdon
plot(x,y)
end
else
disp('这不是一个泊松过程!
')
end
⑤%二维poisson2d
lambda=100;
nmb=poissrnd(lambda)
x=rand(1,nmb);
y=rand(1,nmb);
grid
scatter(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);
grid
scatter3(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);%生成到达序列模型,比如,几何分布
forn=2:
m
N(n)=N(n-1)+A(n)-(N(n-1)+A(n)>=1);
end
stairs((0:
m-1),N);
②M/M/1模型
%[tjump,systsize]=simmm1(n,lambda,mu)
%Inputs:
n-numberofjumps
%lambda-arrivalintensity
%mu-intensityoftheservicetimes
%Outputs:
tjump-cumulativejumptimes
%systsize-systemsize
%setdefaultparametervaluesifommited
Ⅰ:
nargin=0
nargin=0;
if(nargin==0)
n=500;
lambda=0.8;
mu=1;
end
i=0;%initialvalue,startonleveli
tjump
(1)=0;%startattime0
systsize
(1)=i;%attime0:
leveli
fork=2:
n
ifi==0
mutemp=0;
else
mutemp=mu;
end
time=-log(rand)/(lambda+mutemp);%Inter-steptimes:
%Exp(lambda+mu)-distributed
ifrand i=i+1;%jumpup: acustomerarrives else i=i-1;%jumpdown: acustomerisdeparting end%if systsize(k)=i;%systemsizeattimei tjump(k)=time; end%fori tjump=cumsum(tjump);%cumulativejumptimes stairs(tjump,systsize); Ⅱ: nargin不为0时 nargin=2; if(nargin==0) n=500; lambda=0.8; mu=1; end i=0;%initialvalue,startonleveli tjump (1)=0;%startattime0 systsize (1)=i;%attime0: leveli fork=2: n ifi==0 mutemp=0; else mutemp=mu; end time=-log(rand)/(lambda+mutemp);%Inter-steptimes: %Exp(lambda+mu)-distributed ifrand i=i+1;%jumpup: acustomerarrives else i=i-1;%jumpdown: acustomerisdeparting end%if systsize(k)=i;%systemsizeattimei tjump(k)=time; end%fori tjump=cumsum(tjump);%cumulativejumptimes stairs(tjump,systsize); ③M/D/1系统 %function[jumptimes,systsize]=simmd1(tmax,lambda) %SIMMD1simulateaM/D/1queueingsystem.Poissonarrivals %ofintensitylambda,deterministicservicetimesS=1. %[jumptimes,systsize]=simmd1(tmax,lambda) %Inputs: tmax-simulationinterval %lambda-arrivalintensity %Outputs: jumptimes-timepointsofarrivalsordepartures %systsize-systemsizeinM/D/1queue %systtime-systemtimes %Authors: %v1.207-Oct-02 %setdefaultparametervaluesifommited Ⅰ: nargin=0 nargin=0; if(nargin==0) tmax=1500; lambda=0.95; end arrtime=-log(rand)/lambda;%Poissonarrivals i=1; while(min(arrtime(i,: ))<=tmax) arrtime=[arrtime;arrtime(i,: )-log(rand)/lambda]; i=i+1; end n=length(arrtime);%arrivaltimest_1,...t_n arrsubtr=arrtime-(0: n-1)';%t_k-(k-1) arrmatrix=arrsubtr*ones(1,n); deptime=(1: n)+max(triu(arrmatrix));%departuretimes %u_k=k+max(t_1,..,t_k-k+1) B=[ones(n,1)arrtime;-ones(n,1)deptime']; Bsort=sortrows(B,2);%sortjumpsinorder jumps=Bsort(: 1); jumptimes=[0;Bsort(: 2)]; systsize=[0;cumsum(jumps)];%M/D/1process systtime=deptime-arrtime';%systemtimes %plotahistogramofsystemtimes hist(systtime,30); Ⅱ: nargin不为0 %function[jumptimes,systsize]=simmd1(tmax,lambda) %SIMMD1simulateaM/D/1queueingsystem.Poissonarrivals %ofintensitylambda,deterministicservicetimesS=1. % %[jumptimes,systsize]=simmd1(tmax,lambda) % %Inputs: tmax-simulationinterval %lambda-arrivalintensity %Outputs: jumptimes-timepointsofarrivalsordepartures %systsize-systemsizeinM/D/1queue %systtime-systemtimes %Authors: %v1.207-Oct-02 %setdefaultparametervaluesifommited nargin=2; if(nargin==0) tmax=1500; lambda=0.95; end arrtime=-log(rand)/lambda;%Poissonarrivals i=1; while(min(arrtime(i,: ))<=tmax) arrtime=[arrtime;arrtime(i,: )-log(rand)/lambda]; i=i+1; end n=length(arrtime);%arrivaltimest_1,...t_n arrsubtr=arrtime-(0: n-1)';%t_k-(k-1) arrmatrix=arrsubtr*ones(1,n); deptime=(1: n)+max(triu(arrmatrix));%departuretimes %u_k=k+max(t_1,..,t_k-k+1) B=[ones(n,1)arrtime;-ones(n,1)deptime']; Bsort=sortrows(B,2);%sortjumpsinorder jumps=Bsort(: 1); jumptimes=[0;Bsort(: 2)]; systsize=[0;cumsum(jumps)];%M/D/1process systtime=deptime-arrtime';%systemtimes %plotahistogramofsystemtimes hist(systtime,30); ④M/G/infinity系统 %function[jumptimes,systsize]=simmginfty(tmax,lambda) %SIMMGINFTYsimulateaM/G/infinityqueueingsystem.Arrivalsare %ahomogeneousPoissonprocessofintensitylambda.Servicetimes %Paretodistributed(canbemodified). %[jumptimes,systsize]=simmginfty(tmax,lambda)% %Inputs: tmax-simulationinterval %lambda-arrivalintensity %Outputs: jumptimes-timesofstatechangesinthesystem %systsize-numberofcustomersinsystem %SeeSIMSTMGINFTY,SIMGEOD1,SIMMM1,SIMMD1,SIMMG1. %setdefaultparametervaluesifommited Ⅰ: nargin=0 nargin=0; if(nargin==0) tmax=1500; lambda=1; end %generatePoissonarrivals %thenumberofpointsisPoisson-distributed npoints=poissrnd(lambda*tmax); %conditionedthatnumberofpointsisN, %thepointsareuniformlydistributed if(npoints>0) arrt=sort(rand(npoints,1)*tmax); else arrt=[]; end %uncommentifnotavailablePOISSONRND %generatePoissonarrivals %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);%arrivaltimest_1,...,t_n %servt=50.*rand(n,1);%uniformservicetimess_1,...,s_k alpha=1.5;%Paretoservicetimes servt=rand^(-1/(alpha-1))-1;%stationaryrenewalprocess servt=[servt;rand(npoints-1,1).^(-1/alpha)-1];servt=10.*servt;%arbitrarychoiceofmean dept=arrt+servt;%departuretimes %OutputissystemsizeprocessN. B=[ones(npoints,1)arrt;-ones(npoints,1)dept]; Bsort=sortrows(B,2);%sortjumpsinorder jumps=Bsort(: 1); jumptimes=[0;Bsort(: 2)]; systsize=[0;cumsum(jumps)];%M/G/infinitysystemsize %process stairs(jumptimes,systsize); xmax=max(systsize)+5; axis([0tmax0xmax]); grid Ⅱ: nargin不为0时 %function[jumptimes,systsize]=simmginfty(tmax,lambda) %SIMMGINFTYsimulateaM/G/infinityqueueingsystem.Arrivalsare %ahomogeneousPoissonprocessofintensitylambda.Servicetimes %Paretodistributed(canbemodified). %[jumptimes,systsize]=simmginfty(tmax,lambda)% %Inputs: tmax-simulationinterval %lambda-arrivalintensity %Outputs: jumptimes-timesofstatechangesinthesystem %systsize-numberofcustomersinsystem %SeeSIMSTMGINFTY,SIMGEOD1,SIMMM1,SIMMD1,SIMMG1. %setdefaultparametervaluesifommited nargin=2; if(nargin==0) tmax=1500; lambda=1; end %generatePoissonarrivals %thenumberofpointsisPoisson-distributed npoints=poissrnd(lambda*tmax); %conditionedthatnumberofpointsisN, %thepointsareun
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 随机 过程 上机 实验 报告 讲解