隨機(jī)過(guò)程matlab程序

1、基本操作-5/(4.8+5.32)2area=pi*2.52x1=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.2 3.3 -4.4 5.5arr1(3)arr1(1 4)arr1(1:2:5)arr2=1 2 3; -2 -3 -4

2、;3 4 5arr2(1,:)arr2(:,1:2:3)arr3=1 2 3 4 5 6 7 8arr3(5:end)arr3(end)繪圖x=0:1:10;y=x.2-10*x+15;plot(x,y)x=0:pi/20:2*piy1=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(-x2-y2)為例講解基本操

3、作， 首先需要利用 meshgrid函數(shù)生成 X-Y 平面的網(wǎng)格數(shù)據(jù)，如下所示：xa = -2:0.2:2;ya = xa;x,y = meshgrid(xa,ya);z = x.*exp(-x.2 - y.2);mesh(x,y,z);建立M文件functionfenshu( grade )ifgrade > 95.0disp('The grade is A.');elseifgrade > 86.0disp('The grade is B.');elseifgrade > 76.0disp('The grade is C.'

4、elseifgrade > 66.0disp('The grade is D.'););elsedisp('The grade is F.');endendendendendfunctiony=func(x)ifabs(x)<1y=sqrt(1-x2);elsey=x2-1;endfunctionsumm( n)i = 1;sum = 0;while( i <= n )sum = sum+i;i = i+1;endstr = '? á 1?a o',num2str(sum);disp(str)end求極限syms xl

