随机过程matlab程序

1、基本操作-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=1 2 3;4 5 6;7 8 9 a=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)+3 arr1=1.1 -2.2 3.3 -4.4 5.5 arr1(3) arr1(1 4) arr1(1:2:5

2、) arr2=1 2 3; -2 -3 -4;3 4 5 arr2(1,:) arr2(:,1:2:3) arr3=1 2 3 4 5 6 7 8 arr3(5:end) arr3(end)绘图x=0:1:10;y=x.A2-10*x+15;plot(x,y)x=0:pi/20:2*pi y1=sin(x);y2=cos(x);plot(x,y1,b-);hold on; plot(x,y2, -)k; legend ( sin xcos x );x=0:pi/20:2*pi;y=sin(x); figure(1) plot(x,y, r-) grid on以二元函数图z = xexp(-xA

3、2-yA2)为例讲解基本操作，首先需要利用meshgrid函数生成 X-Y 平面的网格数据，如下所示：xa = -2:0.2:2;ya = xa;x,y = meshgrid(xa,ya);z = x.*ex p(-x.A2 - y.A2);mesh(x,y,z);建立 M 文件function fenshu( grade ) if grade 95.0disp( The grade is A. else);if grade 86.0disp( The grade is B. else);The grade is C.);if grade 76.0 disp(elseifdisp(elsegr

4、ade 66.0The grade is D.);disp(The grade is F.);endendendend end function y=func(x) if abs(x)1y=sqrt(1-x2);else y=x2-1;end function summ( n) i = 1;sum = 0;while ( i = n ) sum = sum+i;i = i+1;endstr = disp(str)?a 1?a o ,num2str(sum);end求极限syms xlimit(1+x)A(1 /x),x,0,right)求导数syms x; f=(sin(x)/x); diff

5、(f)diff(log(sin(x)求积分syms x;in t(x2*log(x)syms x;int(abs(x-1),0,2)常微分方程求解dsolve(Dy+2*x*y=x*ex p(帜人2),乂)计算偏导数x/(x2 + yA2 + zA2)A(1 12) diff(xA2+yA2+zA2)A(1 /2),x,2)重积分int(in t(x*y,y,2*x,x2+1),x,0,1)级数syms n; symsum(1/2A n,1,i nf)Taylor 展开式求y=exp(x)在x=0处的5阶Taylor展开式 taylor(exp(x),0,6)矩阵求逆A=0 -6 -1; 6

6、2 -16; -5 20 -10 det(A) inv(A)特征值、特征向量和特征多项式A=0 -6 -1; 6 2 -16; -5 20 -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:3y=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)人2)例子ezplot(-2/3*exp(

7、-t)+5 /3*exp(2*t),-2 /3*exp(-t)+2/3*exp(2*t),0,1) grid on; axis(0, 12, 0, 5)密度函数和概率分布设 x b(20,0.1), binopdf(2,20,0.1)分布函数设 x N(1100,502) ， y N(1 1 50,80 2) ，则有normcdf(1000,1100,50)=0.0228 ，1-0.0228=0.9772normcdf(1000,1150,80)=0.0304, 1-0.0304=0.9696统计量数字特征x=29.8 27.6 28.3mean(x)max(x)min(x)std(x)sym

8、s p k;Ex=symsum(k* p*(1- p)A(k-1),k,1,i nf)syms x y; 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.728.4 27.229.528.528.0 30.029.129.8 29.6 26.9设行驶里程服从正态分布，试用最大似然估计法求总体的均值和方差。x1=29.8 27.6 28.3 27.9 30.1 28.7 29.9 28.0 27.9 2

9、8.7;x2=28.4 27.2 29.5 28.5 28.0 30.0 29.1 29.8 29.6 26.9; x=x1 x2;p=mle(norm,x);muhatmle=p(1), sigma2hatmle=卩(2)人2 m,s,mci,sci=normfit(x,0.5)假设检验例：下面列出的是某工厂随机选取的20 只零部件的装配时间（分）9.8 10.4 10.6 9.6 9.79.910.911.1 9.6 10.210.3 9.6 9.9 11.2 10.6 设装配时间总体服从正态分布， 不小于 10。解 ： ：在正态总体的方差已知时 MATLAB均值检验程序:10.510.1

10、 10.59.79.8标准差为 0.4，是否认定装配时间的均值在0.05 的水平x1=9.8 10.4 10.6 9.6 9.7 9.9 10.9 11.1 9.6 10.2;x2=10.3 9.6 9.9 11.2 10.6 9.8 10.5 10.1 10.5 9.7;x=x1 x2;m=10;sigma=0.4;a=0.05;h,sig,muci=ztest(x,m,sigma,a,1)得到： h =1， sig =0.01267365933873， muci = 10.05287981908398 Inf% PPT例2 一维正态密度与二维正态密度syms x y;s=1;t=2;mu1

