隨機過程matlab程序

1、一口操-5/+八2 area=p 產(chǎn)八2x1=1+1/2+1/3+1/4+1/5+1/6 exp(acos) a=1 2 3;4 5 6;7 8 9 a=1:3,4:6,7:9 a1=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= arr1(3) arr1(1 4) arr1(1:2:5) arr2=1 2 3; -2 -3 -4;3 4 5 arr2(1,:)

2、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 x ', cos x '); x=0:pi/20:2*pi;y=sin(x); figure(1) plot(x,y, 'r-') grid on以二元函數(shù)圖z = xexp(-xA2-

3、yA2)為例講解基本操作，首先需要利用meshgrid函數(shù)生成 X-Y 平面的網(wǎng)格數(shù)據(jù)，如下所示：xa = -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 >);disp( 'The grade is A.' elseif grade >disp( 'The grade is B.' elseif grade >disp(elseif'The grade is

4、 C.'););disp(elsegrade >'The grade is D.');disp('The grade is F.');endendendend endfunction 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;end'? a 1?a £ o' ,num2str(sum);str = disp(str) e

5、nd求極限syms xlimit(1+x)A(1 /x),x,0,'right')求導數(shù)syms x; f=(sin(x)/x);diff(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(-xA2)','x')計算偏導數(shù)x/(x2 + y2 + zA2)A(1 /2) diff(x2+y2+z2)(1 /2),x,2) 重積分in t(i nt(x*y,y,2*x,xA2+1),x,0,1)級數(shù)syms

6、 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)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:3'y=,'A=2*x, 2*y, ones(size(x);B

7、=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)密度函數(shù)和概率分布設(shè) x b(20,binopdf(2,20,分布函數(shù)設(shè) x N(1100,502) ， y N(1150,802) ，則有normcdf(1000,1100,50)= ，=normcdf(1000,1150,80)=, =統(tǒng)計量數(shù)

8、字特征x= mean(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=int(int(x*y*f,y,0,1),0,1)參數(shù)估計例：對某型號的 20 輛汽車記錄其 5L 汽油的行駛里程(公里)觀測數(shù)據(jù)如下：設(shè)行駛里程服從正態(tài)分布，試用最大似然估計法求總體的均值和方差。x1=;x2=;x=x1 x2'p=mle('norm',x);muhatmle=p(1),sigma2hatmle=卩(2)人2m,s,mci,sci=normfit(x,假設(shè)

9、檢驗例：下面列出的是某工廠隨機選取的 20 只零部件的裝配時間(分)10。設(shè)裝配時間總體服從正態(tài)分布，標準差為，是否認定裝配時間的均值在的水平下不小于 解 ： ：在正態(tài)總體的方差已知時 MATLAB 均值檢驗程序：x1=;x2=;x=x1 x2'm=10;sigma=;a=;h,sig,muci=ztest(x,m,sigma,a,1)得到： h =1 ， ，% PPT 例 2 一維正態(tài)密度與二維正態(tài)密度 syms x y;s=1; t=2;mu1=0; mu2=0; sigma1=sqrt(s2); sigma2=sqrt(t2);x=-6:6;f1=1/sqrt(2* pi*sig

10、ma1)*ex p(-(x-mu1)A2/(2*sigma1A2);f2=1/sqrt(2* pi*sigma2)*ex p(-(x-mu2)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)*ex p(-1/(2*(1-rhoA2)*(x-mu1)A2/sigma1A2-2*rho*(x-mu1)*(y-mu2)/(sigma1*sigma2)+(y-mu2)A2/sigma2A2);ezsurf

11、(f)例輸出% P34p1 =p2 = ans =例 輸出% P40 p4 =p3 =% P40例 輸出% 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');% 輸出： ans =ex p(-lamda*t), exp (-lamda*t)*t*

12、lamda, 1/2*ex p(-lamda*t)*tA2*lamdaA2% P40 泊松過程仿真% simulate 10 timesclear;m=10; lamda=1; x=;for i=1:ms=exprnd(lamda, 'seed',1 );% seed 是用來控制生成隨機數(shù)的種子 , 使得生成隨機數(shù)的個數(shù)是一樣的 x=x,exprnd(lamda);t1=cumsum(x);endx',t1'%輸出： ans =%輸入：N=;for t=0:(t1(m)+1)if t<t1(1)N=N,0;elseif t<t1(2)N=N,1;el

13、seif t<t1(3)N=N,2;elseif t<t1(4)N=N,3;elseif t<t1(5)N=N,4;elseif t<t1(6)N=N,5;elseif t<t1(7)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;end end plot(0:(t1(m)+1),N,'r-') %輸出：% simulate 100 times clear;m=100; lamda=1; x=; for i=1:m s= r

14、and( 'seed' ); x=x,exprnd(lamda); t1=cumsum(x);endx',t1'N=;for t=0:(t1(m)+1)if t<t1(1)N=N,0;end for i=1:(m-1)if t>=t1(i) & t<t1(i+1) N=N,i;endendif t>t1(m) N=N,m; end end plot(0:(t1(m)+1),N,'r-')% 輸出：% P48 非齊次泊松過程仿真% simulate 10 times clear;m=10; lamda=1; x=;%

15、 exprnd(lamda,'seed',1 ); set seedsfor i=1:m s=rand( 'seed' ); x=x,exprnd(lamda);t1=cumsum(x);endx',t1'N=; T=;for t=0:(t1(m)+1)T=T,t.A3;% time is adjusted, cumulative inten sity function is t3.if t<t1(1)N=N,0;endfor i=1:(m-1)if t>=t1(i) & t<t1(i+1) N=N,i;endendif

16、 t>t1(m) N=N,m; endend plot(T,N, 'r-' ) % outputans =10 times simulation% P50 復合泊松過程仿真% simulate 100 times100 times simulationclear; niter=100;% iterate numberlamda=1;t=input( 'Input a time:' for i=1:niter% arriving rate, 's' )rand( 'state',sum(clock) x=exprnd(lamd

17、a);t1=x;while t1<t);% interval timex=x,exprnd(lamda);t1=sum(x);% arriving timeendt1=cumsum(x);% rand(1,length(t1);y=trnd(4,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); for t=0:(t1(m)+1) if t<t1(1) X=X,0;endfor i

18、=1:(m-1)if t>=t1(i) & t<t1(i+1) X=X,t2(i);endendif t>t1(m) X=X,t2(m); endendplot(0:(t1(m)+1),X,r-跳躍度服從 0,1均勻分布情形跳躍度服從 (1,1/2)分布情形t ( 10 )分布情形% Simulate the probability that sales revenue falls in some interval. . example in teaching material) clear;niter=; lamda=6;t=720; above=repmat(0,

19、1,niter);for i=1:niterrand( 'state',sum(clock)跳躍度服從% number of iterations% arriving rate (unit:minute)% 12 hours=720 minutes% set up storage);x=exprnd(lamda);n=1;% interval timewhile x<tx=x+exprnd(1/lamda); if x>=t% arriving timeendendn=n;elsen=n+1;z=binornd(200,1,n);% generate n sales

20、y=sum(z);above(i)=sum(y>432000);endpro=mean(above)Output: pro =% Simulate the loss pro. For a Compound Poisson process clear;niter=;lamda=1; t=input('Input a time:','s') below=repmat(0,1,niter);for i=1:niterrand( 'state',sum(clock)x=exprnd(lamda);n=1;while x<t x=x+exprn

21、d(lamda); if x>=t% number of iterations% arriving rate% set up storage);% interval time% arriving timen=n;elsen=n+1;endend% generate n random jumpsr=normrnd253,sqrt(253),1,n); y=log+cumsum(r);minX=min(y);% minmum return over next n jumpsbelow(i)=sum(minX<log(950000);endpro=mean(below)Output: t

22、=50, pro=% P75 (Example 5.1.5) 馬氏鏈 chushivec0=0 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,;0,0,1/2,0,0,1/2;0 ,0,0,0,1,0jueduivec1=chushivec0*P jueduivec2=chushivec0*(卩人2)% 計算 1 到 n 步后的分布 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

23、,1,0,0,0,;0,0,1/2,0,0,1/2;0 ,0,0,0,1,0;n=10t=1/6*ones(1 6); jueduivec=repmat(t,n 1);for k=1:njueduiveck=chushivecO*( Pk);jueduivec(k,1:6)=jueduiveckend % 比較相鄰的兩行n=70;jueduivec n=chushivecO* (Pn)n=71;jueduivec n=chushivecO* (Pn)% Replace the first distribution, Comparing two neighbour absolute-vector

24、s once more chushivecO=1/6 1/6 1/6 1/6 1/6 1/6;P=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,O;n=7O;jueduivec n=chushivecO* (Pn)n=71;jueduivec n=chushivecO* (Pn)1次游走所花的時間及終止狀態(tài)）% 賭博問題模擬(帶吸收壁的隨機游走：結(jié)束a=5; p=1/2;m=O;while m<1OOm=m+1; r=2*binornd(1,p)-1;

25、if r=-1a=a-1;elsea=a+1;endif a=O|a=1O break ;end endm aN次游走所花的平均時間及終止狀態(tài)分布規(guī)律）% 賭博問題模擬(帶吸收壁的隨機游走：結(jié)束% p=q=1/2p=1/2;m1=O; m2=O; N=1OOO;t1=O;t2=O;for n=1:1:Nm=O; a=5;while a>O & a<1Om=m+1; r=2*binornd(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

26、(fprintf(% verify:fprintf('The average times of arriving 0 and 10 respectively 'The frequencies of arriving 0 and 10 respectively are %d,%d.n'are %d,%d.n' ,t1/m1,t2/m2);,m1/N, m2/N);'The probability of arriving 0 and its approximate respectively are %d,%d.n'5/10, m1/N);fprint

27、f( 'The expectation of arriving 0 or 10 and its approximate respectively5*(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);for n=1:1:Nm=0;a=5;while a>0 & a<15are %d,%d.n' ,m=m+1; r=2*binornd(1,p)-1;if r=-1a=a-1;elsea=a+1;endendif a=0t1(1,n)=m; m1=

28、m1+1; elset2(1,n)=m; m2=m2+1;endendfprintf( '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' % verify:fprintf( 'The probability of arriving 0 and its approximate

29、respectively are %d,%d.n' (p 人10*(1- p)A5-pA5*(1- p)A10)/(p 人5*( pA10-(1- p)A10), 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-卩)人10/卩人10),(sum(t1,2)+sum(t2,2)/N);,m1/N, m2/N);%連續(xù)時間馬爾可夫鏈 通過 K

30、olmogorov 微分方程求轉(zhuǎn)移概率 輸入： clear;syms p00 p01 p10 p11 lamda mu;P=p00,p01;p10,p11; Q=-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+p1

31、1*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+lamda)-exp(-t*mu-t*lamda)*lamda/(mu+lamda)*mu)/mu p10 =mu/(mu+lamda)-exp(-t*mu-t*lamda)*mu/(mu+lamda)p11 =(Iamda*mu/(mu+lamda)+ex p(-t

32、*mu-t*lamda)*muA2/(mu+lamda)/mu% BP ATH1 Brow nian path simulatio n: forend% set the state of randnrandn( 'state' ,100)T = 1; N = 500; dt = T/N;% preallocate arrays .% for efficiency% first approximation outside the loop .% since W(0) = 0 is not alloweddW = zeros(1,N);W = zeros(1,N);dW(1) = sqrt(dt)*randn;W(1) = dW(1);for j = 2:N% general incrementdW(j) = sqrt(dt)*randn;W(j) = W(j-1) + dW(j);endplot(0:dt:T,0,W, 'r-' ) xlabel( 't' , 'FontSize' ,16) ylabel( 'W(t)' , 'FontSize' ,16, % BPATH2 Brownian path randn( 'st