随机过程matlab程序.docx

随机过程matlab程序.docx

《随机过程matlab程序.docx》由会员分享，可在线阅读，更多相关《随机过程matlab程序.docx（22页珍藏版）》请在冰豆网上搜索。

随机过程matlab程序

基本操作

-5/（4.8+5.32）八2

area=pi*2.5A2

x1=1+1/2+1/3+1/4+1/5+1/6

exp（acos（0.3））

a=[123;456;789]

a=[1:

3,4:

6,7:

9]

a1=[6:

-1:

1]

a=eye（4）a1=eye（2,3）b=zeros（2,10）c=ones（2,10）c1=8*ones（3,5）

d=zeros（3,2,2）；

r1=rand（2,3）

r2=5-10*rand（2,3）

r4=2*randn（2,3）+3

arr1=[1.1-2.23.3-4.45.5]

arr1（3）arr1（[14]）arr1（1:

2:

5）

arr2=[123;-2-3-4;345]

arr2（1,:

）

arr2（:

1:

2:

3）

arr3=[12345678]

arr3（5:

end）arr3（end）

绘图

x=[0:

1:

10];

y=x42-10*x+15;

plot（x,y）

x=0:

pi/20:

2*piy1=sin（x）;y2=cos（x）;

plot（x,y1,'b-'）;

holdon;

plot（x,y2,‘k--'）;

legend（‘sinx',‘cosx'）;x=0:

pi/20:

2*pi;

y=sin（x）;figure

（1）plot（x,y,'r-'）gridon

meshgrid函数生

以二元函数图z=xexp（-xA2-yA2）为例讲解基本操作，首先需要利用

成X-Y平面的网格数据，如下所示：

xa=-2:

0.2:

2;

ya=xa;

[x,y]=meshgrid（xa,ya）;

z=x.*exp（-x.A2-y.A2）;

mesh（x,y,z）;

建立M文件

functionfenshu（grade）

ifgrade>95.0

disp（'ThegradeisA.'）;

else

ifgrade>86.0

disp（'ThegradeisB.'）;

else

ifgrade>76.0

disp（'ThegradeisC.'）;

else

ifgrade>66.0

disp（'ThegradeisD.'）;

else

disp（'ThegradeisF.'）;endendendendendfunctiony=func（x）ifabs（x）<1

y=sqrt（1-xA2）;

elsey=xA2-1;

endfunctionsumm（n）

i=1;sum=0;

while（i<=n）sum=sum+i;

i=i+1;

end

str=['a1a£o',num2str（sum）];

disp（str）

end求极限symsxlimit（（1+x）A（1/x）,x,0,'right'）求导数symsx;f=（sin（x）/x）;

diff（f）diff（log（sin（x）））求积分symsx;int（xA2*log（x））symsx;

int（abs（x-1）,0,2）常微分方程求解dsolve（'Dy+2*x*y=x*exp（-xA2）','x'）计算偏导数

x/（xA2+yA2+zA2）A（1/2）diff（（xA2+yA2+zA2）A（1/2）,x,2）重积分int（int（x*y,y,2*x,xA2+1）,x,0,1）级数symsn;

symsum（1/2An,1,inf）

Taylor展开式

求y=exp（x）在x=0处的5阶Taylor展开式taylor（exp（x）,0,6）

矩阵求逆

A=[0-6-1;62-16;-520-10]

det（A）

inv（A）

特征值、特征向量和特征多项式

A=[0-6-1;62-16;-520-10];

lambda=eig（A）

[v,d]=eig（A）

poly（A）

多项式的根与计算

p=[10-2-5];

r=roots（p）

p2=poly（r）

y1=polyval（p,4）

例子：

x=[-3:

3]'

y=[3.03,3.90,4.35,4.50,4.40,4.02,3.26]';

A=[2*x,2*y,ones（size（x））];

B=x.A2+y.A2;

c=inv（A'*A）*A'*B;

r=sqrt（c（3）+c

（1）A2+c

（2）A2）

例子ezplot（'-2/3*exp（-t）+5/3*exp（2*t）','-2/3*exp（-t）+2/3*exp（2*t）',[0,1]）gridon;axis（[0,12,0,5]）

密度函数和概率分布

设x~b（20,0.1）,

binopdf（2,20,0.1）