11、=0; mu2=0; sigma1=sqrt(s2); sigma2=sqrt(t2);x=-6:0.1:6;f1=1/sqrt(2* pi*sigma1)*ex p(-(x-mu1).A2/(2*sigma1A2); f2=1/sqrt(2* pi*sigma2)*ex p(-(x-mu2).A2/(2*sigma2A2); p Iot(x,f1,r-,x,f2,k-.) rho=(1+s*t)/(sigma1*sigma2);f=1/(2* pi*sigma1*sigma2*sqrt(1-rhoA2)*ex p(-1/(2*(1-rhoA2)*(x-mu1)A2/sigma1A2-2*rh

12、o*(x-mu1)*(y-mu2)/(sigma1*sigma2)+(y-mu2)2/sigma22);ezsurf(f)0.350-60.30.250.20.150.10.05-4-2024644798133900177/281474976710656 exp(-5/2 x 2+3 x y-y 2)5x0.2 . -0.15 .0.1 .0.05 .0 ,-5y% P34 例 3.1.1 p1= poisscdf(5,10) p2=poiss pdf(0,10) p1,p2%输出p1 =0.0671p2 =4.5400e-005 ans =0.0671 0.0000 % P40 例 3.2.

13、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) 解微分方程% 输入：syms p0 p1 p2 ;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 =ex p(-lamda*t), exp(-lamda*t)*t*lamda, 1/2*ex

14、p(-lamda*t)*tA2*lamdaA2 % P40 泊松过程仿真% simulate 10 timesclear;m=10; lamda=1; x=;for i=1:ms=exprnd(lamda, seed,1 );% seed 是用来控制生成随机数的种子 , 使得生成随机数的个数是一样的x=x,exprnd(lamda);t1=cumsum(x);end x,t1 ans =0.65093.05700.65092.4061%输出：0.10023.15720.12293.28000.82334.10330.24634.34961.90746.25700.47836.73531.344

15、78.08000.80828.8882%输入：N=;for t=0:0.1:(t1(m)+1)if tt1(1)N=N,0;elseif tt1(2)N=N,1;elseif tt1(3)N=N,2;elseif tt1(4)N=N,3;elseif tt1(5)N=N,4;elseif tt1(6)N=N,5;elseif tt1(7)N=N,6;elseif tt1(8)N=N,7;elseif tt1(9)N=N,8;elseif tt1(10)N=N,9;elseN=N,10;endendplot(0:0.1:(t1(m)+1),N,r-)%输出：1098765432100123456

16、78910% simulate 100 times clear;m=100; lamda=1; x=;for i=1:m s= rand( seed); x=x,ex prn d(lamda); t仁cumsum(x);endx,t1N=;for t=0:0.1:(t1(m)+1)if t=t1(i) & tt1(m)N=N,m;endendplot(0:0.1:(t1(m)+1),N,r-)%输出:09080706050403020100102030405060708090100% P48非齐次泊松过程仿真% simulate 10 times clear;% exprn d(lamda.s

17、eed,1 ); set seedsm=10; lamda=1; x=; for i=1:ms=rand( seed); x=x,ex prn d(lamda); t仁cumsum(x);end x,t1N=; T=;for t=0:0.1:(t1(m)+1)T=T,t.A3;% time is adjusted, cumulative inten sity fun cti on is t3.if t=t1(i) & tt1(m)N=N,m;endend plot(T,N, r-) % out put ans =0.42200.42203.33233.75430.16353.91780.068

18、33.98610.38754.37360.27744.65100.29694.94790.93595.88380.42246.30621.76508.0712109876543218006007000100200300400500010 times simulati on10090807060504030201024x 10506 8 10100 times simulatio n2% P50复合泊松过程仿真 % simulate 100 times clear;niter=100;lamda=1;% iterate n umber % arrivi ng rate,s)t=inp ut( I

19、nput a time:for i=1: ni terrand( state,sum(clock);x=ex prn d(lamda);t1=x;% in terval timewhile t1tx=x,ex prn d(lamda); t1=sum(x);end% arriv ing timet1=cumsum(x);y=trnd(4,1,le ngth(t1);gamr nd(1,1/2,1,le ngth(t1); frnd(2,10,1,le ngth(t1);t2=cumsum(y);% ran d(1,le ngth(t1);endx,t1,y,t2X=; m=le ngth(t1

20、);for t=0:0.1:(t1(m)+1) if t=t1(i) & tt1(m)X=X,t2(m);endendplot(0:0.1:(t1(m)+1),X,r-)跳跃度服从0,1均匀分布情形跳跃度服从（1, 1/2）分布情形201510500102030405060708090-550Lr45J.45940-40-35-35-30302525E20-20-1-15-15-10105-5-00020406080100120020406080100120跳跃度服从t（ 10）分布情形% Simulate the probability that sales revenue falls in

21、 some interval. (e.g. example 3.3.6 in teaching material) clear;niter=1.0E4; lamda=6;t=720; above=repmat(0,1,niter);% number of iterations% arriving rate (unit:minute) % 12 hours=720 minutes% set up storagefor i=1:niterrand( state,sum(clock) x=exprnd(lamda); n=1;while x=t% arriving timen=n;elsen=n+1

22、;endendz=binornd(200,0.5,1,n);y=sum(z);above(i)=sum(y432000); end% generate n salespro=mean(above)Output: pro =0.3192% Simulate the loss pro. For a Compound Poisson process clear;niter=1.0E3; lamda=1; t=input(Input a time:,s) below=repmat(0,1,niter);% number of iterations% arriving rate% set up stor

23、age);% interval time% arriving timefor i=1:niterrand( state,sum(clock) x=exprnd(lamda); n=1;while x=tn=n; elsen=n+1;endendr=normrnd(O.O5/253,O.23/sqrt(253),1,n); % generate n random jumps y=log(1.OE6)+cumsum(r);minX=min(y);% minmum return over next n jumpsbelow(i)=sum(minXlog(95OOOO);end pro=mean(be

24、low)Output: t=50, pro=0.45% P75 (Example 5.1.5) 马氏链 chushivecO=O O 1 O O OP=O,1/2,1/2,O,O,O;1/2,O,1/2,O,O,O;1/4,1/4,O,1/4,1/4,O;O,O,1,O,O,O,;O,O,1/2,O, O,1/2;O,O,O,O,1,Ojueduivec1=chushivecO*P jueduivec2=chushivec0*(卩人2) % 计算 1 到 n 步后的分布 chushivecO=O O 1 O O O;P=O,1/2,1/2,O,O,O;1/2,O,1/2,O,O,O;1/4,1

25、/4,O,1/4,1/4,O;O,O,1,O,O,O,;O,O,1/2,O, O,1/2;O,O,O,O,1,O;n=1Ot=1/6*ones(1 6); jueduivec=repmat(t,n 1);for k=1:njueduiveck=chushivecO* (P jueduivec(k,1:6)=jueduiveck end % 比较相邻的两行n=70;jueduivec n=chushivecO* (Pn) n=71;jueduivec n=chushivecO* (Pn) % Replace the first distribution, Comparing two neighb

26、our absolute-vectors once more chushivec0=1/6 1/6 1/6 1/6 1/6 1/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;jueduivec n=chushivecO* (Pn)n=71;jueduivec n=chushivecO* (Pn) % 赌博问题模拟（带吸收壁的随机游走：结束 1次游走所花的时间及终止状态） a=5; p=1/2;m=O;while mO & a0 & a1

27、5m=m+1;r=2*binornd(1,p)-1;if r=-1a=a-1;elsea=a+1;endendif a=0t1(1,n)=m; m1=m1+1; elset2(1,n)=m; m2=m2+1;end endfprintf( The average times of arriving 0 and 10 respectively are %d,%d.n ,sum(t1,2)/m1,sum(t2,2)/m2);fprintf( The frequencies of arriving 0 and 10 respectively are %d,%d.n m2/N);,m1/N,m1/N,

28、% verify:fprintf( The probability of arriving 0 and its approximate respectivelyare %d,%d.n, (p 0*(1-卩)人5-p 人5*(1-卩)人10)/(卩人5*(卩人10-(1- p)。)，m1/N);fprintf( The expectation of arriving 0 or 10 and its approximate respectivelyare %d,%d.n,5/(1-2* p)-10/(1-2* p)*(1-(1-卩)人5/卩人5)/(1-(1- plO/pTO),(sum(t1,2

29、)+sum(t2,2)/N);.1通过Kolmogorov微分方程求转移概率%i续时间马尔可夫链输入：clear;syms pOO pO1 p10 p11 lamda mu; P= pOO, p01; p10, p11;Q=-lamda,lamda;mu,-muP*Q输出:ans =-p 00*lamda+p01*mu, p 00*lamda-p01*mu-p 10*lamda+p11*mu, p 10*lamda-p 11*mu输入：p OO, pO1, p1O, p11=dsolve(D pO O=-pOO*lamda+pO1*mu,D pO仁 pOO*lamda-pO1*mu,D p1O

30、=-p1O*lamda+p11*mu,D p1仁 p1O*lamda-p11*mu, pOO(O)=1, pO1(O)=O, p1O(O)=O, p11(O )=1)输出：pOO =mu/(mu+lamda)+ex p(-t*mu-t*lamda)*lamda/(mu+lamda)pO1 =(lamda*mu/(mu+lamda)-ex p(-t*mu-t*lamda)*lamda/(mu+lamda)*mu)/mup10 =mu/(mu+lamda)-ex p(-t*mu-t*lamda)*mu/(mu+lamda)p11 =(lamda*mu/(mu+lamda)+ex p(-t*mu-t

31、*lamda)*muA2/(mu+lamda)/muend% set the state of randn% p reallocate arrays .% for efficie ncydW(1) = sqrt(dt)*ra ndn;W(1) = dW(1);for j = 2:NdW(j) = sqrt(dt)*ra ndn;W(j) = W(j-1) + dW(j); end% first app roximati on outside the loop % si nee W(0) = 0 is n ot allowed% gen eral in creme ntplot(0:dt:T,0,W,r-)% plot W agai nst txlabel( t, FontSize ,16)ylabel( W(t),FontSize,16,Rotatio n,0)% BP ATH2 Brow nian p ath simulatio n: vectorizedrandn( state ,100) T = 1; N = 500; dt = T/N;% set the state of ra ndndW = sqrt(dt)