基本操作-5/(4.8+5.32)^2area=pi*2.5^2x1=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:9]a1=[6: -1:1]a=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.5]arr1(3) arr1([1 4]) arr1(1:2:5)arr2=[1 2 3; -2 -3 -4;3 4 5]arr2(1,:)arr2(:,1:2:3)arr3=[1 2 3 4 5 6 7 8]arr3(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以二元函数图 z = xexp(-x^2-y^2) 为例讲解基本操作，首先需要利用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.0disp('The grade is A.');elseif grade > 86.0disp('The grade is B.');elseif grade > 76.0disp('The grade is C.');elseif grade > 66.0disp('The grade is D.');elsedisp('The grade is F.');endendendendendfunction y=func(x)if abs(x)<1y=sqrt(1-x^2);else y=x^2-1;endfunction summ( n)i = 1;sum = 0;while ( i <= n )sum = sum+i;i = i+1;endstr = ['¼ÆËã½á¹ûΪ£º',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(x^2*log(x))syms x;int(abs(x-1),0,2)常微分方程求解dsolve('Dy+2*x*y=x*exp(-x^2)','x')计算偏导数x/(x^2 + y^2 + z^2)^(1/2)diff((x^2+y^2+z^2)^(1/2),x,2)重积分int(int(x*y,y,2*x,x^2+1),x,0,1)级数syms n;symsum(1/2^n,1,inf)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=[1 0 -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(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(20,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.9772normcdf(1000,1150,80)=0.0304, 1-0.0304=0.9696统计量数字特征x=[29.8 27.6 28.3]mean(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)参数估计例：对某型号的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),sigma2hatmle=p(2)^2[m,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 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((s^2)); sigma2=sqrt((t^2));x=-6:0.1:6;f1=1/sqrt(2*pi*sigma1)*exp(-(x-mu1).^2/(2*sigma1^2));f2=1/sqrt(2*pi*sigma2)*exp(-(x-mu2).^2/(2*sigma2^2));plot(x,f1,'r-',x,f2,'k-.')rho=(1+s*t)/(sigma1*sigma2);f=1/(2*pi*sigma1*sigma2*sqrt(1-rho^2))*exp(-1/(2*(1-rho^2))*((x-mu1)^2/sigma1^2-2*rho*(x-mu1)*(y-mu2)/(sigma1*sigma2)+(y-mu2)^2/sigma2^2));ezsurf(f)-6-4-20246x 44798133900177/281474976710656 exp(-5/2 x 2+3 x y-y 2)y% P34 例3.1.1p1=poisscdf(5,10)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', '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)*t^2*lamda^2]% 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);end[x',t1']%输出：ans =0.6509 0.65092.40613.05700.1002 3.15720.1229 3.28000.8233 4.10330.2463 4.34961.9074 6.25700.4783 6.73531.3447 8.08000.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)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:(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);end[x',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('seed'); % exprnd(lamda,'seed',1); set seedsx=[x,exprnd(lamda)];t1=cumsum(x);end[x',t1']N=[]; T=[];for t=0:0.1:(t1(m)+1)T=[T,t.^3]; % time is adjusted, cumulative intensity function is t^3. 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(T,N,'r-')% outputans =0.4220 0.42203.3323 3.75430.1635 3.91780.0683 3.98610.3875 4.37360.2774 4.65100.2969 4.94790.9359 5.88380.4224 6.30621.7650 8.0712x 10510 times simulation 100 times simulation% P50 复合泊松过程仿真% simulate 100 timesclear;niter=100; % iterate numberlamda=1; % arriving ratet=input('Input a time:','s')for i=1:niterrand('state',sum(clock));x=exprnd(lamda); % interval timet1=x;while t1<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);end[x',t1',y',t2']X=[]; m=length(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) & 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-')跳跃度服从[0,1]均匀分布情形跳跃度服从)2/1,1( 分布情形跳跃度服从t（10）分布情形%% 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 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(200,0.5,1,n); % generate n salesy=sum(z);above(i)=sum(y>432000);endpro=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 storagefor 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 jumpsy=log(1.0E6)+cumsum(r);minX=min(y); % minmum return over next n jumpsbelow(i)=sum(minX<log(950000));endpro=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,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]jueduivec1=chushivec0*Pjueduivec2=chushivec0*(P^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,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*(P^k);jueduivec(k,1:6)=jueduiveckend% 比较相邻的两行n=70;jueduivecn=chushivec0*(P^n)n=71;jueduivecn=chushivec0*(P^n)% 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,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*(P^n)n=71;jueduivecn=chushivec0*(P^n)% 赌博问题模拟（带吸收壁的随机游走：结束1次游走所花的时间及终止状态）a=5; p=1/2;m=0;while m<100m=m+1;r=2*binornd(1,p)-1;if r==-1a=a-1;elsea=a+1;endif a==0|a==10break;endend[m a]% 赌博问题模拟（带吸收壁的随机游走：结束N次游走所花的平均时间及终止状态分布规律）% p=q=1/2p=1/2;m1=0; m2=0; N=1000;t1=0;t2=0;for n=1:1:Nm=0; a=5;while a>0 & a<10m=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('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',[m1/N, m2/N]);% verify:fprintf('The probability of arriving 0 and its approximate respectivelyare %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)/(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<15m=m+1;r=2*binornd(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 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 approximate respectivelyare %d,%d.\n', [(p^10*(1-p)^5-p^5*(1-p)^10)/(p^5*(p^10-(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/p^5)/(1-(1-p)^10/p^10),(sum(t1,2)+sum(t2,2))/N]);pt a 0 t a%连续时间马尔可夫链 通过Kolmogorov 微分方程求转移概率 输入： clear;syms p00 p01 p10 p11 lamda mu; P=[p00,p01;p10,p11]; Q=[-lamda,lamda;mu,-mu] P*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,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 =(lamda*mu/(mu+lamda)+exp(-t*mu-t*lamda)*mu^2/(mu+lamda))/mu% BPATH1 Brownian path simulation: for…endrandn('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 approximation 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: vectorizedrandn('state',100) % set the state of randnT = 1; N = 500; dt = T/N;dW = sqrt(dt)*randn(1,N); % incrementsW = cumsum(dW); % cumulative sumplot([0:dt:T],[0,W],'r-') % plot W against txlabel('t','FontSize',16)ylabel('W(t)','FontSize',16,'Rotation',0)W(t)t%BPATH3 Function along a Brownian pathrandn('state',100) % set the state of randn T = 1; N = 500; dt = T/N; t = [dt:dt:1];M = 1000; % M paths simultaneously dW = sqrt(dt)*randn(M,N); % incrementsW = cumsum(dW,2); % cumulative sumU = exp(repmat(t,[M 1]) + 0.5*W);Umean = mean(U);plot([0,t],[1,Umean],'b-'), hold on% plot mean over M paths plot([0,t],[ones(5,1),U(1:5,:)],'r--'), hold off% plot 5 individual pathsxlabel('t','FontSize',16)ylabel('U(t)','FontSize',16,'Rotation',0,'HorizontalAlignment','right') legend('mean of 1000 paths','5 individual paths',2)aveerr = norm((Umean - exp(9*t/8)),'inf') % sample error% 输出：U(t)taveerr = 0.0504。