分布函数

设x~N（1100,502），y~N（1150,802），则有

normcdf（1000,1100,50）=0.0228，1-0.0228=0.9772

normcdf（1000,1150,80）=0.0304,1-0.0304=0.9696统计量数字特征

x=[29.827.628.3]

mean（x）

max（x）

min（x）

std（x）

symspk;

Ex=symsum（k*p*（1-p）A（k-1）,k,1,inf）

symsxy;

f=x+y;

Ex=int（int（x*y*f,y,0,1）,0,1）

参数估计

例：

对某型号的20辆汽车记录其

5L汽油的行驶里程（

公里）

观测数据如下：

29.827.628.327.9

30.128.729.9

28.0

27.9

28.7

28.427.229.528.5

28.030.029.1

29.8

29.6

26.9

设行驶里程服从正态分布，试用最大似然估计法求总体的均值和方差。

x1=[29.827.628.327.930.128.729.928.027.928.7];

x2=[28.427.229.528.528.030.029.129.829.626.9];

x=[x1x2]';

p=mle（'norm',x）;

muhatmle=p

（1）,

sigma2hatmle=p

（2）^2

[m,s,mci,sci]=normfit（x,0.5）

假设检验

例：

下面列出的是某工厂随机选取的20只零部件的装配时间（分）：

9.810.410.69.69.79.910.911.19.610.2

10.39.69.911.210.69.810.510.110.59.7

设装配时间总体服从正态分布，标准差为0.4，是否认定装配时间的均值在0.05的水平下不小于10。

解：

：

在正态总体的方差已知时MATLAB匀值检验程序：

x1=[9.810.410.69.69.79.910.911.19.610.2];

x2=[10.39.69.911.210.69.810.510.110.59.7];

x=[x1x2]';m=10;sigma=0.4;a=0.05;[h,sig,muci]=ztest（x,m,sigma,a,1）

得到：

h=1，，

%PPT例2一维正态密度与二维正态密度

symsxy;

s=1;t=2;

mu1=0;mu2=0;sigma1=sqrt（（sA2））;sigma2=sqrt（（tA2））;

x=-6:

0.1:

6;

f1=1/sqrt（2*pi*sigma1）*exp（-（x-mu1）.A2/（2*sigma1A2））;

f2=1/sqrt（2*pi*sigma2）*exp（-（x-mu2）.A2/（2*sigma2A2））;

plot（x,f1,'r-',x,f2,'k-.'）

rho=（1+s*t）/（sigma1*sigma2）;

f=1/（2*pi*sigma1*sigma2*sqrt（1-rhoA2））*exp（-1/（2*（1-rhoA2））*（（x-mu1）A2/sigma1A2-2*rho*（x-mu1）*（y-mu2）/（sigma1*sigma2）+（y-mu2）A2/sigma2A2））;

ezsurf（f）

%P34例3.1.1

p1=poisscdf（5,10）

p2=poisspdf（0,10）

[p1,p2]

%输出

p1=0.0671

p2=4.5400e-005

ans=0.06710.0000

%P40例3.2.1p3=poisspdf（9,12）

%输出p3=0.0874

%P40例3.2.2p4=poisspdf（0,12）

%输出

p4=6.1442e-006

%P35-37（Th3.1.1）解微分方程

%输入：

symsp0p1p2;

S=dsolve（'Dp0=-lamda*p0','Dp1=-lamda*p1+lamda*p0','Dp2=-lamda*p2+lamda*p1','p0（0）=1','p1（0）=0','p2（0）=0'）;

[S.p0,S.p1,S.p2]

%输出：

ans=

[exp（-lamda*t）,exp（-lamda*t）*t*lamda,1/2*exp（-lamda*t）*tA2*lamdaA2]

%P40泊松过程仿真

%simulate10times

clear;

m=10;lamda=1;x=[];

fori=1:

m

s=exprnd（lamda,'seed',1）;

%seed是用来控制生成随机数的种子,使得生成随机数的个数是一样的.x=[x,exprnd（lamda）];

t1=cumsum（x）;

end

[x',t1']

%输出：

ans=

0.6509

0.6509

2.4061

3.0570

0.1002

3.1572

0.1229

3.2800

0.8233

4.1033

0.2463

4.3496

1.9074

6.2570

