随机过程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=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)+3arr1=1.1 -2.2 3.3 -4.4 5.5arr1(3) arr1(1 4) arr1(1:2:5)arr2=1 2 3

2、; -2 -3 -4;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*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以二元函数图 z =

3、xexp(-x2-y2) 为例讲解基本操作，首先需要利用meshgrid函数生成X-Y平面的网格数据，如下所示：xa = -2:0.2:2;ya = xa;x,y = meshgrid(xa,ya);z = x.*exp(-x.2 - y.2);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 disp

4、('The grade is C.'); else if grade > 66.0 disp('The grade is D.'); else disp('The grade is F.'); end end end end endfunction y=func(x) if abs(x)<1 y=sqrt(1-x2); else y=x2-1; endfunction summ( n)i = 1; sum = 0; while ( i <= n ) sum = sum+i; i = i+1; end str = '

5、88;ÆËã½á¹ûΪ£º',num2str(sum); disp(str) end求极限syms xlimit(1+x)(1/x),x,0,'right')求导数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')计算偏导数

6、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)级数syms n;symsum(1/2n,1,inf)Taylor展开式求y=exp(x)在x=0处的5阶Taylor展开式taylor(exp(x),0,6)矩阵求逆A=0 -6 -1; 6 2 -16; -5 20 -10det(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=

7、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(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)密度函数和概率分布设x b(2

8、0,0.1),binopdf(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)(k-1),k,1,inf)syms x y; f=x+y; Ex=int(int(x*y*f,y,0,1),0,1)参数估计例：对某型号的

9、20辆汽车记录其5L汽油的行驶里程（公里），观测数据如下：29.827.628.327.930.128.729.928.027.928.728.427.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), sig

10、ma2hatmle=p(2)2m,s,mci,sci=normfit(x,0.5)假设检验例：下面列出的是某工厂随机选取的20只零部件的装配时间（分）：9.8 10.4 10.6 9.6 9.7 9.9 10.9 11.1 9.6 10.210.3 9.6 9.9 11.2 10.6 9.8 10.5 10.1 10.5 9.7设装配时间总体服从正态分布，标准差为0.4，是否认定装配时间的均值在0.05的水平下不小于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

11、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=0; mu2=0; sigma1=sqrt(s2); sigma2=sqrt(t2); x=-6:0.1:6;f1=1/sqrt(2*pi*sigma1)*exp(-(x-mu1).2/

12、(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*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)% P34 例3.1.1p1=poisscdf(5,10

13、)p2=poisspdf(0,10)p1,p2%输出p1 =0.0671p2 =4.5400e-005ans =0.0671 0.0000 % P40 例3.2.1p3=poisspdf(9,12)% 输出p3 = 0.0874 % P40 例3.2.2p4=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

14、','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% P40 泊松过程仿真% simulate 10 timesclear;m=10; lamda=1; x=; for i=1:ms=exprnd(lamda,'seed',1);% seed是用来控制生成随机数的种子, 使得生成随机数的个数是一样的.x=x,exprnd

15、(lamda);t1=cumsum(x);endx',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 t<t1(1) N=N,0;elseif t<t1(2) N=N,1;elseif t<t1(3) N=N,2;elseif t<t1(4

16、) 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;else N=N,10;endendplot(0:0.1:(t1(m)+1),N,'r-') %输出：% simulate 100 timesclear;m=100; lamda=1; x=; for i=1:ms= rand('seed');x=x,expr

