matlab数学建模数据分析ppt课件.ppt
第四单元 数据分析,Matlab统计工具箱,一:统计工具箱简介二:概率分布三:参数估计四:描述性统计五:假设检验六:统计绘图,一.matlab统计工具箱(statistics toolbox)简介,统计学是处理数据的艺术和科学,通过收集,分析,解释和表达数据来探索事物中蕴含的规律.随着科技水平的迅猛发展,知识经济的时代来临,海量的数据需要人们处理.matlab统计工具箱为人们提供了一个强有力的统计分析工具. 统计工具箱基于matlab数值计算环境,支持范围广泛的统计计算任务.它包括200多个处理函数(m文件)主要应用于以下几方面:,1.1 统计工具箱的几大功能,*概率分布 *参数估计 *描述性统计 *假设检验 *统计绘图,统计工具箱提供了20种概率分布类型,其中包括离散型分布:(如binomial二项分布, 即n次贝努里试验中出现k次成功的概率.poisson 分布, 和 分布等).,1.1.1概率分布-离散型,1.1.2 概率分布连续型,连续型分布如正态分布F(x)=beta分布,uniform平均分布等.每种分布提供5类函数: 1 概率密度 2 (累积)分布函数 3 逆累积分布函数 4 随机数产生器 5 均值和方差函数.,1.1.3另外4大功能,*参数估计-依据原始数据计算参数估计值置信区域.*描述性统计-方差,期望等数字特征.*假设检验-提供最通用的假设检验函数t-检验,z-检验.*统计绘图- box图函数,正态概率图函数等.注意:统计工具箱中的说有函数都可用 type function_name语句查看其代码,也可进行修改,从而变为己用,加入到工具箱中.,二 概率分布,随机变量的统计行为取决于其概率分布,而分布函数常用连续和离散型分布。统计工具箱提供20种分布。每种分布有五类函数。1: 概率密度(pdf) ; 2: 累积分布函数(cdf); 3:逆累积分布函数(icdf);4: 随机数产生器 5: 均值和方差函数;一:离散型概率密度函数:为观察到的特定值的概率。连续型概率密度函数定义为:如存在非负函数p(x) 0,使对任意ba,X 在(a,b)上取值概率为paXb= ;则称p(x)为随机变量X的概率密度函数。二:累积分布 (cdf):它取决于pdf. 表达式为F(x)= . 逆累积分布(icdf):实际上是cdf的逆,它返回给定显著概率条件下假设检验的 临界值。,2.1,三:随机数产生器 所有随机数产生方法都派生于均匀分布随机数。产生方法有:直接法、反演法、拒绝法。四:均值和方差 均值和方差是分布函数的简单函数。在Matlab里用“stat”结尾的函数可计算得到给定参数的分布的均值和方差。以下以正态分布为例说明在Matlab里的实现。一:概率密度函数 X=-3:0.5:3;f=normpdf(x,0,1);(其中normpdf为正态分布的Matlab分布实现函数,可由以下介绍的函数代替。),功能:可选分布的概率密度函数。格式:Y=pdf(name,X,A1,A2,A3)说明:name为特定分布的名称,如Normal,Gamma等。X为分布函数的自变量X的取值矩阵,而A1,A2,A3分别为相应的分布参数值。Y给出结果,为概率密度值矩阵。举例:p=pdf(Normal,-2:2,0,1) 给出标准正态分布在-2到2的分布函数值。 而p=pdf(Poisson,0:4,1:5)给出Poisson分布函数。,2.2,累积分布函数与逆累积分布函数同样地,累积分布和逆累积分布对每个分布都有特定地Matlab实现函数,这里只介绍通用的cdf,icdf.- cdf, icdf功能:计算可选分布的累积分布函数和逆累积分布函数。格式:P=cdf(name,X,A1,A2,A3) X=icdf(name,X,A1,A2,A3)说明:cdf和icdf中的参数使用和pdf中的相同。只是计算结果不同。举例:p=cdf(Normal,0:5,1:6) X=icdf(Normal,0.1:0.2:0.9,0,1),2.3,随机数产生器在Matlab里和pdf,cdf与icdf一样,随机数的产生也有通用函数random.- random功能:产生可选分布的随机数。格式:y=random(name,A1,A2,A3,m,n)说明:random函数产生统计工具箱中任一分布的随机数。name为相应分布的名称。A1,A2,A3为分布参数,意义同pdf参数。m,n确定了结果y的数量,如果分布参数A1,A2,A3为矢量,则m,n是可选的,但应注意,它们给出的长度或矩阵行列数必须与分布参数的长度相匹配。举例:rn=random(Normal,0,1,2,4),2.4,均值和方差和以上其他函数不同的是均值和方差的运算没有通用的函数,只能用各个分布的函数计算。对应于正态分布的计算函数为normstat();它返回两个参数的向量,分别为均值和方差。举例:m,n=normstat(mu,sigma),2.5,三.参数估计,参数估计: 某分布的数学形式已知,应用子样信息来估计其有限个参数的值本节主要介绍 3.1 最大似然估计(Maximum likelihood estimation) 3.2 对数似然函数,3.1最大似然估计,基本思想: 已知一组观测值,给定这组值出自的某类分布中,求得最有可能出现这组值的一个分布.调用方法: phat,pci=mlsdist,data,alpha phat为参数估计结果,pci为置信区间计算结果dist为用户给定的分布名称,data为数据列表,(1-alpha)置信区域.,3.1.1 最大似然估计(mls)举例,例: rv=binornd(20,0.75) rv= 17 p,pci=mle(binomial,rv,0.05,20) p= 0.8000 pci= 0.5634 0.9427,3.2 对数似然函数,统计工具箱提供了分布,分布,正态分布和威布尔分布的负对数似然函数值的求取函数.正态分布的负对数似然函数调用方法 L=normlike(params,data) Params为正态分布参数:params(1)为,params(2)为,3.2.1其他负对数似然函数,分布的负对数似然函数 logL=betalike(params,data)分布的负对数似然函数 logL=gamlike(params,data)威布尔分布的负对数似然函数 logL=weiblike(params,data) 参数设置与正态分布的负对数似然函数类似,不加冗述.,四 描述性统计,概述:人们希望用少数样本来体现样本总体的规律。描述性统计就是收集、整理、加工和分析统计数据,使之系统化、条理化,以显示出数据资料的趋势、特征和数量关系。根据统计量特征性质的不同,工具箱提供了位置度量、散布度量、自助法以及在缺失数据情况下处理方法等方面的描述性统计工具函数。,4.1中心趋势(位置)度量,数据样本中心度量的目的在于对数据样本的数据分布线上分布的中心予以定位,即中心位置的度量。均值是对位置的简单和通常的估计量。但野值的存在往往影响位置的确定。而中位数和修正的均值则受野值的干扰很小。中位数是样本的50%分位点。而修正的均值所蕴涵的思想则是剔除样本中最高值和最低值来确定样本的中心位置。几何均值和调和均值对野值都较敏感。当样本服从对数正态分布或偏斜程度很大时,它们也都是有效的方法。以下介绍位置度量有关函数。,4.2.1: 几何平均数(geomean),功能:样本的几何均值。格式:m=geomean(X)说明:几何均值的定义为 m= (1.4.1) geomean 函数计算样本的几何均值 。X若为矢量,它返回X中元素的几何均值;X若为矩阵,它的结果为一个行矢量,每个元素为X对应列元素的几何均值。举例:x=exprnd(1,10,6); geometric=geomean(X); average=mean(X);,4.2.2: (调和均值)harmmean,功能:样本数据的调和均值。格式:m=harmmean(X)说明:调和均值定义为举例:样本均值大于或等于调和均值。 X=exprnd(1,10,6); harmonic=harmmean(X) average=mean(X),4.2.3(平均值)mean,功能:样本数据的平均值。说明:平均值定义为举例:x=normrnd(0,1,100,5); xbar=mean(X),4.2.4:median,功能:样本数据的中值。说明:中值即数据样本的50%中位数。中位数对野值出现的影响较小。举例:xodd=1:5; modd=median(xodd) meven=median(xeven),4.2.5:trimmean,功能:剔除极端数据的样本均值。格式:m=trimmean(X,percent)说明:函数计算剔除观测量中最高百分比和最低百分比数据后的均值。 函数中percent代表百分比。举例:X=normrnd(0,1,100,100); m=mean(X) trim=trimmean(X,10) sm=std(m) strim=std(trim) efficiency=(sm/strim).2,4.3散布度量,散布度量可以理解为样本中的数据偏离其数值中心的程度,也称离差。极差,定义为样本最大观测值与最小观测值之差。标准差和方差为常用的散布度量,对正态分布的样本描述是最优的。但抗野值干扰能力较小。平均绝对值偏差对野值也敏感。四分位数间距为随机变量的上四分位数 和下四分位之差。,在Matlab里,有关散布度量计算的函数为:1:计算样本的内四分位数间距的 iqr(X).2:求样本数据的平均绝对偏差的 mad(X).3:计算样本极差的 range(X).4: 计算样本方差的 var(X,w).5: 求样本的标准差的 std(X).6: 求协方差矩阵的cov(X).这些函数的详细说明可以参见Matlab的帮助文档。,4.4 Matlab里有关散布度量计算的函数,4.5处理缺失数据的函数,在对大量的数据样本进行处理分析时,常会遇到一些数据无法找到或不能确定的情况。这时可用NaN标注这个数据。而工具箱中有一些函数自动处理它们。如 :忽视NaN, 求其他数据的最大值的nanmax.格式:m=nanmax(X)举例:m=magic(3); m(1 6 8)=NaN NaN NaN nmax,maxidx=nanmax(m),4.6中心矩,中心矩是关于数学期望的矩。对于任意的r 0,称 为随机变量X的r阶中心矩。一阶中心矩为0,二阶中心矩为方差: 函数moment计算任意阶中心矩。 格式:m=moment(X,order) 说明:order确定阶。,4.7相关系数,相关系数是两个随机变量间线性相依程度的度量。可用函数corrcoef计算它。格式:R=corrcoef(X)说明:输入矩阵X的行元素为观测值,列元素为变量,R=corrcoef(X)返回相关系数矩阵R.,五.假设检验,假设检验 是统计的基本问题.旨在应用得到的少量信息,判断整体是否满足给定条件或达到给定的标准. 回顾一下我们以前在统计学中所学的假设检验.其步骤为:,5.1 假设检验步骤,1.设: 零假设.(成立则h=0,否则h=1).2.取得一组观测值(子样).3.给定显著型水平(一般取0.05). 4.应用子样的某些统计量特征.5.在 成立前提下,若出现已知观测值的概率小于5%,则拒绝,否则认为观测值与假设无显著差别.,5.2 ranksum函数,调用方法: p,h=ranksum(x,y,alpha) p返回x,y的母体一致的显著性水平,h为假设检验的返回值.x,y为两组观测值,alpha为显著性水平. 请参考下面例子,5.2.1 Ranksum函数举例,例:检验两组服从poisson分布的随机数样本的均值是否相同. x=poissrnd(5,10,1); y=poissrnd(2,10,1); p,h=ranksum(x,y,0.05) p= 0.0028 h= 1,5.3 signrank函数,调用方法: p,h=signrank(x,y,alpha) 参数与ranksum函数类似. 例:检验两个正态分布的样本子样均值是否相等. x=normrnd(0,1,20,1); y=normrnd(0,2,20,1); p,h=signrank(x,y,0.05) p= 0.2568 h= 0,5.4 ttest-t检验,调用方法: h,sig,ci=ttest(x,m,alpha) h为假设检验的返回值.sig与T统计量有关,T= . ci为均值的(1-alpha)置信区域.m为假设的样本均值.,5.4.1 ttest函数举例,例:给出理论均值为0、标准差为1的100个正态随机数样本。当然,观测样本的均值和标准差与理论值不同的,但假设检验的结果却还原其本质规律。 x=normrnd(0,1,1,100); h,sig,ci=ttest(x,0); h= 0 sig= 0.4474 ci= -0.1165 0.2620 结果h=0,意味着我们不能拒绝零假设。,5.5 ztest函数,已知方差的单样本均值的检验假设.调用方法: h,sig,ci=ztest(x,m,sigma,alpha,tail) ztest(x,m,sigma)是在0.05显著性水平下检验正态分布的样本是否具有均值m和标准差sigma. h=ztest(x,m,sigma,alpha)则可由您确定显著性水平alpha值,并返回检验结果h。 Sig、ci与ttest函数中相应的意义相同。,5.5.1函数ztest举例,例: x=normrnd(0,1,100,1); m=mean(x); m=0.0727 h,sig,ci=ztest(x,0,1); h= 0 sig= 0.4669 ci= -0.1232 0.2687,六 统计绘图,概述统计工具箱在Matlab丰富的绘图功能上又添加了图形表现函数,box图用于展现样本及其统计量的内在规律,也用于通过图形来比较多个样本的均值。正态概率图是确定样本是否为正态分布的图形。分位数-分位数图用于比较两个样本的分布。,6.1 Box图,-boxplot功能:数据样本的box图。格式:boxplot(X) boxplot(X,notch,sym,vert,whis) 举例:x1=normrnd(5,1,100,1); x2=normrnd(6,1,100,1); x=x1 x2; boxplot(x,1),6.2误差条图,-errorbar功能:误差条图。格式:errorbar(X,Y,L,U,symbol)举例:lambda=(0.1:0.2:0.5); r=poissrnd(lambda(ones(50,1),:)); p,pci=poissfit(r,0.001); L=p-pci(1,: ) U=pci(2,: )-p errorbar(1:3,p,L,U,+),还有其他函数:1: fsurfht 画交互轮廓图2: gline 绘制交互3:gname 用实例名称或实例号来标记图中的点4: lsline 绘制数据的最小二乘拟合线5: normplot 图形化正态检验的正态概率图6: pareto 帕累托图7: qqplot 两个样本的分位数-分位数图8 :rcoplot 回归残差图9: refcurve 在当前图形中给出多项式拟合曲线,6.3,几个统计绘图的例子,画正态概率图Normplot(x)画数据的正态概率图X=normrnd(0,1,50,1)H=normplot(x);,areto图,Pareto(y,names)defects=pits ;cracks;holes ;dents ;quantify5,3,19,25;quantity=5,3,19,25;,用实例名来标记图中的点,Gname(case)功能:用实例名来标记图中的点Load citiesEduation=rating(:,6);arts=ratings(:,7);Plot(eduation,artsk,+)Gname(names),在大量的应用领域中,人们经常面临用一个解析函数描述数据(通常是测量值)的任务。对这个问题有两种方法。一种是插值法,数据假定是正确的,要求以某种方法描述数据点之间所发生的情况。另一种方法是曲线拟合或回归。人们设法找出某条光滑曲线,它最佳地拟合数据,但不必要经过任何数据点。,本专题的主要目的是:了解插值和拟合的基本内容;掌握用Matlab求解插值与拟合问题的基本命令。,内容提纲,1.拟合问题引例及基本理论2.Matlab求解拟合问题3.应用实例4.插值问题引例及基本理论5.Maltab求解插值问题6.应用实例,拟合问题,拟 合 问 题 引 例 1,温度t(0C) 20.5 32.7 51.0 73.0 95.7电阻R() 765 826 873 942 1032,已知热敏电阻数据:,求600C时的电阻R。,设 R=at+ba,b为待定系数,拟 合 问 题 引 例 2,求血药浓度随时间的变化规律c(t).,作半对数坐标系(semilogy)下的图形,曲 线 拟 合 问 题 的 提 法,已知一组(二维)数据,即平面上 n个点(xi,yi) i=1,n, 寻求一个函数(曲线)y=f(x), 使 f(x) 在某种准则下与所有数据点最为接近,即曲线拟合得最好。,y=f(x),i 为点(xi,yi) 与曲线 y=f(x) 的距离,线性最小二乘拟合 f(x)=a1r1(x)+ +amrm(x)中函数r1(x), rm(x)的选取,1. 通过机理分析建立数学模型来确定 f(x);,2. 将数据 (xi,yi) i=1, n 作图,通过直观判断确定 f(x):,曲线拟合问题最常用的解法线性最小二乘法的基本思路,第一步:先选定一组函数 r1(x), r2(x), rm(x), mn, 令 f(x)=a1r1(x)+a2r2(x)+ +amrm(x) (1)其中 a1,a2, am 为待定系数。,第二步: 确定a1,a2, am 的准则(最小二乘准则):使n个点(xi,yi) 与曲线 y=f(x) 的距离i 的平方和最小 。,记,问题归结为,求 a1,a2, am 使 J(a1,a2, am) 最小。,线性最小二乘法的求解:预备知识,超定方程组:方程个数大于未知量个数的方程组,超定方程一般是不存在解的矛盾方程组。,如果有向量a使得 达到最小,则称a为上述超定方程的最小二乘解。,线性最小二乘法的求解,定理:当RTR可逆时,超定方程组(3)存在最小二乘解, 且即为方程组 RTRa=RTy -正则(正规)方程组的解:a=(RTR)-1RTy,所以,曲线拟合的最小二乘法要解决的问题,实际上就是求以下超定方程组的最小二乘解的问题。,用MATLAB解拟合问题,1、线性最小二乘拟合,2、非线性最小二乘拟合,用MATLAB作线性最小二乘拟合,1. 作多项式f(x)=a1xm+ +amx+am+1拟合,可利用已有命令:,a=polyfit(x,y,m),3.对超定方程组,2.多项式在x处的值y的计算命令:y=polyval(a,x),例 对下面一组数据作二次多项式拟合,1)输入命令:x=0:0.1:1; y=-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.30,11.2; R=(x.2), x, ones(11,1); A=Ry,解法1解超定方程的方法,2)计算结果: = -9.8108, 20.1293, -0.0317,2)计算结果: = -9.8108, 20.1293, -0.0317,解法2用多项式拟合的命令,MATLAB(zxec2),1)输入命令:x=0:0.1:1; y=-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.30,11.2;A=polyfit(x,y,2)z=polyval(A,x);plot(x,y,k+,x,z,r) %作出数据点和拟合曲线的图形,1. lsqcurvefit已知数据点:xdata=(xdata1,xdata2,xdatan) ydata=(ydata1,ydata2,ydatan),用MATLAB作非线性最小二乘拟合,两个求非线性最小二乘拟合的函数:lsqcurvefit、lsqnonlin。相同点和不同点:两个命令都要先建立M-文件fun.m,定义函数f(x),但定义f(x)的方式不同,请参考例题。,lsqcurvefit用以求含参量x(向量)的向量值函数F(x,xdata)=(F(x,xdata1),F(x,xdatan)T中的参变量x(向量),使得,输入格式: (1) x = lsqcurvefit (fun,x0,xdata,ydata); (2) x =lsqcurvefit (fun,x0,xdata,ydata,lb, ub); (3) x =lsqcurvefit (fun,x0,xdata,ydata, lb, ub, options); (4) x, resnorm = lsqcurvefit (fun,x0,xdata,ydata,); (5) x, resnorm, residual = lsqcurvefit (fun,x0,xdata,ydata,);,说明:x = lsqcurvefit (fun,x0,xdata,ydata,options);,lsqnonlin用以求含参量x(向量)的向量值函数 f(x)=(f1(x),f2(x),fn(x)T 中的参量x,使得 最小。 其中 fi(x)=f(x,xdatai,ydatai) =F(x,xdatai)-ydatai,2. lsqnonlin,已知数据点: xdata=(xdata1,xdata2,xdatan) ydata=(ydata1,ydata2,ydatan),输入格式: 1) x=lsqnonlin(fun,x0); 2) x= lsqnonlin (fun,x0,lb,ub); 3) x= lsqnonlin (fun,x0, ,lb,ub,options); 4) x, resnorm= lsqnonlin (fun,x0,); 5) x, resnorm , residual= lsqnonlin (fun,x0,);,说明:x= lsqnonlin (fun,x0,options);,例2 用下面一组数据拟合 中的参数a,b,k,该问题即解最优化问题:,MATLAB(fzxec1),1)编写M-文件 curvefun1.m function f=curvefun1(x,tdata) f=x(1)+x(2)*exp(-0.02*x(3)*tdata) %其中 x(1)=a; x(2)=b;x(3)=k;,2)输入命令tdata=100:100:1000cdata=1e-03*4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59; x0=0.2,0.05,0.05; x=lsqcurvefit (curvefun1,x0,tdata,cdata) f= curvefun1(x,tdata),F(x,tdata)= ,x=(a,b,k),解法1. 用命令lsqcurvefit,3)运算结果:f =0.0043 0.0051 0.0056 0.0059 0.0061 0.0062 0.0062 0.0063 0.0063 0.0063 x =0.0063 -0.0034 0.2542,4)结论:a=0.0063, b=-0.0034, k=0.2542,MATLAB(fzxec2),1)编写M-文件 curvefun2.m function f=curvefun2(x) tdata=100:100:1000; cdata=1e-03*4.54,4.99,5.35,5.65,5.90, 6.10,6.26,6.39,6.50,6.59; f=x(1)+x(2)*exp(-0.02*x(3)*tdata)- cdata,2)输入命令: x0=0.2,0.05,0.05;x=lsqnonlin(curvefun2,x0)f= curvefun2(x),函数curvefun2的自变量是x,cdata和tdata是已知参数,故应将cdata tdata的值写在curvefun2.m中,解法 2 用命令lsqnonlin,x=(a,b,k),3)运算结果为 f =1.0e-003 *(0.2322 -0.1243 -0.2495 -0.2413 -0.1668 -0.0724 0.0241 0.1159 0.2030 0.2792) x =0.0063 -0.0034 0.2542,可以看出,两个命令的计算结果是相同的。,4)结论:即拟合得a=0.0063 b=-0.0034 k=0.2542,说明:拟合与统计回归,区别与联系,统计回归线性回归 (regress命令)非线性回归,非线性回 归,(1)确定回归系数的命令: beta,r,J=nlinfit(x,y,model, beta0),(2)非线性回归命令:nlintool(x,y,model, beta0,alpha),1、回归:,例题的求解:,2、输入数据: x=2:16; y=6.42 8.20 9.58 9.5 9.7 10 9.93 9.99 10.49 10.59 10.60 10.80 10.60 10.90 10.76; beta0=8 2;,3、求回归系数: beta,r ,J=nlinfit(x,y,volum,beta0); beta,得结果:beta = 11.6036 -1.0641,即得回归模型为:,To MATLAB(liti41),4、预测及作图: YY,delta=nlpredci(volum,x,beta,r ,J); plot(x,y,k+,x,YY,r),插值问题,拟合与插值的关系,说明:函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似,由于近似的要求不同,二者的数学方法上完全不同。,实例:下面数据是某次实验所得,希望得到x和 f之间的关系?,MATLAB(cn),问题:给定一批数据点,需确定满足特定要求的曲线或曲面,解决方案:,若不要求曲线(面)通过所有数据点,而是要求它反映对象整体的变化趋势,就是数据拟合,又称曲线拟合或曲面拟合。,若要求所求曲线(面)通过所给所有数据点,就是插值问题;,最临近插值、线性插值、样条插值与曲线拟合结果:,拉格朗日插值,分段线性插值,三次样条插值,一 维 插 值,一、插值的定义,二、插值的方法,三、用Matlab解插值问题,返回,返回,二维插值,一、二维插值定义,二、网格节点插值法,三、用Matlab解插值问题,最邻近插值,分片线性插值,双线性插值,网格节点数据的插值,散点数据的插值,一维插值的定义,节点可视为由,产生,,,表达式复杂,,,或无封闭形式,,,或未知.。,返回,称为拉格朗日插值基函数。,已知函数f(x)在n+1个点x0,x1,xn处的函数值为 y0,y1,yn 。求一n次多项式函数Pn(x),使其满足: Pn(xi)=yi,i=0,1,n.,解决此问题的拉格朗日插值多项式公式如下,其中Li(x) 为n次多项式:,拉格朗日(Lagrange)插值,拉格朗日(Lagrange)插值,特别地:,两点一次(线性)插值多项式:,三点二次(抛物)插值多项式:,拉格朗日多项式插值的这种振荡现象叫 Runge现象,采用拉格朗日多项式插值:选取不同插值节点个数n+1,其中n为插值多项式的次数,当n分别取2,4,6,8,10时,绘出插值结果图形.,例,返回,To Matlablch(larg1),分段线性插值,计算量与n无关;n越大,误差越小.,To MATLABxch11,xch12,xch13,xch14,返回,例,用分段线性插值法求插值,并观察插值误差.,1.在-6,6中平均选取5个点作插值(xch11),4.在-6,6中平均选取41个点作插值(xch14),2.在-6,6中平均选取11个点作插值(xch12),3.在-6,6中平均选取21个点作插值(xch13),比分段线性插值更光滑。,在数学上,光滑程度的定量描述是:函数(曲线)的k阶导数存在且连续,则称该曲线具有k阶光滑性。 光滑性的阶次越高,则越光滑。是否存在较低次的分段多项式达到较高阶光滑性的方法?三次样条插值就是一个很好的例子。,三次样条插值,g(x)为被插值函数。,三次样条插值,例,用三次样条插值选取11个基点计算插值(ych),返回,To MATLABych(larg1),用MATLAB作插值计算,一维插值函数:,yi=interp1(x,y,xi,method),nearest :最邻近插值linear : 线性插值;spline : 三次样条插值;cubic : 立方插值。缺省时: 分段线性插值。,注意:所有的插值方法都要求x是单调的,并且xi不能够超过x的范围。,例:在1-12的11小时内,每隔1小时测量一次温度,测得的温度依次为:5,8,9,15,25,29,31,30,22,25,27,24。试估计每隔1/10小时的温度值。,To MATLAB(temp),hours=1:12;temps=5 8 9 15 25 29 31 30 22 25 27 24;h=1:0.1:12;t=interp1(hours,temps,h,spline); (直接输出数据将是很多的)plot(hours,temps,+,h,t,hours,temps,r:) %作图xlabel(Hour),ylabel(Degrees Celsius),例 已知飞机下轮廓线上数据如下,求x每改变0.1时的y值。,To MATLAB(plane),返回,二维插值的定义,第一种(网格节点):,已知 mn个节点,第二种(散乱节点):,返回,注意:最邻近插值一般不连续。具有连续性的最简单的插值是分片线性插值。,最邻近插值,二维或高维情形的最邻近插值,与被插值点最邻近的节点的函数值即为所求。,返回,将四个插值点(矩形的四个顶点)处的函数值依次简记为:,分片线性插值,f (xi, yj)=f1,f (xi+1, yj)=f2,f (xi+1, yj+1)=f3,f (xi, yj+1)=f4,插值函数为:,第二片(上三角形区域):(x, y)满足,插值函数为:,注意:(x, y)当然应该是在插值节点所形成的矩形区域内。显然,分片线性插值函数是连续的;,分两片的函数表达式如下:,第一片(下三角形区域): (x, y)满足,返回,双线性插值是一片一片的空间二次曲面构成。双线性插值函数的形式如下:,其中有四个待定系数,利用该函数在矩形的四个顶点(插值节点)的函数值,得到四个代数方程,正好确定四个系数。,双线性插值,返回,要求x0,y0单调;x,y可取为矩阵,或x取行向量,y取为列向量,x,y的值分别不能超出x0,y0的范围。,z=interp2(x0,y0,z0,x,y,method),用MATLAB作网格节点数据的插值,nearest 最邻近插值linear 双线性插值cubic 双三次插值缺省时, 双线性插值,例:测得平板表面3*5网格点处的温度分别为: 82 81 80 82 84 79 63 61 65 81 84 84 82 85 86 试作出平板表面的温度分布曲面z=f(x,y)的图形。,输入以下命令:x=1:5;y=1:3;temps=82 81 80 82 84;79 63 61 65 81;84 84 82 85 86;mesh(x,y,temps),1.先在三维坐标画出原始数据,画出粗糙的温度分布曲图.,2以平滑数据,在x、y方向上每隔0.2个单位的地方进行插值.,再输入以下命令:xi=1:0.2:5;yi=1:0.2:3;zi=interp2(x,y,temps,xi,yi,cubic);mesh(xi,yi,zi)画出插值后的温度分布曲面图.,To MATLAB(wendu),通过此例对最近邻点插值、双线性插值方法和双三次插值方法的插值效果进行比较。,To MATLAB (moutain),返回,插值函数griddata格式为:,cz =griddata(x,y,z,cx,cy,method),用MATLAB作散点数据的插值计算,要求cx取行向量,cy取为列向量。,nearest 最邻近插值linear 双线性插值cubic 双三次插值v4- Matlab提供的插值方法缺省时, 双线性插值,最优化,最佳水槽断面问题,(矩形断面) 用宽 l = 24 cm的长方铁板折成一个断面为矩形的水槽,问怎样的折法可使水槽的断面面积达到最大,最佳水槽断面问题,(梯形断面) 将问题推广等腰梯形的水槽,问怎样折法可使水槽断面面积达到最大?,最佳水槽断面问题,(对称五边形断面)将铁板折成如图所示的对称五边形,问怎样的折法可使水槽的断面面积达到最大?,最佳水槽断面问题(五边),最佳水槽断面问题(五边),运行zxy6_6s,fval = fmincon(zxy6_6S,x0,A,b,lb,ub)求解,最佳水槽断面问题(多边和无限边),优化变量数与最大断面面积的关系断面形状 优化变量数 最大断面积 cm2矩形断面 1 72梯形断面 2 83.14对称五边形 4 88.637将铁板折成对称7边形、9边形,一般为对称2n+1边形可以期望最大断面面积得到进一步的增加随之而来是计算代价也随之增加。,最佳水槽断面问题(无限边),最佳水槽断面问题(无限边),最佳水槽断面问题(无限边),微分法求最大和最小(高等数学),微分法求最大和最小(高等数学),微分法求最大和最小(高等数学),运行zxy6_1.m syms x1 x2 %定义符号变量。f=x13-x23+3*x12+3*x22-9*x1; % 函数z。v=x1 x2;df=jacobian(f,v) %计算雅可比。 X,Y=solve(df(1),df(2) % 用指令solve求驻点。运行zxy6_2画等值线图并将点标注在图上,微分法求最大和最小(高等数学),jacobian(f,v):计算函数的符号梯度、雅可比矩阵如:若f(v),v=v1 v2则df=df/dv1 df/dv2如:若f(v)=f1(v) f2(v) ,v=v1 v2则df=df1/dv1 df1/dv1 df2/dv1 df2/dv2,微分法求最大和最小(高等数学),solve指令:solve(eqn1,eqn2,eqnn)求n个方程eqn1,eqn2, ,eqnn所构成的方程组的根(符号解),盲人下山与迭代寻优,一个盲人处于山上的某一点x0要下到谷底,他应如何做?,盲人下山与迭代寻优,Matlab优化工具箱简介,多元函数无约束优化指令fminunc、fminsearch的剖析,Matlab优化工具箱简介,观察:在命令窗口键入bandemo选择不同方法观察对香蕉函数的优化结果和过程。,Matlab优化工具箱简介,x,fval,exitflag,output,grad,hessian= fminunc(fun,x0,options,P1,)其中有些项可以缺省,如exitflag,output,grad,hessian,options,P1,P2, 等等。x0是初始点;fun是目标函数,可以用inline指令或建立M文件的方法生成目标函数;,Matlab优化工具箱简介,参数的传递:使用变量P1,P2,可在目标函数和主程序之间需要传递某些参数也可使用全局变量Gobal说明来进行传递。输出:x,fval,exitflag,output,grad,hessian为输出信息,最优点、最优函数值、算法结束的状态 (exitflag 0 算法收敛;=0达到最大步骤而停止;0算法不收敛)、算法结束后的某些信息(如迭代次数、所使用的优化方法等,可在命令窗口查看output的具体内容)、梯度值和海森