随机过程matlab程序

1、根本操作-5/(4.8+5.32)八2area=pi*2.5A2x1=1+1/2+1/3+1/4+1/5+1/6exp(acos(0.3)a=1 2 3;4 5 6;7 8 9a=1:3,4:6,7:9a仁6: -1:1a=eye(4) a仁eye(2,3)b=zeros(2,10) c=o nes(2,10) c仁 8*o nes(3,5)d=zeros(3,2,2)；r1=ra nd(2, 3)r2=5-10*rand(2, 3)r4=2*ra ndn (2,3)+3arr仁1.1 -2.2 3.3 -4.4 5.5arr1(3) arr1(1 4) arr1(1:2:5)arr2=1 2

2、 3; -2 -3 -4;3 4 5arr2(1,:)arr2(:,1:2:3)arr3=1 2 3 4 5 6 7 8arr3(5:e nd) arr3(e nd)绘图x=0:1:10;y=x42-10*x+15;Plot(x,y)x=0:pi/20:2*piy1=si n( x);y2=cos(x);plot(x,y1,'b-');hold on;plot(x,y2, -' ')klege nd ( sin x ' , cos x ');x=0:pi/20:2*pi;y=si n( x);figure(1)plot(x,y, 'r-&

3、#39;)grid on以二元函数图z = xexp(-xA2-yA2)为例讲解根本操作，首先需要利用meshgrid函数生成X-Y平面的网格数据，如下所示：xa = -2:022;ya = xa;x,y = meshgrid(xa,ya);z = x.*exp(-x.A2 - y.A2);mesh(x,y,z);建立M文件fun cti onfen shu( grade )if grade > 95.0disp( 'The grade is A.');elseif grade > 86.0disp( 'The grade is B.');elsew

4、ordif grade > 76.0disp('The grade is C.'）；elseifgrade > 66.0disp('The grade is D.'）；elsedisp('The grade is F.'）；endendendendendfun cti ony=fu nc(x)if abs(x)<1y=sqrt(1-xA2);else y=xA2-1;endfun cti onsumm( n)i = 1;sum = 0;while ( i <= n )sum = sum+i;i = i+1;endstr =

5、 '? a 1?a £ o' ,num2str(sum);disp(str) end求极限syms xlimit(1+x)A(l /x),x,0,'right')求导数syms x; f=(si n( x)/x);diff(f)diff(log(si n(x)求积分syms x;init(xA2*log(x)syms x;in t(abs(x-1),0,2)常微分方程求解dsolve('Dy+2*x*y=x*exp(-xA2)','x')计算偏导数x/(xA2 + yA2 + zA2)A(1 /2) diff(xA2+y

6、A2+zA2F(1 /2),x,2)重积分in t(i nt(x*y,y,2*x,xA2+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 2 -16; -5 20 -10 det(A) in v(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)

7、y仁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=in v(A'*A)*A'*B;r=sqrt(c(3)+c(1)A2+c (2)人2)例子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),bin

8、 opdf(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.8 27.6 28.3mean(x)max(x)min (x)std(x)syms p k;Ex=symsum(k*p*(1-p)A(k-1),k,1,i nf) syms x y; f=x+y;Ex=i nt(i nt(x*y*f,y,O,1),O,1)参数估计例：对某型号的20辆

9、汽车记录其5L汽油的行驶里程公里，观测数据如下：29.8 27.628.327.930.128.729.928.027.928.728.4 27.229.528.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),sigma2

10、hatmle=p(2)A2m,s,mci,sci=no rmfit(x,0.5)假设检验例：下面列出的是某工厂随机选取的20只零部件的装配时间分9.810.410.69.69.79.910.911.19.610.2标准差为0.4，是否认定装配时间的均值在0.05的水平10.39.69.911.210.6设装配时间总体服从正态分布，9.810.510.110.59.7F不小于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.

11、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， ， 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)

12、.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-rh°A2)*(x-mu1)A2/sigma1A 2-2*rho*(x-mu1)*(y-mu2)/(sigma1*sigma2)+(y-mu2)A2/sigma2A2);ezsurf(f)0.350.30.250.20.150.10.050-6-4-202462+3 x y-y 2)0.20.15% P34

13、例 3.1.1 p1=poisscdf(5,10) p2=poisspdf(0,10) p1,p2%俞出pl =0.0671p2 =4.5400e-005ans =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)解微分方程%输入：syms p0 p1 p2 ;S=dsolve('Dp0=-lamda*p0','Dp1=-lamda*p1+lamda*p0',&

14、#39;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 timesclear;m=10; lamda=1; x=;for i=1:ms=expr nd(lamda, 'seed',1 );% seed是用来控制生成随机数的

15、种子，使得生成随机数的个数是一样的.x=x,expr nd(lamda);t仁cumsum(x);endx',t1'%输出：ans =0.6509 0.65092.40613.05700.10023.15720.12293.28000.82334.10330.24634.34961.90746.25700.47836.73531.34478.08000.80828.8882%俞入：N=;for t=0:0.1:(t1(m)+1)if t<t1(1)N=N,0;elseif t<t1(2)N=N,1;elseif t<t1 (3)N=N,2;elseif t&l

16、t;t1 (4)N=N,3;elseif t<t1 (5)N=N,4;elseif t<t1(6)N=N,5;elseif t<t1 N=N,6;elseif t<t1(8)N=N,7;elseif t<t1(9)N=N,8;elseif t<t1(10)N=N,9;elseN=N,10;endendplot(0:0.1:(t1(m)+1),N,'r-')%俞出：10987654321123456780910% simulate 100 times clear;m=100; lamda=1; x=; for i=1:ms= rand( 

17、9;seed'); x=x,expr nd(lamda);t仁cumsum(x);endx',t1'N=；for t=0:0.1:(t1(m)+1)if t<t1(1)N=N,0;endfor i=1:(m-1)if t>=t1(i) & t<t1(i+1)N=N,i;endendif t>t1(m)N=N,m;endendplot(0:0.1:(t1(m)+1),N,'r-')%输出：100 |.rE!ETL-90 _80 -70 -60 L50 -40 .30 -20 -10 -1000 |LEI|ZI01020304

18、05060708090% P48非齐次泊松过程仿真% simulate 10 timesclear;m=10; lamda=1; x=;for i=1:mset seedss=rand( 'seed' );% exprnd(lamda,'seed',1 );x=x,expr nd(lamda);t仁cumsum(x);endx',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

19、(1)N=N,0;endfor i=1:(m-1)if t>=t1(i) & t<t1(i+1)N=N,i;endendif t>t1(m)N=N,m;endendplot(T,N, 'r-')% output ans =0.42200.42203.33233.75430.16353.91780.06833.98610.38754.37360.27744.65100.29694.94790.93595.88380.42246.30621.76508.0712109876543218006007000100200300400500010 times si

20、mulati on10090807060504030201040 26 8 10 125x 10100 times simulatio n% P50复合泊松过程仿真 % simulate 100 times clear;niter=100;lamda=1;% iterate n umber % arrivi ng rate,'s')t=in put( 'In put a time:' fori=1: ni terrand( 'state',sum(clock);x=expr nd(lamda);t1=x;% in terval timewhile

21、 t1<tx=x,expr nd(lamda); t1=sum(x);% arrivi ng timeendt1=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);% ran d(1,le ngth(t1);endx',t1',y',t2'X=; m=le ngth(t1);for t=0:0.1:(t1(m)+1)if t<t1(1)X=X,0;endfor i=1:(m-1)if t>=t1(i) &

22、amp; t<t1(i+1)X=X,t2(i);endendif t>t1(m)X=X,t2(m);endendplot(0:0.1:(t1(m)+1),X,'r-')L1ii4545-40-40-35-35-30-j'30-25-25-20-20-r15-15-10口r*10-50ii1K*|50厂rii.i50跳跃度服从0,1均匀分布情形跳跃度服从（1, 1/2）分布情形201510500102030405060708090-5跳跃度服从t 10分布情形% Simulate the probability that sales revenue falls

23、 in some in terval. (e.g. example 3.3.6 in teach ing material) clear;n iter=1.0E4; lamda=6;t=720;above=repmat(0,1, niter);% nu mber of iterati ons% arriving rate (unit:minute)% 12 hours=720 min utes% set up storagefor i=1: ni terrand( 'state',sum(clock)x=expr nd(lamda);n=1;while x<t);% in

24、 terval timex=x+expr nd(1/lamda); if x>=t% arriving timeendendn=n;elsen=n+1;z=bi nornd(200,0.5,1, n);y=sum(z);above(i)=sum(y>432000);end% gen erate n salespro=mea n( above)Output: pro =0.3192% Simulate the loss pro. For a pound Poiss on process clear;n iter=1.0E3;lamda=1;t=in put(' In put

25、a time:','s') below=repmat(0,1, niter);% nu mber of iterati ons% arrivi ng rate% set up storage);% in terval time% arriving timefor i=1: ni terrand( 'state',sum(clock)x=expr nd(lamda);n=1;while x<tx=x+expr nd(lamda); if x>=tn=n;elsen=n+1;endendr=no rmrnd(0.05/253,0.23/sqrt(

26、253),1, n); % gen erate n ran dom jumps y=log(1.0E6)+cumsum(r);mi nX=mi n(y);% minmum return over n ext n jumpsbelow(i)=sum(mi nX<log(950000);endpro=mea n( below)Output: t=50, pro=0.45% P75 (Example 5.1.5) 马氏链chushivecO=O 0 1 0 0 0P=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,

27、;0,0,1/2,0,0,1/2;0,0,0,0,1,0jueduivec仁chushivecO*Pjueduivec2=chushivec0*(PA2)%计算1到n步后的分布chushivecO=O 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;n=10t=1/6* on es(1 6);jueduivec=repmat(t, n 1);for k=1: njueduiveck=chushivec0*(pAk);jueduivec

28、(k,1:6)=jueduiveckend%比拟相邻的两行n=70;jueduive=chushivec0*(PA n)n=71;jueduive=chushivec0*(PA n)%Replace the firstdistributio n,pari ng twon eighbour absolute-vectors once morechushivec0=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

29、;n=70;jueduive=chushivecO*(pA n)n=71;jueduive=chushivec0*(PA n)%赌博问题模拟带吸收壁的随机游走：完毕1次游走所花的时间与终止状态a=5; p=1/2;m=0;while m<100m=m+1;r=2*bi nornd(1,p)-1;if r=-1a=a_1;elsea=a+1;endif a=0|a=10break ;endendm a%赌博问题模拟带吸收壁的随机游走：完毕N次游走所花的平均时间与终止状态分布规律% p=q=1/2p=1/2;m1= 0; m2=0; N=1000;t1=0;t2=0;for n=1:1:Nm

30、=0; a=5;while a>0 & a<10m=m+1;r=2*bi nornd(1,p)-1;if r=-1a=a-1;elsea=a+1;endendif a=0t1=t1+m; m1= m1+1;elset2=t2+m; m2=m2+1;endendfprintf('The average times of arrivi ng 0 and 10 respectivelyare %d,%d.n',t1/m1,t2/m2);fprintf('The frequencies of arriving 0 and 10 respectively ar

31、e %d,%d.n'm2/N);% verify:fprintf('The probability of arriv ing 0 and its approximate respectivelyare %d,%d.n', 5/10, m1/N);fprintf( 'The expectati on of arrivi ng 0 or 10 and its approximate respectively are %d,%d.n', 5*(10-5)/(2*p), (t1+ t2)/N );% p=qp=1/4;m1= 0; m2=0; N=1000;t仁

32、 zeros(1,N);t2=zeros(1,N);for n=1:1:Nm=0;a=5;while a>0 & a<15m=m+1;r=2*bi nornd(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;endendfprintf('The average times of arrivi ng 0 and 10 respectivelyare %d,%d.n',sum(t1,2)/m1,sum(t2,2)/m2);fprintf(&#

33、39;The frequencies of arriving 0 and 10 respectively are %d,%d.n'm2/N);% verify:fprintf('The probability of arriv ing 0 and its approximate respectivelyare %d,%d.n', (pA10*(1-p)A5-pA5*(1-p)A10)/(pA5*(pA10-(1-p)A10), m1/N);fprintf( 'The expectati on of arrivi ng 0 or 10 and its approx

34、imate respectively are %d,%d.n',5/(1-2*p)-10/(1-2*p)*(1-(1-p)人5巾人 5)/(1-(1-p)人10巾A10),m1/N,m1/N,(sum(t1,2)+sum(t2,2)/N);5045403530a0 252015105甲的预期输光时间 赌博平均持续时间00.050.10.150.20.250.30.350.40.450.5P%!续时间马尔可夫链通过Kolmogorov微分方程求转移概率输入:clear;syms p00 p01 p10 p11 lamda mu; P=p00,p01;p10,p11;Q=-lamda,la

35、mda;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','Dp1 0=-p10*lamda+p11*mu','Dp11= p10*lamda-p11*mu','p00(0)=1,p01(0)=0,p10(0)=0,p11(0 )=1')输出：

36、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)/muend% BPATH1 Brow nian path simulatio n: forrandn( 'state'

37、; ,100)T = 1; N = 500; dt = T/N; dW = zeros(1,N);W = zeros(1,N);% set the state of randn% preallocate 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 approximation outside the loop .% si nee W(0) = 0 is n ot allowed% gen eral in creme nt% plot W agai nst tplot(0:dt:T,0,W,'r-')xlabel( 't', 'FontSize' ,16)ylabel( 'W(t)