0.4783

6.7353

1.3447

8.0800

0.8082

8.8882

%输入：

N=[];

fort=0:

0.1:

（t1（m）+1）ift

（1）

N=[N,0];

elseift

（2）

N=[N,1];

elseift

N=[N,2];

elseift

N=[N,3];

elseift

N=[N,4];

elseift

N=[N,5];

elseift

N=[N,6];

elseift

N=[N,7];

elseift

N=[N,8];

elseift

N=[N,9];

else

N=[N,10];

end

end

plot（0:

0.1:

（t1（m）+1）,N,'r-'）

%输出：

%simulate100timesclear;

m=100;lamda=1;x=[];fori=1:

ms=rand（'seed'）;x=[x,exprnd（lamda）];t1=cumsum（x）;

end

[x',t1']

N=[];

fort=0:

0.1:

（t1（m）+1）

ift

（1）

N=[N,0];

end

fori=1:

（m-1）

ift>=t1（i）&t

N=[N,i];

end

end

ift>t1（m）

N=[N,m];

endend

plot（0:

0.1:

（t1（m）+1）,N,'r-'）

%输出：

%P48非齐次泊松过程仿真

%simulate10times

clear;

m=10;lamda=1;x=[];

fori=1:

m

setseeds

s=rand（'seed'）;%exprnd（lamda,'seed',1）;x=[x,exprnd（lamda）];

t1=cumsum（x）;

end

[x',t1']

N=[];T=[];

fort=0:

0.1:

（t1（m）+1）

T=[T,t43];%timeisadjusted,cumulativeintensityfunctionistA3.

ift

（1）

N=[N,0];

end

fori=1:

（m-1）

ift>=t1（i）&t

N=[N,i];

end

end

ift>t1（m）N=[N,m];end

end

plot（T,N,

'r-'）

%output

ans=

0.4220

0.4220

3.3323

3.7543

0.1635

3.9178

0.0683

3.9861

0.3875

4.3736

0.2774

4.6510

0.2969

4.9479

0.9359

5.8838

0.4224

6.3062

1.7650

8.0712

10timessimulation100timessimulation

%P50复合泊松过程仿真

%simulate100times

clear;

niter=100;

lamda=1;

t=input（'Inputatime:

','s'）

fori=1:

niter

rand（'state',sum（clock））;

x=exprnd（lamda）;

t1=x;

whilet1

x=[x,exprnd（lamda）];

t1=sum（x）;

end

t1=cumsum（x）;

y=trnd（4,1,length（t1））;

frnd（2,10,1,length（t1）））;

t2=cumsum（y）;

end

[x',t1',y',t2']

X=[];m=length（t1）;

fort=0:

0.1:

（t1（m）+1）

ift

（1）

X=[X,0];

end

fori=1:

（m-1）

ift>=t1（i）&t

X=[X,t2（i）];

end

end

ift>t1（m）

X=[X,t2（m）];

end

end

plot（0:

0.1:

（t1（m）+1）,X,'r-'）

%iteratenumber

%arrivingrate

%intervaltime

%arrivingtime

%rand（1,length（t1））;gamrnd（1,1/2,1,length（t1）））;

跳跃度服从[0,1]均匀分布情形

跳跃度服从（1,1/2）分布情形

跳跃度服从t（10）分布情形

%%Simulatetheprobability

thatsalesrevenuefallsinsomeinterval.（e.g.example3.3.6

inteachingmaterial）clear;

niter=1.0E4;

lamda=6;

t=720;above=repmat（0,1,niter）;

fori=1:

niter

%numberofiterations

%arrivingrate（unit:

minute）

%12hours=720minutes%setupstorage

rand（'state',sum（clock）

）;

%intervaltime

x=exprnd（lamda）;

n=1;

whilex

x=x+exprnd（1/lamda）;%arrivingtimeifx>=t

n=n;

elsen=n+1;

end

end

z=binornd（200,0.5,1,n）;%generatensales

y=sum（z）;

above（i）=sum（y>432000）;

end

pro=mean（above）

Output:

pro=0.3192

%%Simulatethelosspro.ForaCompoundPoissonprocess

clear;

niter=1.0E3;

lamda=1;

t=input（'Inputatime:

','s'）below=repmat（0,1,niter）;

fori=1:

niter

