随机过程matlab程序资料

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

2、:2:5) 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=x42-10*x+15; plot(x,y) x=0:pi/20:2*pi y1=sin(x);y2=cos(x); plot(x,y1,b-); cos x ); hold on; plot(x,y2, -)k; legend ( sin x x=0:pi/20:2*pi; y=sin(x); figure(1) plot(x,y, r-) grid on 以二元函

3、数图z = xexp(-xA2-yA2)为例讲解基本操作，首先需要利用meshgrid 函数生成 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 文件 function fenshu( grade ) if grade 95.0 disp( The grade is A.); else if grade 86.0 disp( The grade is B.); else if grade 76.0 ); ); ); else

4、if grade 66.0 disp( The grade is D. else disp( The grade is F. end end disp( The grade is C. end end end function y=func(x) if abs(x)1 y=sqrt(1-xA2); else y=xA2-1; end function summ( n) i = 1; sum = 0; while ( i = n ) sum = sum+i; i = i+1; end str = ? a 1?a o ,num2str(sum); disp(str) end 求极限 syms x

5、limit(1+x)A(1 /x),x,0,right) 求导数 syms x; f=(sin(x)/x); diff(f) diff(log(sin(x) 求积分 syms x; in t(xA2*log(x) syms x; 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) 级数 syms n; symsum(1/2An,1,inf

6、) Taylor 展开式 求y=exp(x)在x=0处的5阶Taylor展开式 taylor(exp(x),0,6) 矩阵求逆 A=0 -6 -1; 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=1 0 -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*

7、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) grid on; axis(0, 12, 0, 5) 密度函数和概率分布 设 x b(20,0.1), binopdf(2,20,0.1) 分布函数 22 设 x N(1100,502) ， y N(1150,802) ，则有 normcdf(1000,1100,50)=0.0228 ，1-0.0228=0

8、.9772 normcdf(1000,1150,80)=0.0304, 1-0.0304=0.9696 统计量数字特征 x=29.8 27.6 28.3 mean(x) max(x) min(x) std(x) syms p k; Ex=symsum(k*p*(1-p)A(k-1),k,1,inf) syms x y; 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.7 28.427.229.528

9、.528.030.029.129.829.626.9 设行驶里程服从正态分布，试用最大似然估计法求总体的均值和方差。 x1=29.8 27.6 28.3 27.9 30.1 28.7 29.9 28.0 27.9 28.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=p(2)2 m,s,mci,sci=normfit(x,0.5) 假设检验 例：下面列出的是某工厂随机选取的 20 只零部件的装配时间（分） 9.8 10.4 1

10、0.6 9.6 9.7 10.3 9.6 9.9 11.2 10.6 设装配时间总体服从正态分布， 9.9 10.9 11.1 9.6 10.2 9.8 10.5 10.1 10.5 9.7 0.05 的水平 标准差为 0.4，是否认定装配时间的均值在 不小于 10。 解 ： ： 在正态总体的方差已知时 MATLAB均值检验程序: 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,si

11、g,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=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); p

12、lot(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/sigma1A 2-2*rho*(x-mu1)*(y-mu2)/(sigma1*sigma2)+(y-mu2)A2/sigma2A2); ezsurf(f) 0.35 ,1C1 0.3 -r 0.25 $ fil 0.2 0.15 f1 0.1- 0.05- 0 rrIr -6-4-20246 44798133900177/281474976710656 ex

13、p(-5/2 x 2+3 x y-y 2) 0.2 . % P34 例 3.1.1 p1=poisscdf(5,10) p2=poisspdf(0,10) p1,p2 %输出 p1 =0.0671 p2 =4.5400e-005 ans =0.0671 0.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) 解微分方程 % 输入： syms p0 p1 p2 ; S=dsolve(Dp0=-lam

14、da*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 泊松过程仿真 % simulate 10 times clear; m=10; lamda=1; x=; for i=1:m s=exprnd(lamda, seed,1 ); % seed 是用来控制生成随机数的种子 , 使得

15、生成随机数的个数是一样的 . 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=; 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; el

16、seif 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; else N=N,10; end end plot(0:0.1:(t1(m)+1),N,r-) %输出： 10 9 8 7 6 5 4 3 2 1 1 2 3 4 5 6 7 8 0 910 % simulate 100 times clear; m=100; lamda=1; x=; for i=1:m

17、s= rand( seed); x=x,expr nd(lamda); t仁cumsum(x); end x,t1 N=； for t=0:0.1:(t1(m)+1) if t=t1(i) end end plot(0:0.1:(t1(m)+1),N,r-) %输出: 100 |rE-，tL- 90 _ 80 - .1 70 - 60 L 50 - 40 . 30 - 20 - 10 100 0 |jI 0102030405060708090 % P48非齐次泊松过程仿真 % simulate 10 times clear; m=10; lamda=1; x=; for i=1:m set s

18、eeds s=rand( seed );% exprnd(lamda,seed,1 ); x=x,expr nd(lamda); t仁cumsum(x); end x,t1 N=； T=; for t=0:0.1:(t1(m)+1) T=T,t.A3;% time is adjusted, cumulative inten sity fun cti on is tA3. if t=t1(i) 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.

19、3736 0.2774 4.6510 0.2969 4.9479 0.9359 5.8838 0.4224 6.3062 1.7650 8.0712 10 9 8 100 1 r r 90 80 - - 70 - 60 - 50 - - - - 40 d J - - 30 - - - 20 - - 10 - ” a 0 IT I 7 6 5 4 3 2 1 0 10 times simulati on 100 times simulatio n % P50复合泊松过程仿真 % simulate 100 times clear; niter=100; lamda=1; t=in put( In

20、put a time: for i=1: ni ter % iterate n umber % arrivi ng rate ,s) rand( state,sum(clock) ); x=expr nd(lamda); t1=x; while t1t % in terval time x=x,expr nd(lamda); t1=sum(x); % arrivi ng time end t1=cumsum(x); y=trnd(4,1,le ngth(t1); gamrnd(1,1/2,1,length(t1); frnd(2,10,1,length(t1); t2=cumsum(y); %

21、 ran d(1,le ngth(t1); end x,t1,y,t2 X=; m=le ngth(t1); for t=0:0.1:(t1(m)+1) if t=t1(i) end r-) end plot(0:0.1:(t1(m)+1),X, 50 50 45 40 35 30 25 20 15 10 5 002040 60 80 100 120 跳跃度服从0,1均匀分布情形 跳跃度服从丨(1,1/2)分布情形 跳跃度服从t (10)分布情形 % Simulate the probability that sales revenue falls in some interval. (e.g

22、. 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 storage for i=1:niter rand( state,sum(clock) x=exprnd(lamda); n=1; ); % interval time while x=t % arriving time n=n;

23、 else n=n+1; end end z=binornd(200,0.5,1,n); y=sum(z); above(i)=sum(y432000); % generate n sales end pro=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 % arr

24、iving rate % set up storage for i=1:niter rand( state,sum(clock) x=exprnd(lamda); n=1; ); % interval time while x=t % arriving time n=n; else n=n+1; end end r=normrnd(0.05/253,0.23/sqrt(253),1,n); % generate n random jumps y=log(1.0E6)+cumsum(r); minX=min(y); % minmum return over next n jumps below(

25、i)=sum(minXlog(950000); end pro=mean(below) Output: t=50, pro=0.45 % P75 (Example 5.1.5) 马氏链 chushivec0=0 0 1 0 0 0 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=0 0 1 0

26、0 0; 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(1 6); jueduivec=repmat(t,n 1); for k=1:n jueduiveck=chushivec0*(PAk); jueduivec(k,1:6)=jueduiveck end % 比较相邻的两行 n=70; jueduivecn=chushivec0*(PAn) n=71; jueduivecn=chushivec0*(PAn) %

27、Replace the first distribution, Comparing two neighbour 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*(P n) n=71; jueduivec n=chushivec0*(PA n) % 赌博问题模拟(带吸收壁

28、的随机游走：结束1次游走所花的时间及终止状态) a=5; p=1/2; m=0; while m0 r=2*binornd(1,p)-1; if r=-1 a=a-1; else a=a+1; end end if a=0 t1(1,n)=m; m1=m1+1; else t2(1,n)=m; m2=m2+1; end end fprintf( 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 arriv

29、ing 0 and 10 respectively are %d,%d.n m2/N); % verify: fprintf( The probability of arriving 0 and its approximate respectively are %d,%d.n, (pA10*(1-p)A5-pA5*(1-p)A10)/(pA5*(pA10-(1-p)A10), m1/N); ,m1/N, ,m1/N, fprintf( The expectation of arriving 0 or 10 and its approximate respectively are %d,%d.n

30、,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); 0.05 0.15 0.25 0.35 甲的预期输光时间 赌博平均持续时间 0.450.5 P 40 35 30 25 20 5 o o 5 31 O3L %连续时间马尔可夫链通过Kolmogorov微分方程求转移概率 输入： clear; syms p00 p01 p10 p11 lamda mu; P=p00,p01;p10,p11; Q=-lamda,lamda;mu,-mu P*Q 输出： ans = -p00*lamda+

31、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,Dp1 0=-p10*lamda+p11*mu,Dp11= p10*lamda-p11*mu,p00(0)=1,p01(0)=0,p10(0)=0,p11(0 )=1) 输出： p00 = mu/(mu+lamda)+exp(-t*mu-t*lamda)*lamda/(mu+lamda) p01 = (lamda*mu/(mu+la

32、mda)-exp(-t*mu-t*lamda)*lamda/(mu+lamda)*mu)/mu p10 = mu/(mu+lamda)-exp(-t*mu-t*lamda)*mu/(mu+lamda) p11 = (lamda*mu/(mu+lamda)+exp(-t*mu-t*lamda)*muA2/(mu+lamda)/mu end % set the state of randn % preallocate arrays . % for efficie ncy dW(1) = sqrt(dt)*ra ndn; W(1) = dW(1); for j = 2:N dW(j) = sqrt(

33、dt)*ra ndn; W(j) = W(j-1) + dW(j); end % first approximation outside the loop . % si nee W(0) = 0 is n ot allowed % gen eral in creme nt plot(0:dt:T,0,W, r-) % plot W agai nst t xlabel( t, FontSize ,16) ylabel( W(t), FontSize,16, Rotatio n,0) % BPATH2 Brow nian path simulatio n: vectorized randn( state ,100) T