随机过程maab程序

1、基本操作-5/(4.8+5.32)A2area=pi*2.5A2x1=1+1/2+1/3+1/4+1/5+1/6exp(acos(0.3)a=123;456;789a=1:3,4:6,7:9a1=6:-1:1a=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)+3arr1=1.1-2.23.3-4.45.5arr1(3)arr1(14)arr1(1:2:5)arr2=123;-2-3-4;345arr2(1,:)a

2、rr2(:,1:2:3)arr3=12345678arr3(5:end)arr3(end)绘图x=0:1:10;y=xA2-10*x+15;plot(x,y)x=0:pi/20:2*piy1=sin(x);y2=cos(x);plot(x,y1,'b-');holdon;cos x );plot(x,y2,-k);legend(sinxx=0:pi/20:2*pi;y=sin(x);figure(1)plot(x,y,'r-')gridonmeshgrid 函数生以二元函数图z=xexp(-xA2-yA2)为例讲解基本操作，首先需要利用成X-Y平面的网格数据，如

3、下所示：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.0disp('ThegradeisA.');elseifgrade>86.0disp('ThegradeisB.');elseifgrade>76.0disp('ThegradeisC.');elseifgrade>66.0disp('ThegradeisD.');elsedisp(&#

4、39;ThegradeisF.);endendendendendfunctiony=func(x)ifabs(x)<1y=sqrt(1-xA2);elsey=xA2-1;endfunctionsumm(n)i=1;sum=0;while(i<=n)sum=sum+i;i=i+1;endstr='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

5、(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-10det(A)inv(A)特征值、

6、特征向量和特征多项式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)',

7、'-2/3*exp(-t)+2/3*exp(2*t)',0,1)gridon;axis(0,12,0,5)密度函数和概率分布设xb(20,0.1),binopdf(2,20,0.1)分布函数设xN(1100,502)，yN(1150,802)，则有normcdf(1000,1100,50)=0.0228，1-0.0228=0.9772normcdf(1000,1150,80)=0.0304,1-0.0304=0.9696统计量数字特征x=29.827.628.3mean(x)max(x)min(x)std(x)symspk;Ex=symsum(k*p*(1-p)A(k-1),k

8、,1,inf)symsxy;f=x+y;Ex=int(int(x*y*f,y,0,1),0,1)参数估计例：对某型号的20辆汽车记录其5L汽油的行驶里程(公里)，观测数据如下：29.827.628.327.930.128.729.928.027.928.728.427.229.528.528.030.029.129.829.626.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=

9、mle('norm',x);muhatmle=p(1),sigma2hatmle=p(242m,s,mci,sci=normfit(x,0.5)假设检验例：下面列出的是某工厂随机选取的20只零部件的装配时间(分)：9.810.410.69.69.79.910.911.19.610.210.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;x

10、2=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*sigma

11、2A2);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.1p1=poisscdf(5,10)p2=poisspdf(0,10)p1,p2%输出p1=0.0671p2=4.5400e-005ans=0.

12、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%输出

13、：ans=exp(-lamda*t),exp(-lamda*t)*t*lamda,1/2*exp(-lamda*t)*tA2*lamdaA2%P40泊松过程仿真%simulate10timesclear;m=10;lamda=1;x=;fori=1:ms=exprnd(lamda,'seed',1);%seed是用来控制生成随机数的种子,使得生成随机数的个数是一样的.x=x,exprnd(lamda);t1=cumsum(x);endx',t1'%输出：ans=0.65090.65092.40613.05700.10023.15720.12293.28000.8

14、2334.10330.24634.34961.90746.25700.47836.73531.34478.08000.80828.8882%输入：N=;fort=0:0.1:(t1(m)+1)ift<t1(1)N=N,0;elseift<t1(2)N=N,1;elseift<t1(3)N=N,2;elseift<t1(4)N=N,3;elseift<t1(5)N=N,4;elseift<t1(6)N=N,5;elseift<t1(7)N=N,6;elseift<t1(8)N=N,7;elseift<t1(9)N=N,8;elseift<

15、;t1(10)N=N,9;elseN=N,10;endendplot(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);endx',t1'N=;fort=0:0.1:(t1(m)+1)ift<t1(1)N=N,0;endfori=1:(m-1)ift>=t1(i)&t<t1(i+1)N=N,i;endendift>t1(m)N

16、=N,m;endendplot(0:0.1:(t1(m)+1),N,'r-')%输出：%P48非齐次泊松过程仿真%simulate10timesclear;m=10;lamda=1;x=;fori=1:ms=rand('seed');%exprnd(lamda,'seed',1);setseedsx=x,exprnd(lamda);t1=cumsum(x);endx',t1'N=;T=;fort=0:0.1:(t1(m)+1)T=T,t.A3;%timeisadjusted,cumulativeintensityfunctioni

17、stA3.ift<t1(1)N=N,0;endfori=1:(m-1)ift>=t1(i)&t<t1(i+1)N=N,i;endendift>t1(m)N=N,m;endendplot(T,N,'r-')%outputans=0.42200.42203.33233.75430.16353.91780.06833.98610.38754.37360.27744.65100.29694.94790.93595.88380.42246.30621.76508.071210timessimulation100timessimulation%P50复合泊松