5、imit(1+x)(1 /x),x,0,'right')求導(dǎo)數(shù)syms x;f=(sin(x)/x);diff(f)diff(log(sin(x)求積分syms x;int(x2*log(x)syms x;int(abs(x-1),0,2)常微分方程求解dsolve('Dy+2*x*y=x*exp(-x2)','x')計(jì)算偏導(dǎo)數(shù)x/(x2 + y2 + z2)(1 /2)diff(x2+y2+z2)(1 /2),x,2)重積分int(int(x*y,y,2*x,x2+1),x,0,1)級(jí)數(shù)syms n;symsum(1/2n,1,inf)Tayl

6、or 展開(kāi)式求 y=exp(x)在 x=0 處的 5 階 Taylor 展開(kāi)式taylor(exp(x),0,6)矩陣求逆A=0 -6 -1; 6 2 -16; -5 20 -10det(A)inv(A)特征值、特征向量和特征多項(xiàng)式A=0 -6 -1; 6 2 -16; -5 20 -10;lambda=eig(A)v,d=eig(A)poly(A)多項(xiàng)式的根與計(jì)算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(

7、size(x);B=x.2+y.2;c=inv(A'*A)*A'*B;r=sqrt(c(3)+c(1)2+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)密度函數(shù)和概率分布設(shè) x b(20,0.1), binopdf(2,20,0.1)分布函數(shù)設(shè) x N(1100,50 2) ， y N(1150,80 2) ，則有normcdf(1000,1100,50)=0.0228 ， 1-0.022

8、8=0.9772normcdf(1000,1150,80)=0.0304, 1-0.0304=0.9696統(tǒng)計(jì)量數(shù)字特征x=29.8 27.6 28.3mean(x)max(x)min(x)std(x)syms p k;Ex=symsum(k*p*(1-p)(k-1),k,1,inf)syms x y;f=x+y;Ex=int(int(x*y*f,y,0,1),0,1)參數(shù)估計(jì)例：對(duì)某型號(hào)的 20 輛汽車記錄其 5L 汽油的行駛里程（公里） ，觀測(cè)數(shù)據(jù)如下：29.827.628.327.930.128.729.928.027.928.728.427.229.528.528.030.029.12

9、9.829.626.9設(shè)行駛里程服從正態(tài)分布，試用最大似然估計(jì)法求總體的均值和方差。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)2m,s,mci,sci=normfit(x,0.5)假設(shè)檢驗(yàn)例：下面列出的是某工廠隨機(jī)選取的20只零部件的裝配時(shí)間（分）：9.810.410.69.69.79.910.

10、911.19.610.210.39.69.911.210.69.810.510.110.59.7設(shè)裝配時(shí)間總體服從正態(tài)分布，標(biāo)準(zhǔn)差為0.4，是否認(rèn)定裝配時(shí)間的均值在0.05 的水平下不小于10。解 ：在正態(tài)總體的方差已知時(shí)MATLAB 均值檢驗(yàn)程序：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， muc

11、i = 10.05287981908398Inf% 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:0.1:6;f1=1/sqrt(2*pi*sigma1)*exp(-(x-mu1).2/(2*sigma12);f2=1/sqrt(2*pi*sigma2)*exp(-(x-mu2).2/(2*sigma22);plot(x,f1,'r-',x,f2,'k-.')rho=(1+s*t)/(sigma1*sigma2);f=1/(2*

12、pi*sigma1*sigma2*sqrt(1-rho2)*exp(-1/(2*(1-rho2)*(x-mu1)2/sigma12-2*rho*(x-mu1)*(y-mu2)/(sigma1*sigma2)+(y-mu2)2/sigma22);ezsurf(f)0.350.30.250.20.150.10.050-6-4-202462+3 x y-y 2 )0.20.150.10.0505500-5y-5x% P34 例%輸出p1 =0.0671p2 =4.5400e-005ans =0.06710.0000% P40例 p3=poisspdf(9,12)% 輸出p3 = 0.0874% P4

13、0 例% 輸出p4 = 6.1442e-006解微分方程% 輸入：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 =exp(-lamda*t), exp(-lamda*t)*t*lamda, 1/2*exp(-lamda*t)*t2*lamda2% P

14、40 泊松過(guò)程仿真% simulate 10 times clear;m=10; lamda=1; x=; for i=1:m s=exprnd(lamda, 'seed',1 );% seed 是用來(lái)控制生成隨機(jī)數(shù)的種子 , 使得生成隨機(jī)數(shù)的個(gè)數(shù)是一樣的 . x=x,exprnd(lamda);t1=cumsum(x); endx',t1'%輸出：ans =0.65090.65092.40613.05700.10023.15720.12293.28000.82334.10330.24634.34961.90746.25700.47836.73531.34478

15、.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<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;endendplot(0:0.1:

16、(t1(m)+1),N,'r-')%輸出：109876543210012345678910% simulate 100 timesclear;m=100; lamda=1; x=;for i=1:ms= rand('seed');x=x,exprnd(lamda);t1=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,

17、m;endendplot(0:0.1:(t1(m)+1),N,'r-')% 輸出：10090807060504030201000102030405060708090100% P48 非齊次泊松過(guò)程仿真% simulate 10 timesclear;m=10; lamda=1; x=;fori=1:ms=rand( 'seed');% exprnd(lamda,'seed',1); set seedsx=x,exprnd(lamda);t1=cumsum(x);endx',t1'N=; T=;fort=0:0.1:(t1(m)+1

18、)T=T,t.3;% time is adjusted, cumulative intensity function is t3.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-')% output ans =0.42200.42203.33233.75430.16353.91780.06833.98610.38754.37360.27744.65100.29694.94790.93595.88380.

19、42246.30621.76508.0712101009908807706605504403302201100100200300400500600700800024681012005x 1010 times simulation100 times simulation% P50 復(fù)合泊松過(guò)程仿真% simulate 100 timesclear;niter=100;% iterate numberlamda=1;% arriving ratet=input('Input a time:', 's')fori=1:niterrand('state'

20、,sum(clock);x=exprnd(lamda);% interval timet1=x;whilet1<tx=x,exprnd(lamda);t1=sum(x);% arriving timeendt1=cumsum(x);y=trnd(4,1,length(t1);% rand(1,length(t1);gamrnd(1,1/2,1,length(t1); frnd(2,10,1,length(t1);t2=cumsum(y);endx',t1',y',t2'X=; m=length(t1);fort=0:0.1:(t1(m)+1)ift<

21、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-')5050454540403535303025252020151510105502040608010012002040608010012000跳躍度服從 0,1均勻分布情形跳躍度服從(1, 1/ 2) 分布情形20151050-50102030405060708090跳躍度服從t （ 10）分布情形%clear;niter=

22、1.0E4;% number of iterationslamda=6;% arriving rate (unit:minute)t=720;% 12 hours=720 minutesabove=repmat(0,1,niter);% set up storagefor i=1:niterrand('state',sum(clock);x=exprnd(lamda);% interval timen=1;while x<tx=x+exprnd(1/lamda);% arriving timeif x>=tn=n;elsen=n+1;endendz=binornd(

23、200,0.5,1,n);y=sum(z);above(i)=sum(y>432000);end% generate n salespro=mean(above)Output: pro =0.3192% Simulate the loss pro. For a Compound Poisson process clear;niter=1.0E3;% number of iterationslamda=1;% arriving ratet=input('Input a time:','s')below=repmat(0,1,niter);% set up s

24、toragefor i=1:niterrand('state',sum(clock);x=exprnd(lamda);% interval timen=1;while x<tx=x+exprnd(lamda);% arriving timeif x>=tn=n;elsen=n+1;endendr=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 bel

25、ow(i)=sum(minX<log(950000);endpro=mean(below)Output: t=50, pro=0.45馬氏鏈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*Pjueduivec2=chushivec0*(P2)% 計(jì)算 1 到 n 步后的分布chushivec0=0 0 1 0 0 0;P=0,1/2,1/2,0,0,0;1/2,0

26、,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(1 6);jueduivec=repmat(t,n 1);for k=1:njueduiveck=chushivec0*(Pk);jueduivec(k,1:6)=jueduiveckend% 比較相鄰的兩行n=70;jueduivecn=chushivec0*(Pn)n=71;jueduivecn=chushivec0*(Pn)% Replace the first distribution, Comparing two

27、 neighbour 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;n=70;jueduivecn=chushivec0*(Pn)n=71;jueduivecn=chushivec0*(Pn)% 賭博問(wèn)題模擬（帶吸收壁的隨機(jī)游走：結(jié)束1次游走所花的時(shí)間及終止?fàn)顟B(tài)）a=5; p=1/2;m=0;whilem<100m=m+1;r

28、=2*binornd(1,p)-1;ifr=-1a=a-1;elsea=a+1;endifa=0|a=10break ;endendm a% 賭博問(wèn)題模擬 （帶吸收壁的隨機(jī)游走： 結(jié)束 N次游走所花的平均時(shí)間及終止?fàn)顟B(tài)分布規(guī)律）% 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=0t1=t1+m; m1=m1+1;elset2=t2+m; m2=m2+1;en

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

30、N);fprintf('The expectation of arriving 0 or 10 and its approximate respectivelyare %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;

31、m1=m1+1;elset2(1,n)=m; m2=m2+1;endendfprintf('The average times of arriving 0 and 10 respectivelyare %d,%d.n',sum(t1,2)/m1,sum(t2,2)/m2);fprintf('The frequencies of arriving 0 and 10 respectively are %d,%d.n',m1/N,m2/N);% verify:fprintf('The probability of arriving 0 and its appr

32、oximate respectivelyare %d,%d.n', (p10*(1-p)5-p5*(1-p)10)/(p5*(p10-(1-p)10), 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-p)5/p5)/(1-(1-p)10/p10),(sum(t1,2)+sum(t2,2)/N);50甲的預(yù)期輸光時(shí)間45 賭博平均持續(xù)時(shí)間403530at0 25a t2015

33、10500.050.10.150.20.250.30.350.40.450.5p%連續(xù)時(shí)間馬爾可夫鏈 通過(guò) Kolmogorov 微分方程求轉(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','D

34、p01=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+lamda)-exp(-t*mu-t*lamda)*lamda/(mu+lamda)*mu)/mup10 =mu/(mu+lamda)-exp(-t*mu-t*lamda

35、)*mu/(mu+lamda)p11 =(lamda*mu/(mu+lamda)+exp(-t*mu-t*lamda)*mu2/(mu+lamda)/mu% BPATH1 Brownian path simulation: forendrandn('state',100)% set the state of randnT = 1; N = 500; dt = T/N;dW = zeros(1,N);% preallocate arrays .W = zeros(1,N);% for efficiencydW(1) = sqrt(dt)*randn;% first approxi

36、mation outside the loop .W(1) = dW(1);% since W(0) = 0 is not allowedfor j = 2:NdW(j) = sqrt(dt)*randn;% general incrementW(j) = W(j-1) + dW(j);endplot(0:dt:T,0,W,'r-')% plot W against txlabel('t' , 'FontSize',16)ylabel('W(t)', 'FontSize',16,'Rotation',0)% BPATH2 Brownian path simulation