随机过程上机实验报告讲解.docx

随机过程上机实验报告讲解.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