随机过程matlab程序.docx

随机过程matlab程序.docx

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

随机过程matlab程序

基本操作

-5/（4.8+5.32）^2

area=pi*2.5^2

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=x.^2-10*x+15;

plot（x,y）

x=0:

pi/20:

2*pi

y1=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

以二元函数图z=xexp（-x^2-y^2）为例讲解基本操作，首先需要利用meshgrid函数生成X-Y平面的网格数据，如下所示：

xa=-2:

0.2:

2;

ya=xa;

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

z=x.*exp（-x.^2-y.^2）;

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.'）;

end

end

end

end

end

functiony=func（x）

ifabs（x）<1

y=sqrt（1-x^2）;

elsey=x^2-1;

end

functionsumm（n）

i=1;

sum=0;

while（i<=n）

sum=sum+i;

i=i+1;

end

str=['?

?

?

?

?

á1?

?

a￡o',num2str（sum）];

disp（str）

end

求极限

symsx

limit（（1+x）^（1/x）,x,0,'right'）

求导数

symsx;

f=（sin（x）/x）;

diff（f）

diff（log（sin（x）））

求积分

symsx;

int（x^2*log（x））

symsx;

int（abs（x-1）,0,2）

常微分方程求解

dsolve（'Dy+2*x*y=x*exp（-x^2）','x'）

计算偏导数

x/（x^2+y^2+z^2）^（1/2）

diff（（x^2+y^2+z^2）^（1/2）,x,2）

重积分

int（int（x*y,y,2*x,x^2+1）,x,0,1）

级数

symsn;

symsum（1/2^n,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.^2+y.^2;

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

r=sqrt（c（3）+c

（1）^2+c

（2）^2）

例子

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）^（k-1）,k,1,inf）

symsxy;

f=x+y;

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

参数估计

例：

对某型号的20辆汽车记录其5L汽油的行驶里程（公里），

观测数据如下：

29.8　27.6　28.3　27.9　30.1　28.7　29.9　28.0　27.9　28.7

28.4　27.2　29.5　28.5　28.0　30.0　29.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（（s^2））;sigma2=sqrt（（t^2））;

x=-6:

0.1:

6;

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

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

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

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

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

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.1

p3=poisspdf（9,12）

%输出

p3=0.0874

%P40例3.2.2

p4=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）*t^2*lamda^2]

%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.65090.6509

2.40613.0570

0.10023.1572

0.12293.2800

0.82334.1033

0.24634.3496

1.90746.2570

0.47836.7353

1.34478.0800

0.80828.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-'）

%输出：

%simulate100times

clear;

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

fori=1:

m

s=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];

end

end

plot（0:

0.1:

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

%输出：

%P48非齐次泊松过程仿真

%simulate10times

clear;

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

fori=1:

m

s=rand（'seed'）;%exprnd（lamda,'seed',1）;setseeds

x=[x,exprnd（lamda）];

t1=cumsum（x）;

end

[x',t1']

N=[];T=[];

fort=0:

0.1:

（t1（m）+1）

T=[T,t.^3];%timeisadjusted,cumulativeintensityfunctionist^3.

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.42200.4220

3.33233.7543

0.16353.9178

0.06833.9861

0.38754.3736

0.27744.6510

0.29694.9479

0.93595.8838

0.42246.3062

1.76508.0712

10timessimulation100timessimulation

%P50复合泊松过程仿真

%simulate100times

clear;

niter=100;%iteratenumber

lamda=1;%arrivingrate

t=input（'Inputatime:

','s'）

fori=1:

niter

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

x=exprnd（lamda）;%intervaltime

t1=x;

whilet1

x=[x,exprnd（lamda）];

t1=sum（x）;%arrivingtime

end

t1=cumsum（x）;

y=trnd（4,1,length（t1））;%rand（1,length（t1））;gamrnd（1,1/2,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-'）

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

分布情形

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

%%Simulatetheprobabilitythatsalesrevenuefallsinsomeinterval.（e.g.example3.3.6inteachingmaterial）

clear;

niter=1.0E4;%numberofiterations

lamda=6;%arrivingrate（unit:

minute）

t=720;%12hours=720minutes

above=repmat（0,1,niter）;%setupstorage

fori=1:

niter

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

x=exprnd（lamda）;%intervaltime

n=1;

whilex

x=x+exprnd（1/lamda）;%arrivingtime

ifx>=t

n=n;

else

n=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;%numberofiterations

lamda=1;%arrivingrate

t=input（'Inputatime:

','s'）

below=repmat（0,1,niter）;%setupstorage

fori=1:

niter

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

x=exprnd（lamda）;%intervaltime

n=1;

whilex

x=x+exprnd（lamda）;%arrivingtime

ifx>=t

n=n;

else

n=n+1;

end

end

r=normrnd（0.05/253,0.23/sqrt（253）,1,n）;%generatenrandomjumps

y=log（1.0E6）+cumsum（r）;

minX=min（y）;%minmumreturnovernextnjumps

below（i）=sum（minX

end

pro=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=chushivec0*（P^2）

%计算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*（P^k）;

jueduivec（k,1:

6）=jueduiveck

end

%比较相邻的两行

n=70;

jueduivecn=chushivec0*（P^n）

n=71;

jueduivecn=chushivec0*（P^n）

%Replacethefirstdistribution,Comparingtwoneighbourabsolute-vectorsoncemore

chushivec0=[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*（P^n）

n=71;

jueduivecn=chushivec0*（P^n）

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

结束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

end

ifa==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

end

ifa==0

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

else

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

end

end

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

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

%verify:

fprintf（'Theprobabilityofarriving0anditsapproximaterespectivelyare%d,%d.\n',[（p^10*（1-p）^5-p^5*（1-p）^10）/（p^5*（p^10-（1-p）^10））,m1/N]）;

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

%连续时间马尔可夫链通过Kolmogorov微分方程求转移概率

输入：

clear;

symsp00p01p10p11lamdamu;

P=[p00,p01;p10,p11];

Q=[-lamda,lamda;mu,-mu]

P*Q

输出：

ans=

[-p00*lamda+p01*mu,p00*lamda-p01*mu]

[-p10*lamda+p11*mu,p10*lamda-p11*mu]

输入：

[p00,p01,p10,p11]=dsolve（'Dp00=-p00*lamda+p01*mu','Dp01=p00*lamda-p01*mu','Dp10=-p10*lamda+p11*mu','Dp11=p10*lamda-p11*mu','p00（0）=1,p01（0）=0,p10（0）=0,p11（0）=1'）