随机过程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=510*rand(2， 3）r4=2randn(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.210*x+15； plot（x，y) x=0:pi/20：2pi 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 = xexp（x2y2) 为例讲解基本操作，首先需要利用

3、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(The grade is C.)； else if grade 66.0 disp（The grade is D.）； else

4、disp（The grade is F。）； end end end end endfunction y=func(x) if abs(x）1 y=sqrt(1x2); else y=x2-1; endfunction summ（ n）i = 1； sum = 0； while ( i = n ) sum = sum+i； i = i+1; end str = ，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

5、x；int(x2log（x）)syms x;int(abs(x-1），0，2）常微分方程求解dsolve（Dy+2*xy=x*exp(-x2）,x）计算偏导数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

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:3y=3.03,3。90,4。35,4.50，4.40，4。02,3.26;A=2x， 2*y, ones（size（x); B=x。2+y。2;c=inv(AA）*A*B;r=sqrt(c（3)+c（1)2+c（2)2）例子ezplot（2/3*exp(t）+5/3*exp（2t）,-2/3exp(t)+2/3*exp(2*t），0,1) grid on； axis(

7、0， 12， 0, 5)密度函数和概率分布设x b(20,0。1），binopdf（2,20,0。1)分布函数设x N（1100，502) ， y N（1150,802) ，则有normcdf（1000，1100,50)=0.0228，10。0228=0。9772 normcdf(1000，1150，80）=0.0304, 10.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(xyf

8、，y，0,1）,0,1)参数估计例：对某型号的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=

9、p(1）， sigma2hatmle=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、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, 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(2pi*sigma1）*exp（(xmu1）。2

11、/（2*sigma12）)；f2=1/sqrt（2pisigma2）exp（-（x-mu2）。2/（2*sigma22）)；plot（x,f1，r-,x,f2,k。） rho=（1+st)/(sigma1*sigma2）； f=1/（2pi*sigma1sigma2*sqrt(1-rho2）exp(1/(2*（1-rho2)*（x-mu1）2/sigma122rho（xmu1）(ymu2)/(sigma1*sigma2)+（ymu2）2/sigma22）)；ezsurf（f）% P34 例3。1.1p1=poisscdf（5，10）p2=poisspdf（0,10）p1，p2%输出p1 =0.

12、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.1442e006% P35-37（Th3。1.1） 解微分方程 输入：syms p0 p1 p2 ；S=dsolve（Dp0=-lamda*p0，Dp1=lamdap1+lamda*p0，Dp2=-lamda*p2+lamda*p1,p0(0） = 1，p1(0) = 0，p2(0） = 0）；S.p0，S.p1，S。p2% 输出：ans =exp（lam

13、dat), exp(lamda*t）*tlamda， 1/2*exp(-lamdat）*t2lamda2% P40 泊松过程仿真 simulate 10 timesclear；m=10； lamda=1； x=； for i=1：ms=exprnd(lamda,seed,1）;% seed是用来控制生成随机数的种子, 使得生成随机数的个数是一样的.x=x,exprnd（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

14、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；elseif 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；

15、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,exprnd（lamda)；t1=cumsum（x）;endx,t1 N=；for t=0:0。1：(t1(m)+1）if tt1（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(seed

16、)； 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（i） tt1（i+1) N=N,i； endend if tt1（m） N=N,m; endendplot(T，N,r) % outputans = 0.4220 0。4220 3。3323 3.7543 0.1635 3。9178 0。0683 3.9

17、861 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； arriving ratet=input（Input a time：，s)for i=1:niter rand(state，sum(clock)）; x=exprnd(lamda）; inter

18、val time t1=x； while t1t 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,10，1,length(t1); t2=cumsum(y);endx，t1，y,t2 X=； m=length（t1）;for t=0:0.1：(t1（m)+1) if tt1(1) X=X，0; end for i=1:（m1) if t=t1(i） t4320

19、00)； end pro=mean（above)Output: pro =0.3192% Simulate the loss pro. For a Compound Poisson 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; w

20、hile xt x=x+exprnd（lamda); arriving time if x=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（minXlog（950000）)； end pro=mean（below)Output： t=50, pro=0.45% P75 （Example 5。1。5）

21、 马氏链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，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

22、，0;n=10t=1/6*ones（1 6）；jueduivec=repmat（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 absolutevectors once morechushivec0=1/6 1/6 1/6 1/6 1/6 1/6

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

24、m a 赌博问题模拟（带吸收壁的随机游走：结束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 a0 a10 m=m+1； r=2*binornd(1,p）-1; if r=1 a=a1; 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,

25、%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 approximate respectively are %d，%d.n, 5（10-5

26、）/(2p）， (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 a0 a15 m=m+1; r=2*binornd（1,p)1； if r=-1 a=a1； else a=a+1； endendif a=0 t1(1，n）=m; m1=m1+1；else t2(1，n）=m; m2=m2+1；endendfprintf（The average times of arriving 0 and 10 respectively are d，%d。n，sum（

27、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)5p5*（1p）10)/（p5(p10（1-p)10）), m1/N）; fprintf（The expectation of arriving 0 or 10 and its appr

28、oximate respectively are d,d.n,5/(12p）-10/（1-2*p)*（1-(1-p）5/p5)/(1-（1p)10/p10)， (sum（t1,2）+sum（t2，2）)/N)； 连续时间马尔可夫链 通过Kolmogorov微分方程求转移概率输入：clear;syms p00 p01 p10 p11 lamda mu; P=p00,p01；p10，p11；Q=-lamda，lamda;mu，muP*Q输出：ans = -p00*lamda+p01mu, p00*lamda-p01*mu -p10*lamda+p11mu, p10*lamda-p11mu输入：p0

29、0,p01，p10,p11=dsolve(Dp00=p00*lamda+p01*mu,Dp01=p00*lamda-p01*mu，Dp10=p10*lamda+p11mu,Dp11=p10*lamda-p11mu,p00(0）=1，p01(0）=0，p10(0）=0，p11（0)=1）输出：p00 =mu/（mu+lamda）+exp(-tmu-t*lamda）lamda/（mu+lamda)p01 =(lamdamu/(mu+lamda)-exp(-t*mu-t*lamda）lamda/(mu+lamda)mu）/mup10 =mu/(mu+lamda)exp(-t*mut*lamda)mu

30、/(mu+lamda）p11 =（lamdamu/(mu+lamda）+exp（-tmutlamda)mu2/（mu+lamda）/mu BPATH1 Brownian path simulation： forend randn（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）*ran