18、过程仿真%simulate100timesclear;niter=100;lamda=1;t=input('Inputatime:','s')fori=1:niterrand('state',sum(clock);x=exprnd(lamda);t1=x;whilet1<tx=x,exprnd(lamda);t1=sum(x);endt1=cumsum(x);y=trnd(4,1,length(t1);gamrnd(1,1/2,1,length(t1);frnd(2,10,1,length(t1);t2=cumsum(y);endx'

19、;,t1',y',t2'X=;m=length(t1);fort=0:0.1:(t1(m)+1)ift<t1(1)X=X,0;endfori=1:(m-1)ift>=t1(i)&t<t1(i+1)X=X,t2(i);endendift>t1(m)X=X,t2(m);endendplot(0:0.1:(t1(m)+1),X,'r-')跳跃度服从0,1均匀分布情形%iteratenumber%arrivingrate%intervaltime%arrivingtime%rand(1,length(t1);跳跃度服从(1,1/2

20、)分布情形跳跃度服从t(10)分布情形%Simulatetheprobabilitythatsalesrevenuefallsinsomeinterval.(e.g.example3.3.6inteachingmaterial)clear;niter=1.0E4;%numberofiterationslamda=6;%arrivingrate(unit:minute)t=720;%12hours=720minutesabove=repmat(0,1,niter);%setupstoragefori=1:niterrand('state',sum(clock);x=exprnd(

21、lamda);n=1;whilex<tx=x+exprnd(1/lamda);ifx>=t%intervaltime%arrivingtimen=n;elsen=n+1;endendz=binornd(200,0.5,1,n);%generatensalesy=sum(z);above(i)=sum(y>432000);endpro=mean(above)Output:pro=0.3192%Simulatethelosspro.ForaCompoundPoissonprocessclear;niter=1.0E3;lamda=1;t=input('Inputatime

22、:','s')below=repmat(0,1,niter);fori=1:niterrand('state',sum(clock)x=exprnd(lamda);n=1;whilex<tx=x+exprnd(lamda);ifx>=t%numberofiterations%arrivingrate%setupstorage);%intervaltime%arrivingtimeendendn=n;elsen=n+1;r=normrnd(0.05/253,0.23/sqrt(253),1,n);%generatenrandomjumpsy=l

23、og(1.0E6)+cumsum(r);minX=min(y);%minmumreturnovernextnjumpsbelow(i)=sum(minX<log(950000);endpro=mean(below)Output:t=50,pro=0.45%P75(Example5.1.5)马氏链chushivec0=001000P=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,0jueduivec1=chushivec0*Pjueduivec2=chus

24、hivec0*(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=10t=1/6*ones(16);jueduivec=repmat(t,n1);fork=1:njueduiveck=chushivec0*(P");jueduivec(k,1:6)=jueduiveckend%比较相邻的两行n=70;jueduivecn=chushivec0*(PAn)n=71;jueduiv

25、ecn=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)%赌博问题模拟(带吸收壁的随机游走

26、：结束1次游走所花的时间及终止状态)a=5;p=1/2;m=0;whilem<100m=m+1;r=2*binornd(1,p)-1;ifr=-1a=a-1;elsea=a+1;endifa=0|a=10break;endendma%赌博问题模拟(带吸收壁的随机游走：结束N游走所花的平均时间及终止状态分布规律)%p=q=1/2p=1/2;m1=0;m2=0;N=1000;t1=0;t2=0;forn=1:1:Nm=0;a=5;whilea>0&a<10m=m+1;r=2*binornd(1,p)-1;ifr=-1a=a-1;elsea=a+1;endendifa=0t

27、1=t1+m;m1=m1+1;elset2=t2+m;m2=m2+1;endendfprintf('Theaveragetimesofarriving0and10respectivelyare%d,%d.n',t1/m1,t2/m2);fprintf('Thefrequenciesofarriving0and10respectivelyare%d,%d.n',m1/N,m2/N);%verify:fprintf('Theprobabilityofarriving0anditsapproximaterespectivelyare%d,%d.n',5

28、/10,m1/N);fprintf('Theexpectationofarriving0or10anditsapproximaterespectivelyare%d,%d.n',5*(10-5)/(2*p),(t1+t2)/N);%p=qp=1/4;m1=0;m2=0;N=1000;t1=zeros(1,N);t2=zeros(1,N);forn=1:1:Nm=0;a=5;whilea>0&a<15m=m+1;r=2*binornd(1,p)-1;ifr=-1a=a-1;elsea=a+1;endendifa=0t1(1,n)=m;m1=m1+1;elset

29、2(1,n)=m;m2=m2+1;endendfprintf('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'(pA10*

30、(1-p)A5-pA5*(1-p)A10)/(pA5*(pA10-(1-p)A10),m1/N);fprintf('Theexpectationofarriving0or10anditsapproximaterespectivelyare%d,%d.n',5/(1-2*p)-10/(1-2*p)*(1-(1-p)A5/pA5)/(1-(1-p)A10/pA10),(sum(t1,2)+sum(t2,2)/N);%连续时间马尔可夫链通过Kolmogorov微分方程求转移概率输入：clear;symsp00p01p10p11lamdamu;P=p00,p01;p10,p11;Q=-

31、lamda,lamda;mu,-muP*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')输出

32、：p00=mu/(mu+lamda)+exp(-t*mu-t*lamda)*lamda/(mu+lamda)p01=(lamda*mu/(mu+lamda)-exp(-t*mu-t*lamda)*lamda/(mu+lamda)*mu)/mup10=mu/(mu+lamda)-exp(-t*mu-t*lamda)*mu/(mu+lamda)p11=(lamda*mu/(mu+lamda)+exp(-t*mu-t*lamda)*muA2/(mu+lamda)/muendrandn( 'state' ,100)T = 1; N = 500; dt = T/N;dW = zeros(1,N);W = zeros(1,N);dW(1) = sqrt(dt)*randn;W(1) = dW(1);for j = 2:NdW(j) = sqrt(dt)*randn;W(j) = W(j-1) + dW(j);% set the state of randn%BPATH1Brownianpathsimulation:for%preallocatearrays.%forefficiency%firstapproximationoutsidetheloop.% general increment%sinceW(0)=0isnotallowedendplot(0