《第4MATLAB的数值计算.ppt》由会员分享，可在线阅读，更多相关《第4MATLAB的数值计算.ppt（62页珍藏版）》请在三一文库上搜索。
1、2019/4/3,第1页,第4章 MATLAB 的数值计算,4.1 数值微积分 4.2 矩阵和代数方程 4.3 概率分布和统计分析 4.4 多项式运算和卷积,2019/4/3,第2页,4.1数值微积分,4.1.1 近似数值极限及导数 4.1.2 数值求和与近似数值积分 4.1.3 计算精度可控的数值积分 4.1.4 函数极值的数值求解 4.1.5 常微分方程的数值解,2019/4/3,第3页,在MATLAB数值计算中，既没有专门的求极限指令，也没有专门的求导指令。但MATLAB提供了与“求导”概念有关的“求差分”指令。 dx=diff(X) %计算向量X的前向差分 FX=gradient(F)
2、 %求一元(函数)梯度 FX, FY =gradient(F) %求二元（函数）梯度,对diff而言，当X是向量时，dx= X(2:n)-X(1:n-1) ;当X是矩阵时，dx= X(2:n, :)-X(1:n-1, :) 。 dx的长度比x的长度少1个元素。,4.1.1 近似数值极限及导数,2019/4/3,第4页,dx=diff(X) %计算向量X的前向差分 FX=gradient(F) %求一元(函数)梯度 FX, FY =gradient(F) %求二元（函数）梯度,对gradient而言，当F是向量时，FX(1) = F(2)-F(1), FX(2:end-1) = (F(3:end
3、)-F(1:end-2)/2 , FX(end) = F(end)-F(end-1) ; FX长度与F的长度相同 当F是矩阵时， FX, FY是与F同样大小的矩阵。 FX的每行对应F相应行元素间的梯度 ; FY的每列对应F相应列元素间的梯度 ;,4.1.1 近似数值极限及导数,2019/4/3,第5页,数值极限和导数的应用应十分谨慎,x=eps; L1=(1-cos(2*x)/(x*sin(x), L2=sin(x)/x, L1 = 0 L2 = 1,syms t f1=(1-cos(2*t)/(t*sin(t); f2=sin(t)/t; Ls1=limit(f1,t,0) Ls2=limi
4、t(f2,t,0) Ls1 = 2 Ls2 = 1,x=pi/1000; %可得到正确结果,2019/4/3,第6页,数值极限和导数的应用应十分谨慎,%（1）自变量的增量取得过小（eps数量级） d=pi/100; t=0:d:2*pi; x=sin(t); dt=5*eps; x_eps=sin(t+dt); dxdt_eps=(x_eps-x)/dt; plot(t,x,LineWidth,5) hold on plot(t,dxdt_eps) hold off legend(x(t),dx/dt) xlabel(t),数值导数受计算中有限精度困扰，当增量dt过小时，f(t+dt)与f(t
5、)的数值十分接近，高位有效数字完全相同， df =f(t+dt)-f(t) 造成高位有效数字消失，精度急剧变差。,2019/4/3,第7页,数值极限和导数的应用应十分谨慎,%（2）自变量的增量取得适当 x_d=sin(t+d); dxdt_d=(x_d-x)/d; plot(t,x,LineWidth,5) hold on plot(t,dxdt_d) hold off legend(x(t),dx/dt) xlabel(t),d=pi/100;,2019/4/3,第8页,d=pi/100; t=0:d:2*pi; x=sin(t); dxdt_diff=diff(x)/d; dxdt_gra
6、d=gradient(x)/d;,【例4.1-3】已知 采用diff和gradient计算该函数在区间 中的近似导函数。,subplot(1,2,1); plot(t,x,b);hold on plot(t,dxdt_grad,m,LineWidth,8) plot(t(1:end-1),dxdt_diff,.k,MarkerSize,8) axis(0,2*pi,-1.1,1.1); title(0, 2pi) legend(x(t),dxdt_grad,dxdt_diff,Location,North) xlabel(t), box off;hold off,subplot(1,2,2)
7、kk=(length(t)-10):length(t); hold on; plot(t(kk),dxdt_grad(kk),om,MarkerSize,8) plot(t(kk-1),dxdt_diff(kk-1),.k,MarkerSize,8) title(end-10, end) legend(dxdt_grad,dxdt_diff,Location,SouthEast) xlabel(t),box off; hold off,宏观上, diff和gradient结果大致相同,细节上, diff和gradient数值有差异， diff没有给出最后一点导数,2019/4/3,第9页,Sx
8、=sum(X) % sum按列向求和得(1n)数组 Scs=cumsum(X) %沿X列向求累计和, 仍是数组, 第(i, k)个元素是X数组第k列前i个元素的和。最后一行等于Sx,4.1.2 数值求和与近似数值积分,St=trapz(t,X) 或 St=dt*trapz(X) %梯形法求積分 Sct=cumtrapz(t,X) 或 Sct=dt*cumtrapz(X) %梯形法沿列向求X关于x的累计积分,最后一个值等于St,S=dt*sum(X), S=sum(t,X) %近似矩形法求积分,Scs=cumsum(t,X) =dt*cumsum(X),2019/4/3,第10页,clear;
9、d=pi/8; t=0:d:pi/2; y=0.2+sin(t); s=sum(y); s_sa=d*s; % s_sa=sum(t, y), 近似矩形法积分 s_ta=trapz(t,y); %梯形法积分,disp(sum求得积分,blanks(3),trapz求得积分) disp(s_sa, s_ta) t2=t,t(end)+d; y2=y,nan; stairs(t2,y2,:k);hold on,plot(t,y,r,LineWidth,3) h=stem(t,y,LineWidth,2); set(h(1),MarkerSize,10) axis(0,pi/2+d,0,1.5) h
10、old off; shg sum求得积分 trapz求得积分 1.5762 1.3013,2019/4/3,第11页,4.1.3 计算精度可控的数值积分,一重积分(quadrature精度可控)： S1=quad(fun,a,b,tol) %自适应Simpson法 S2=quadl(fun,a,b,tol) %自适应罗巴托 Lobatlo法 二重积分(精度可控) ： S3=dblquad(fun,xmin,xmax,ymin,ymax,tol) 三重积分(精度可控) ： S4=triplequad(fun,xmin,xmax,ymin,ymax, zmix,zmax,tol) fun:可为字符
11、串,内联对象,匿名函数,函数句柄,注意乘除法运算一定加点.用数组运算 tol: 默认积分的绝对精度为10-6,2019/4/3,第12页,（1） syms x Isym=vpa(int(exp(-x2),x,0,1) % x.2数组乘方亦可 Isym = 0.74682413281242702539946743613185,例4.1-5：求 .,% （2） 梯形法积分 format long d=0.001;x=0:d:1; Itrapz=d*trapz(exp(-x.*x) % x.*x必须为数组乘 Itrapz = 0.746824071499185,（3） fx=exp(-x.2); %
12、一定用数组乘 Ic=quad(fx,0,1,1e-8) Ic = 0.746824132854452 %实际精度控制到1e-10,2019/4/3,第13页,（1）符号计算法 syms x y s=vpa(int(int(xy,x,0,1),y,1,2) Warning: Explicit integral could not be found. s = 0.40546510810816438197801311546435,例4.1-6：求 .,（2）数值积分法 format long s_n=dblquad(x,y)x.y,0,1,1,2) %匿名函数表示被积函数 s_n = 0.40546
13、6267243508,s_n=dblquad(x.y,0,1,1,2) %字符串表示被积函数 s_n=dblquad(inline(x.y),0,1,1,2) %内联函数表示被积函数,2019/4/3,第14页,4.1.4 函数极值的数值求解,x,fval,exitflag,output=fminbnd(fun,x1,x2,options) %一元函数在区间bound(x1,x2)中极小值 x,fval: 极值点坐标和对应目标函数极值 x,fval,exitflag,output=fminsearch(fun,x0,options) %单纯形法求搜索起点x0附近多元函数极值点 % x每列代表一
14、个候选极值点，各列按目标函数极小值递增顺序 %x(:,1)对应的目标函数极小值点由fval给出 fun: 字符串,内联函数,匿名函数,函数句柄,注意乘除法运算一定加点.用数组运算 options: 配置优化参数,可略 exitflag: 给出大于0的数,则成功搜索到极值点 output: 给出具体的优化算法和迭代次数,2019/4/3,第15页,【例4.1-7 】已知 ，在-10x10区间，求函数的最小值。,（1）用“导数为零法”求极值点 syms x y=sin(x)2*exp(-0.1*x)-0.5*sin(x)*(x+0.1); yd=diff(y,x); %求导函数 xs0=solve
15、(yd,x) %求导函数为零的根，即极值点 yd_xs0=vpa(subs(yd,x,xs0) %验证极值点处导函数为零 y_xs0=vpa(subs(y,x,xs0) %求极值点处极值 xs0 = matrix(0.050838341410271656880659496266968) yd_xs0 = 2.2958874039497802890014385492622*10(-41) y_xs0 = -0.001263317776974196724544154344118 无法判断是否最小值,2019/4/3,第16页,【例4.1-7 】已知 ，在-10x10区间，求函数的最小值。,（2）采
16、用优化算法求极小值 x1=-10;x2=10; yx=(x)(sin(x)2*exp(-0.1*x)-0.5*sin(x)*(x+0.1); xn0,fval,exitflag,output=fminbnd(yx,x1,x2) xn0 = 2.514797840754235 fval = %比“导数为零法”求得的极值更小, 更可能是最小值 -0.499312445280039 exitflag = 1 output = iterations: 13 funcCount: 14 algorithm: golden section search, parabolic interpolation m
17、essage: 1x112 char,2019/4/3,第17页,【例4.1-7 】已知 ，在-10x10区间，求函数的最小值。,（4）据图形观察，重设fminbnd的搜索区间 x11=6;x2=10; yx=(x)(sin(x)2*exp(-0.1*x)-0.5*sin(x)*(x+0.1); xn00,fval,exitflag,output=fminbnd(yx,x11,x2) xn00 = 8.023562824723015 fval = -3.568014059128578 %最小值 exitflag = 1 output = iterations: 9 funcCount: 10
18、algorithm: golden section search, parabolic interpolation message: 1x112 char,y=sin(x)2*exp(-0.1*x)-0.5*sin(x)*(x+0.1);,（3）绘图观察最小值 xx=-10:pi/200:10; yxx=subs(y,x,xx); plot(xx,yxx) xlabel(x),grid on,2019/4/3,第18页,例4.1-8: f(x,y)=100(y-x2)2+(1-x)2在区间-5,5的极小值,ff=(x)100*(x(2)-x(1).2).2+(1-x(1).2; %匿名函数 x
19、0=-5,-2,2,5;-5,-2,2,5; %4个搜索起点 sx,sfval,sexit,soutput=fminsearch(ff,x0),其理论极小值点为x=1, y=1,%sx给出一组使优化函数值非减的局部极小点 sx = 0.99998 -0.68971 0.41507 8.0886 0.99997 -1.9168 4.9643 7.8004 sfval = 2.4112e-010,2.4112e-010 5.7525e+002 2.2967e+003 3.3211e+005,format short e %取5位科学计数法 disp(ff(sx(:,1),ff(sx(:,2),ff
20、(sx(:,3),ff(sx(:,4),x0 =-5 -2 2 5 -5 -2 2 5,2019/4/3,第19页,4.1.5 常微分方程Ordinary Differential Equation的数值解,t,Y=ode45(odefun,tspan,y0) % 4阶龙格库塔 数值法 odefun: 待求解一阶微分方程组的函数文件句柄 tspan: 自变量微分二元区间t0, tf y0: 一阶微分方程组的(n*1)初值列向量,matlab为解常微分方程初值问题提供了一组配套齐全,结构严整的指令, 包括: ode45, ode23, ode113, ode23t, ode15s, ode23s
21、, ode23tb.在此只介绍最常用的ode45的基本使用方法.,ode45使用方法:,t-二元区间的点系列,Y-原函数在微分区间点系列上的函数值,2019/4/3,第20页,例4.1-9求解：,%解算微分方程 tspan=0,30; y0=1;0; tt,yy=ode45(DyDt,tspan,y0); figure(1) plot(tt,yy(:,1) xlabel(t),title(x(t),%画相平面图(函数和其导数勾画的曲线称为相轨迹) figure(2) plot(yy(:,1),yy(:,2) xlabel(位移),ylabel(速度),2019/4/3,第21页,4.2矩阵和代
22、数方程,4.2.1 矩阵运算和特征参数 4.2.2 矩阵的变换和特征值分解 4.2.3 线性方程的解 4.2.4 一般代数方程的解,2019/4/3,第22页,4.2.1 矩阵运算和特征参数,矩阵与标量之间的四则运算与数组运算相同 矩阵和矩阵之间的四则运算 矩阵和矩阵之间的加减运算与数组运算相同 设 A 是一个 mn 矩阵，B 是一个 pq 矩阵，当 np 时，两个矩阵可以相乘，乘积为 mq 矩阵。矩阵乘法不可逆。在 MATLAB 中，矩阵乘法由“*”实现。 矩阵除法在实际中主要用于求解线性方程组 矩阵转置： 符号“ ”实现矩阵的转置操作。对于实数矩阵， “ ”表示矩阵转置，对于复数矩阵，“
23、”实现共轭转置。对于复数矩阵，如果想要实现非共轭转置，可以使用符号“ . ”。,2019/4/3,第23页,format rat %有理格式显示 A=magic(2) + j*pascal(2) A = 1 + 1i 3 + 1i 4 + 1i 2 + 2i A1=A A1=1 - 1i 4 - 1i %共轭转置 3 - 1i 2 - 2i A2=A. A2=1 + 1i 4 + 1i %非共轭转置，数组运算操作 3 + 1i 2 + 2i,例4.2-2 矩阵和数组转置操作的差别,2019/4/3,第24页,计算矩阵标量特征参数-秩,迹,行列式的指令,rank(A) %求秩（Rank） det
24、(A) %求行列式（Determinant） trace(A) %求迹（Trace），即矩阵主对角元素的和,2019/4/3,第25页,A=reshape(1:9,3,3); r=rank(A) %求秩 d3=det(A) %非满秩矩阵的行列式一定为零 d2=det(A(1:2,1:2) %求子式的行列式 t=trace(A),【例4.2-3】矩阵标量特征参数计算示例。,A = 1 4 7 2 5 8 3 6 9,r = 2,d3 =0,d2 = -3,t = 15,2019/4/3,第26页,【例4.2-4】矩阵标量特征参数的性质。,format short ; rand(twister,0
25、) A=rand(3,3); B=rand(3,3); C=rand(3,4); D=rand(4,3); tAB=trace(A*B) %任何符合矩阵乘法规则的两个矩阵乘积 tBA=trace(B*A) %的“迹”不变。同阶乘积“迹”不变 tCD=trace(C*D) %两个“內维”相等矩阵的乘积“迹”不变 tDC=trace(D*C),tAB = 2.6030 tBA = 2.6030 tCD = 4.1191 tDC = 4.1191,dCD=det(C*D) dDC=det(D*C) dCD = 0.0424 dDC = -2.6800e-018,d_A_B = 0.0094 dAB
26、= 0.0094 dBA = 0.0094,非同阶矩阵乘积行列式不等于各矩阵行列式之乘积,2019/4/3,第27页,dx=diff(X) %计算向量X的前向差分 FX=gradient(F) %求一元(函数)梯度 FX, FY =gradient(F) %求二元（函数）梯度,4.1.1 近似数值极限及导数,S1=sum(X,1) % sum按列向求和得(1n)数组, =sum(X) (X多行) S2=sum(X,2) % sum按行向求和得(n1)数组 Scs=cumsum(X)%沿X列向求累计和, 仍是数组, 最后一行等于Sx,4.1.2 数值求和与近似数值积分,St=trapz(t,X)
27、 或 St=dt*trapz(X) %梯形法求積分 Sct=cumtrapz(t,X) 或 Sct=dt*cumtrapz(X) %梯形法沿列向求X关于x的累计积分,最后一个值等于St,Review,2019/4/3,第28页,4.1.3 计算精度可控的数值积分,一重积分(quadrature精度可控)： S1=quad(fun,a,b,tol) %自适应Simpson法 S2=quadl(fun,a,b,tol) %自适应罗巴托 Lobatlo法 二重积分(精度可控) ： S3=dblquad(fun,xmin,xmax,ymin,ymax,tol) 三重积分(精度可控) ： S4=trip
28、lequad(fun,xmin,xmax,ymin,ymax, zmix,zmax,tol) fun:可为字符串,内联对象,匿名函数,函数句柄,注意乘除法运算一定加点.用数组运算 tol: 默认积分的绝对精度为10-6,Review,2019/4/3,第29页,4.1.4 函数极值的数值求解,x,fval,exitflag,output=fminbnd(fun,x1,x2,options) %一元函数在区间bound(x1,x2)中极小值 x,fval: 极值点坐标和对应目标函数极值 x,fval,exitflag,output=fminsearch(fun,x0,options) %单纯形法
29、求搜索起点x0附近多元函数极值点 % x每列代表一个候选极值点，各列按目标函数极小值递增顺序 %x(:,1)对应的目标函数极小值点由fval给出 fun: 字符串,内联函数,匿名函数,函数句柄, 注意乘除法运算一定加点.用数组运算 options: 配置优化参数,可略 exitflag: 给出大于0的数,则成功搜索到极值点 output: 给出具体的优化算法和迭代次数,Review,2019/4/3,第30页,rank(A) %求秩（Rank） det(A) %求行列式（Determinant） trace(A) %求迹(Trace)，即矩阵主对角元素的和,Review,t,Y=ode45(o
30、defun,tspan,y0) % 4阶龙格库塔 数值法 odefun: 待求解一阶微分方程组的函数文件句柄 tspan: 自变量微分二元区间t0, tf y0: 一阶微分方程组的(n*1)初值列向量,ode45使用方法:,t-二元区间的点系列,Y-原函数在微分区间点系列上的函数值,4.1.5 常微分方程Ordinary Differential Equation的数值解,4.2.1 矩阵运算和特征参数,2019/4/3,第31页,4.2.2 矩阵的变换和特征值分解,R, ci=rref(A) %借助初等变换将A缩减行变成行阶梯矩阵R。,B=orth(A)可得到矩阵A的正交基，B的列与A的列可
31、张成相同的空间，而且B的列是正交的，因此B*B=eye(rank(A)，B的列数正好是A的秩。,X=null(A) %A矩阵零空间的全部正交基，满足AX=0, X*X=I。,B=orth(A),V, D=eig(A) %A矩阵的特征值、特征向量分解，使AV=VD。,ci 是行数组, 其元素表示A中线性独立列的序号。length(ci)=rank(A),2019/4/3,第32页,【例4.2-5】行阶梯阵简化指令rref计算结果的含义。,A=magic(4) R,ci=rref(A) %行阶梯分解 A = 16 2 3 13 5 11 10 8 9 7 6 12 4 14 15 1 R = 1
32、0 0 1 0 1 0 3 0 0 1 -3 0 0 0 0 ci = 1 2 3,r_A=length(ci) % =rank(A) r_A = 3,aa=A(:,1:3)*R(1:3,4) %A前三列线形组合 err=norm(A(:,4)-aa) aa = 13 8 12 1 err = 0,2019/4/3,第33页,A=reshape(1:15,5,3) X=null(A) S=A*X,S =1.0e-014 * 0 -0.1776 -0.2665 -0.3553 -0.5329,A =1 6 11 2 7 12 3 8 13 4 9 14 5 10 15,n= 3 l= 1 Ran
33、k_A=2,X =0.4082 -0.8165 0.4082,ans = 1,例4.2-6 矩阵零空间及含义。设 是矩阵 的零空间, 即,n=size(A,2) l=size(X,2) Rank_A=rank(A) n-l=rank(A),2019/4/3,第34页,【例4.2-7】简单实阵的特征分解。eig，cdf2rdf Complex Diagonal Form to Real-block Diagonal Form,%（1） A=1,-3;2,2/3 V,D=eig(A) A =1.0000 -3.0000 2.0000 0.6667 V = 0.7746 0.7746 0.0430
34、- 0.6310i 0.0430 + 0.6310i D =0.8333 + 2.4438i 0 0 0.8333 - 2.4438i %（2） VR,DR=cdf2rdf(V,D) VR =0.7746 0 0.0430 -0.6310 DR =0.8333 2.4438 -2.4438 0.8333,%（3） A1=V*D/V %因计算误差有很小虚部 A1_1=real(A1) %去除虚部后等于A A2=VR*DR/VR err1=norm(A-A1,fro) err2=norm(A-A2,fro) A1 =1.0000 + 0.0000i -3.0000 2.0000 - 0.0000i
35、 0.6667 A1_1 =1.0000 -3.0000 2.0000 0.6667 A2 =1.0000 -3.0000 2.0000 0.6667 err1 =6.7532e-016 err2 =4.4409e-016,Frobenius矩阵范数,矩阵Frobenius范数,类似向量2范数，不同于矩阵2范数,2019/4/3,第35页,4.2.2 线性方程的解,对于方程组Amnx=b (m个方程，n个未知数)， 当向量b在矩阵A列向量所张空间中, rank(A,b) =rank(A)= r a) 若n =r，则解唯一。 b) 若n r，则解不唯一。 当向量b不在矩阵A列向量所张空间中, 则
36、无准确解但存在最小二乘解。,1. 线性方程解的一般讨论,matlab定义的左除运算可以很方便地解上述方程组： x=Ab (x=A-1b只适用于A非奇异时),2.除法运算解方程解, 当m=n时, “恰定”方程 当mn时, “超定”方程 当mn时, “欠定”方程,2019/4/3,第36页,% （1）创建待解方程的A和b A=reshape(1:12,4,3); b=(13:16); % （2）若rank(A,b) =rank(A),则b在A的值空间中 ra=rank(A), rab=rank(A,b) ra = 2 rab = 2 % nr, 解不唯一 % （3）求特解和通解，并对由之构成的全解
37、进行验算 xs=Ab; xg=null(A); % xg是齐次方程Ax=0的解 c=rand(1); ba=A*(xs+c*xg) % ba是A乘 “一个随机的全解” Warning: Rank deficient, rank = 2, tol = 1.8757e-014. ba =13.0000 14.0000 15.0000 16.0000,例4.2-8: 求方程 的解,norm(ba-b) ans = 1.1374e-014,2019/4/3,第37页,例4.2-9: “逆阵法”和”左除法”解恰定方程的性能对比。,测试阵,产生指定异常值和特殊带宽的随机阵,3. 矩阵的逆,A_1=inv(
38、A) %求非奇异阵A的逆,randn(state,0); A=gallery(randsvd,300,2e13); %产生条件数 为2e13的300阶随机矩阵; x=ones(300,1); %定义真解 b=A*x; cond(A) %验算矩阵条件数, 结果a1.9978e+013, 值越大, 阵越病态, %求逆法 tic; xi=inv(A)*b; ti=toc eri=norm(x-xi) rei=norm(A*xi-b)/norm(b),%左除法 tic; xd=Ab;td=toc erd=norm(x-xd) red=norm(A*xd-b)/norm(b),ti =0.0185 er
39、i =0.0883 rei =0.0051,td =0.0066 erd =0.0298 red =8.7810e-015,矩阵计算对于误差越敏感, 数值稳定性差,2019/4/3,第38页,求解任意函数f(x)=0(可能无解,单解,多解)的步骤: 作图获取初步近似解 观察f(x)与横轴的交点坐标，用zoom放大，用ginput得较精确些的交点坐标值。,4.3 一般代数方程的解,用”泛函”指令求精确解 x,favl = fzero(fun,x0) %一元函数求零点 x,fval = fsolve(fun,x0) %解非线性方程组 fun: 内联对象,匿名函数,函数句柄,字符串;被解函数自变量一
40、般用x, 注意乘除法运算一定加点.用数组运算 x0: 零点初始猜测值. x0为标量时取与之最靠近的零点; x0取a,b时,在此区间内寻找一个零点. x: 所求零点的自变量值 fval: 函数值,2019/4/3,第39页,例4.2-10: 求f(t)=(sin2t)e-0.1t0.5|t|的零点,y_C=inline(sin(t).2.*exp(-0.1*t)-0.5*abs(t),t); t=-10:0.01:10; Y=y_C(t); clf, plot(t,Y,r); hold on plot(t,zeros(size(t),k); xlabel(t);ylabel(y(t),t4 =
41、0.5993 y4 = 0,tt =-2.0039 -0.5184 -0.0042 0.6052 1.6717,t3 = -0.5198 y3 = 5.5511e-017,zoom on tt,yy=ginput(5),zoom off t4,y4=fzero(y_C,tt(4) t3,y3=fzero(y_C,tt(3) %t=0处没穿越横轴,2019/4/3,第40页,4.3 概率分布和统计分析,4.3.1 概率函数、分布函数、逆分布函数和随机数的发生 4.3.2 随机数发生器和 统计分析指令,2019/4/3,第41页,4.3 概率分布和统计分析,1. 二项分布（Binomial dis
42、tribution ）,4.3.1 概率函数、分布函数、逆分布函数和随机数的发生,pk=binopdf(k,N,p) 事件A发生k次的概率,Fk=binocdf(k,N,p) 事件A发生次数不大于k的概率,R=binornd(N,p) 产生符合二项分布B(N,p)的(mn) 随机数组（元素值为事件A可能发生的次数）。,Binomial Cumulative Distribution Function,Binomial Probability Density Function,2019/4/3,第42页,N=100;p=0.5; k=0:N; pdf=binopdf(k,N,p); cdf=bi
43、nocdf(k,N,p); h=plotyy(k,pdf,k,cdf); set(get(h(1),Children),Color,b,Marker,.,MarkerSize,13) set(get(h(1),Ylabel),String,pdf) set(h(2),Ycolor,1,0,0) set(get(h(2),Children),Color,r,Marker,+,MarkerSize,4) set(get(h(2),Ylabel),String,cdf) xlabel(k) grid on,例4.3-1. 画出N=100，p=0.5情况下的二项分布概率特性曲线。,4.3.1 概率函数
44、、分布函数、逆分布函数和随机数的发生,2019/4/3,第43页,px=normpdf(x,Mu,Sigma),Fx=normcdf(x,Mu,Sigma),R=normrnd(Mu,Sigma,m,n),2. 正态分布（Normal distribution ）,服从 分布的随机变量取值x的概率密度。,服从 分布的随机变量取值不大于x的概率分布。,产生元素服从 分布的(mn) 随机数组。,x可取任何实数，Mu是正态分布的数学期望, Sigma 是正态分布的均方差。,2019/4/3,第44页,例4.3-2. 正态分布标准差的几何表示。,mu=3; sigma=0.5; x=mu+sigma*
45、-3:-1,1:3 yf=normcdf(x,mu,sigma) x = 1.5 2.0 2.5 3.5 4.0 4.5 yf = 0.0013 0.0228 0.1587 0.8413 0.9772 0.9987 P=yf(4)-yf(3),yf(5)-yf(2),yf(6)-yf(1) P = 0.6827 0.9545 0.9973 xd=1:0.1:5; yd=normpdf(xd,mu,sigma);,2019/4/3,第45页,例4.3-2. 正态分布标准差的几何表示。,clf for k=1:3 xx=x(4-k):sigma/10:x(3+k); yy=normpdf(xx,m
46、u,sigma); subplot(3,1,k),plot(xd,yd,b); hold on fill(x(4-k),xx,x(3+k),0,yy,0,g); hold off if k2 text(3.8,0.6,mu-sigma,mu+sigma) else kk=int2str(k); text(3.8,0.6,mu-,kk,sigma,mu+,kk,sigma) end text(2.8,0.3,num2str(P(k);shg end; xlabel(x); shg,x =1.5 2 2.5 3.5 4 4.5,k=1, x(4-k) =2.5, x(3+k) =3.5,2019/4/3,第46页,3.各种概率分布(distiburion)的交互式观察界面,disttool,2019/4/3,第47页,4.3.2 随机数发生器和统计分析指令,2019/4/3,第48页,例4.3-4. 随机数据的统计量。,randn(state,0) % 也可用randn(twister,0) A=randn(1000,4); AMAX=max(A) AMIN=min