17、nd(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,m;endendplot(0:0.1:(t1(m)+1),N,'r-') % 输出：% P48 非齐次泊松过程仿真% simulate 10 timesclear;m=10; lamda=1; x=; for i=1:ms=rand('see

18、d'); % exprnd(lamda,'seed',1); set seedsx=x,exprnd(lamda);t1=cumsum(x);endx',t1' N=; T=;for t=0:0.1:(t1(m)+1)T=T,t.3; % time is adjusted, cumulative intensity 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; endend if t>t1(m) N=N,m;

19、 endendplot(T,N,'r-') % outputans = 0.4220 0.4220 3.3323 3.7543 0.1635 3.9178 0.0683 3.9861 0.3875 4.3736 0.2774 4.6510 0.2969 4.9479 0.9359 5.8838 0.4224 6.3062 1.7650 8.0712 10 times simulation 100 times simulation% P50 复合泊松过程仿真% simulate 100 timesclear;niter=100; % iterate numberlamda=1;

20、% arriving ratet=input('Input a time:','s')for i=1:niter rand('state',sum(clock); x=exprnd(lamda); % interval time t1=x; while t1<t x=x,exprnd(lamda); t1=sum(x); % arriving time end t1=cumsum(x); y=trnd(4,1,length(t1); % rand(1,length(t1); gamrnd(1,1/2,1,length(t1); frnd(2

21、,10,1,length(t1); t2=cumsum(y);endx',t1',y',t2' X=; m=length(t1);for t=0:0.1:(t1(m)+1) if t<t1(1) X=X,0; end for i=1:(m-1) if t>=t1(i) & t<t1(i+1) X=X,t2(i); end end if t>t1(m) X=X,t2(m); endendplot(0:0.1:(t1(m)+1),X,'r-') 跳跃度服从0,1均匀分布情形 跳跃度服从分布情形跳跃度服从t（10）分布情

22、形 % Simulate the probability that sales revenue falls in some interval. (e.g. example 3.3.6 in teaching material)clear; niter=1.0E4; % number of iterationslamda=6; % arriving rate (unit:minute)t=720; % 12 hours=720 minutesabove=repmat(0,1,niter); % set up storage for i=1:niter rand('state',s

23、um(clock); x=exprnd(lamda); % interval time n=1; while x<t x=x+exprnd(1/lamda); % arriving time if x>=t n=n;else n=n+1; end end z=binornd(200,0.5,1,n); % generate n sales y=sum(z); above(i)=sum(y>432000); end pro=mean(above)Output: pro =0.3192% Simulate the loss pro. For a Compound Poisson

24、processclear; niter=1.0E3; % number of iterationslamda=1; % arriving ratet=input('Input a time:','s') below=repmat(0,1,niter); % set up storage for i=1:niter rand('state',sum(clock); x=exprnd(lamda); % interval time n=1; while x<t x=x+exprnd(lamda); % arriving time if x>

25、;=t 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(i)=sum(minX<log(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,

26、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)% 计算 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,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=repma

27、t(t,n 1); for k=1:n jueduiveck=chushivec0*(Pk); jueduivec(k,1:6)=jueduiveckend% 比较相邻的两行n=70;jueduivecn=chushivec0*(Pn)n=71;jueduivecn=chushivec0*(Pn)% Replace the first distribution, Comparing two 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,

28、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)% 赌博问题模拟（带吸收壁的随机游走：结束1次游走所花的时间及终止状态）a=5; p=1/2;m=0; while m<100 m=m+1; r=2*binornd(1,p)-1; if r=-1 a=a-1; else a=a+1; end if a=0|a=10 break; endendm a % 赌博问题模拟（带吸收壁的随机游走：结束

29、N次游走所花的平均时间及终止状态分布规律） % p=q=1/2 p=1/2;m1=0; m2=0; N=1000;t1=0;t2=0;for n=1:1:N m=0; a=5; while a>0 & a<10 m=m+1; r=2*binornd(1,p)-1; if r=-1 a=a-1; else a=a+1; endendif a=0 t1=t1+m; m1=m1+1;else t2=t2+m; m2=m2+1;endendfprintf('The average times of arriving 0 and 10 respectively are %d,

30、%d.n',t1/m1,t2/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 approximate respectively are %d,%d.n', 5/10, m1/N); fprintf('The expectation of arriving 0 or 10 and its approximat

31、e respectively are %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);for n=1:1:N m=0;a=5; while a>0 & a<15 m=m+1; r=2*binornd(1,p)-1; if r=-1 a=a-1; else a=a+1; endendif a=0 t1(1,n)=m; m1=m1+1;else t2(1,n)=m; m2=m2+1;endendfprintf('The

32、 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',m1/N, m2/N); % verify: fprintf('The probability of arriving 0 and its approximate respectively are %d,%d.n', (p10*(1-p)

33、5-p5*(1-p)10)/(p5*(p10-(1-p)10), m1/N); fprintf('The expectation of arriving 0 or 10 and its approximate respectively are %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); %连续时间马尔可夫链 通过Kolmogorov微分方程求转移概率输入：clear;syms p00 p01 p10 p11 lamda mu; P=p00,p01;p1

34、0,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+p11*mu','Dp11=p10*lamda-p11*mu','p00(0)=1,p01(0)=0,p10(0)=0,p1

35、1(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)*mu/(mu+lamda)p11 =(lamda*mu/(mu+lamda)+exp(-t*mu-t*lamda)*mu2/(mu+lamda)/mu% BPATH1 Brownian path simulation: forend randn(&#

36、39;state',100) % set the state of randnT = 1; N = 500; dt = T/N;dW = zeros(1,N); % preallocate arrays .W = zeros(1,N); % for efficiency dW(1) = sqrt(dt)*randn; % first approximation outside the loop .W(1) = dW(1); % since W(0) = 0 is not allowedfor j = 2:N dW(j) = sqrt(dt)*randn; % general increment W(j) = W(j-1) + dW(j); end plot(0:dt:T,0,W,'r-&#