rand（'state',sum（clock）x=exprnd（lamda）;n=1;

whilex

x=x+exprnd（lamda）;

ifx>=t

%numberofiterations%arrivingrate

%setupstorage

）;

%intervaltime

%arrivingtime

n=n;

else

n=n+1;

end

end

r=normrnd（0.05/253,0.23/sqrt（253）,1,n）;%generatenrandomjumpsy=log（1.0E6）+cumsum（r）;

minX=min（y）;%minmumreturnovernextnjumps

below（i）=sum（minX

endpro=mean（below）

Output:

t=50,pro=0.45

%P75（Example5.1.5）马氏链chushivec0=[001000]

P=[0,1/2,1/2,0,0,0;1/2,0,1/2,0,0,0;1/4,1/4,0,1/4,1/4,0;0,0,1,0,0,0,;0,0,1/2,0,0,1/2;0,0,0,0,1,0]

jueduivec1=chushivec0*P

jueduivec2=chushivecO*（PA2）

%计算1到n步后的分布chushivec0=[001000];

P=[0,1/2,1/2,0,0,0;1/2,0,1/2,0,0,0;1/4,1/4,0,1/4,1/4,0;0,0,1,0,0,0,;0,0,1/2,0,0,1/2;0,0,0,0,1,0];

n=10

t=1/6*ones（[16]）;

jueduivec=repmat（t,[n1]）;

fork=1:

n

jueduiveck=chushivec0*（PAk）;jueduivec（k,1:

6）=jueduiveck

end

%比较相邻的两行

n=70;

jueduivecn=chushivec0*（PAn）

n=71;

jueduivecn=chushivec0*（PAn）

%Replacethefirstdistribution,Comparingtwoneighbourabsolute-vectorsoncemorechushivec0=[1/61/61/61/61/61/6];

P=[0,1/2,1/2,0,0,0;1/2,0,1/2,0,0,0;1/4,1/4,0,1/4,1/4,0;0,0,1,0,0,0,;0,0,1/2,0,0,1/2;0,0,0,0,1,0];

n=70;

jueduivecn=chushivec0*（PAn）

n=71;

jueduivecn=chushivec0*（PAn）

%赌博问题模拟（带吸收壁的随机游走：

结束1次游走所花的时间及终止状态）

a=5;p=1/2;

m=0;

whilem<100

m=m+1;

r=2*binornd（1,p）-1;

ifr==-1

a=a-1;

else

a=a+1;

end

ifa==0|a==10

break;

end

end

[ma]

%赌博问题模拟（带吸收壁的随机游走：

结束N次游走所花的平均时间及终止状态分布规律）

%p=q=1/2

p=1/2;

m1=0;m2=0;N=1000;

t1=0;t2=0;

forn=1:

1:

N

m=0;a=5;whilea>0&a<10

m=m+1;

r=2*binornd（1,p）-1;

ifr==-1

a=a-1;

else

a=a+1;

end

endifa==0

t1=t1+m;m1=m1+1;

else

t2=t2+m;m2=m2+1;

end

end

fprintf（'Theaveragetimesofarriving0and10respectivelyare%d,%d.\n',[t1/m1,t2/m2]）;fprintf（'Thefrequenciesofarriving0and10respectivelyare%d,%d.\n',[m1/N,m2/N]）;

%verify:

fprintf（'Theprobabilityofarriving0anditsapproximaterespectivelyare%d,%d.\n',[5/10,m1/N]）;

fprintf（'Theexpectationofarriving0or10anditsapproximaterespectivelyare%d,%d.\n',[5*（10-5）/（2*p）,（t1+t2）/N]）;

%p~=q

p=1/4;

m1=0;m2=0;N=1000;t1=zeros（1,N）;t2=zeros（1,N）;

forn=1:

1:

N

m=0;a=5;

whilea>0&a<15

m=m+1;

r=2*binornd（1,p）-1;

ifr==-1

a=a-1;

else

a=a+1;

end

endifa==0

t1（1,n）=m;m1=m1+1;

else

t2（1,n）=m;m2=m2+1;

endendfprintf（'Theaveragetimesofarriving0and10respectively

are%d,%d.\n',[sum（t1,2）/m1,sum（t2,2）/m2]）;

fprintf（'Thefrequenciesofarriving