数据拟合与模型选择

1、 数据拟合与模型选择数据拟合与模型选择Curve Fit and Model SelectingJie FeiMPTCJune 20, 2012引例引例 美国人口预测美国人口预测 给出美国人口从给出美国人口从1790年到年到1990年间的人口如表年间的人口如表1（每（每10年为年为一个间隔），请估计出美国一个间隔），请估计出美国2010年的人口。年的人口。年 份1790180018101820183018401850人口(106)3.95.37.29.612.917.123.2年 份1860187018801890190019101920人口(106)31.438.650.262.976.09

2、2.0105.8年 份1930194019501960197019801990人口(106)122.8131.7150.7179.3203.2226.5248.7表1 美国人口统计数据 一、数据拟合的方法一、数据拟合的方法1.平面上绘出已知数据的分布图平面上绘出已知数据的分布图 (散点图，散点图，Scatter Plot).2.通过直观观察或经验公式猜测人口随时间的变化规律通过直观观察或经验公式猜测人口随时间的变化规律 (函数关系函数关系).3. 利用函数拟合的方法确定拟合函数中的未知参数利用函数拟合的方法确定拟合函数中的未知参数.4. 利用拟合函数估计出利用拟合函数估计出2010年的美国人口

3、年的美国人口.曲线拟合问题的提法曲线拟合问题的提法什么是曲线拟合什么是曲线拟合 已知一组二元数据，已知一组二元数据， 这组数据形成这组数据形成 平面上的一组散点平面上的一组散点. 在在某一类函数某一类函数中寻找一个函数中寻找一个函数 使得函数曲线在使得函数曲线在某种准则某种准则下与所有数据点最为接近，下与所有数据点最为接近， 这种构造近似函数的方法称为曲线拟合，这种构造近似函数的方法称为曲线拟合， 称为拟合函数称为拟合函数.( ,), 1,2,., ,iix yin( )f x( )f xxyO曲线拟合的图示曲线拟合的图示图图1 数据散点图数据散点图设有设有 m 个数据点（个数据点（xi, y

4、i）i=1, 2, , m. 作作散点图散点图如下：如下：Scatter Plot 要想对图要想对图1所示的数据拟合模型所示的数据拟合模型 f (x) = a x + b. 应如何选应如何选择择 a 和和 b ，使直线最好地拟合数据？即最佳拟合的原则与方，使直线最好地拟合数据？即最佳拟合的原则与方法是什么？法是什么？最佳拟合准则最佳拟合准则 从图上看，存在两个以上点时，不能期望它们精确地处从图上看，存在两个以上点时，不能期望它们精确地处于一直线上。数据点和直线间总存在一些纵向差异。称这于一直线上。数据点和直线间总存在一些纵向差异。称这些纵向差异为些纵向差异为绝对偏差绝对偏差 。最佳拟合的准则：

5、最佳拟合的准则：（1）极小化这些绝对偏差的和）极小化这些绝对偏差的和. 即极小化即极小化 图图2 极小化绝对偏差的和极小化绝对偏差的和解这一问题必须用数值最优化方法解这一问题必须用数值最优化方法解出模型参数解出模型参数a, b的估计的估计.1( )miiiyf x偏差偏差（2）极小化所有点的最大偏差）极小化所有点的最大偏差. 图图3 极小化最大绝对偏差极小化最大绝对偏差（3） 极小化这些绝对偏差的平方和（极小化这些绝对偏差的平方和（最小二乘原则最小二乘原则）即极小化即极小化Max| yi f (xi) |， i =1,2, ,m解这个问题可能需要高级的数学解这个问题可能需要高级的数学方法，或者

6、需要计算机的数值算方法，或者需要计算机的数值算法。简单情形，会转化成一个线法。简单情形，会转化成一个线性规划问题。性规划问题。偏差偏差（Least-Squares Criterion）即极小化即极小化21( )miiiyf x21( )miiiyf x 最小二乘准则最小二乘准则函数中的参数作为自变量，绝对偏差的平方和作为目标函函数中的参数作为自变量，绝对偏差的平方和作为目标函数，利用多元函数的极值理论就可以解决。数，利用多元函数的极值理论就可以解决。给定某一函数类型给定某一函数类型 y = f (x) , 以及以及m个数据点（个数据点（xi, yi）的集）的集合，极小化绝对偏差合，极小化绝对偏

7、差 | yi - f(xi)| 的平方和，即确定函数的平方和，即确定函数 y = f (x)中的参数，极小化中的参数，极小化21( )miiiyf x21( )miiiyf x用最小二乘准则来估计各种类型曲线参数的数学过程如下：用最小二乘准则来估计各种类型曲线参数的数学过程如下：一、拟合直线一、拟合直线设预期模型的形式为设预期模型的形式为 y =Ax+B，用，用a、b记记A、B的最小二乘的最小二乘估计，这时要求极小化：估计，这时要求极小化：2211( )mmiiiiiiSyf xyaxb等价于下述线性方程组的最小二乘解：等价于下述线性方程组的最小二乘解：1122 mmaxbyaxbyaxby根

8、据多元函数的极值理论，最优的必要条件是两个偏导数根据多元函数的极值理论，最优的必要条件是两个偏导数等于等于0.21212()02()0miiiimiiiSx yaxbaSyaxbb 211111mmmiiiiiiimmiiiiaxbxx yaxmby整理得正规方程：整理得正规方程：（正规方程）（正规方程）正规方程组解的矩阵表示：正规方程组解的矩阵表示：1211,1mxxAx 12,myyyy 记记若若 可逆，正规方程的解为可逆，正规方程的解为 .ba TAA().TTAAA y 则正规方程组为则正规方程组为1().TTAAA y 解得解得1112211211112211,mmmiiiiiiim

9、miiiimmmmiiiiiiiiimmiiiimx yxyamxxxyx yxbmxx可以编写计算机程序解可以编写计算机程序解 a 和和 b. 用用Matlab编程，只要输入数据，编程，只要输入数据，再输入再输入 regress命令即可得到系数的最小二乘估计。命令即可得到系数的最小二乘估计。统计上利用最小二乘准则估计拟合直线的参数称为统计上利用最小二乘准则估计拟合直线的参数称为线性回归线性回归。例例1：弹簧的弹性系数的测定：弹簧的弹性系数的测定.下表给出了悬挂不同重量的物体时弹簧的长度，试由这下表给出了悬挂不同重量的物体时弹簧的长度，试由这组数据推测弹簧的弹性系数组数据推测弹簧的弹性系数.已

10、知弹簧的弹力与弹簧长度有线性关系：已知弹簧的弹力与弹簧长度有线性关系： y=ax+b其中常数其中常数 a与与b需要根据数据测定需要根据数据测定. 51015202530长度(cm)7.258.128.959.9010.9011.80 x=5:5:30;y=7.25 8.12 8.95 9.90 10.90 11.80;x1=ones(length(x),1),x;regress(y,x1)输入输入 Matlab 指令如下指令如下： 得到最佳拟合直线为：得到最佳拟合直线为： y= 0.1831 x+ 0.6827 51015202530拟合效果图拟合效果图：二、拟合多项式二、拟合多项式最小二乘估

11、计即要求极小化（以最小二乘估计即要求极小化（以n=2为为例）为为例）S对参数求偏导，令其等于零得方程：对参数求偏导，令其等于零得方程：设预期模型的形式为设预期模型的形式为 n是固定的是固定的. 01( ),nnyf xaa xa x 2 20121()miiSyaa xa x 22012102201211222012102()02()02()0miiiimiiiiimiiiiiSyaa xa xaSx yaa xa xaSxyaa xa xa 正规方程组解的矩阵表示：正规方程组解的矩阵表示：211222211,1mmxxxxAxx 12,myyyy 记记若若 可逆，正规方程的解为可逆，正规方程

12、的解为 012.aaaa TAA().TTAA aA y 则正规方程组为则正规方程组为1().TTaAAA y 多项式拟合的多项式拟合的Matlab指令为：指令为：polyfit 其调用格式为其调用格式为 a=polyfit(adata, ydata, n)其中其中 n拟合多项式的次数拟合多项式的次数. xdata, ydate拟合的数据，以数组方式输入拟合的数据，以数组方式输入. a输出参数，拟合多项式的系数输出参数，拟合多项式的系数.多项式在多项式在 x处的取值处的取值y可用如下命令格式计算：可用如下命令格式计算： y=polyval(a, x)三、经变换的最小二乘拟合三、经变换的最小二乘

13、拟合这个方程组是非线性的，手工难于求解。这个方程组是非线性的，手工难于求解。理论上最小二乘准则可用于其它模型，方法都是对参数求理论上最小二乘准则可用于其它模型，方法都是对参数求导，令其等于零，解得到的方程，求出模型的参数。但在导，令其等于零，解得到的方程，求出模型的参数。但在实践上可能有困难。实践上可能有困难。例如要拟合模型例如要拟合模型 ,作最小二乘估计时需要极小化作最小二乘估计时需要极小化 bxyae2211( )immbxiiiiiSyf xyae120iimbxbxiiSeyaea 120iimbxbxiiiSaexyaeb 处理方法：在模型处理方法：在模型 取对数以后，变成线性模型：

14、取对数以后，变成线性模型：bxyaelnln.yabx 可以利用最小二乘求解。可以利用最小二乘求解。非线性关系变成线性关系，常见形式有一下几种：非线性关系变成线性关系，常见形式有一下几种：111. ,xyyxuabvaxbuv 2. ,ln ,lnln,byax uy vxuabv 3. ,lnln,bxyaeuyuabx 14. ,ln ,ln,bxyaeuy vuabvx 115. ,xxyyveuabvabeu 一般地，变换以后再进行最小二乘拟合与直接进行最小一般地，变换以后再进行最小二乘拟合与直接进行最小二乘拟合的结果是不一样的。二乘拟合的结果是不一样的。 在原始问题中，寻找曲线时，是

15、极小化原始数据的偏差在原始问题中，寻找曲线时，是极小化原始数据的偏差平方和，而在变换后的问题中，极小化变换后的变量的偏差平方和，而在变换后的问题中，极小化变换后的变量的偏差平方和。平方和。非线性拟合的非线性拟合的Matlab指令为：指令为：lsqcurvefit 其调用格式为其调用格式为 p= lsqcurvefit(Fun, p0 , xdata, ydata)其中其中 Fun表示拟合函数的表示拟合函数的M文件或内联函数文件或内联函数. xdata, ydate拟合的数据，以数组方式输入拟合的数据，以数组方式输入. p0 拟合参数的初值拟合参数的初值. p拟合的参数拟合的参数.若要求在若要求

16、在 x处的取值处的取值y可用如下命令格式计算：可用如下命令格式计算： y=Fun(p, x) 2. 按照一个或一些选出的模型类型对数据进行拟合按照一个或一些选出的模型类型对数据进行拟合 (确定模型中的最佳参数确定模型中的最佳参数).3. 对拟合模型的评价对拟合模型的评价 (不同模型之间的优劣不同模型之间的优劣).4. 根据收集的数据做出预报根据收集的数据做出预报, 评价预报的效果评价预报的效果.在分析一个数据集合时，有四个可能需要解决的任务：在分析一个数据集合时，有四个可能需要解决的任务：1. 模型选择模型选择 , 线性模型线性模型, 多项式模型多项式模型, 指数模型或是其它模型指数模型或是其

17、它模型.怎样评价一个好模型？怎样评价一个好模型？评估模型是否很好地拟合了数据的途径是：计算残差。即评估模型是否很好地拟合了数据的途径是：计算残差。即计算模型点与实际点间的偏差。计算模型点与实际点间的偏差。要回答哪个模型最好时要以考虑模型的目的、实际情况要要回答哪个模型最好时要以考虑模型的目的、实际情况要求的精度、数据的准确性以及使用模型时独立变量的值的求的精度、数据的准确性以及使用模型时独立变量的值的范围等因素。范围等因素。1. 单项模型的选择；单项模型的选择；2. 通过数据点的高阶多项式；通过数据点的高阶多项式；3. 低阶多项式对数据光滑化；低阶多项式对数据光滑化；4. 三次样条插值。三次样

18、条插值。当实际问题中不能根据一些假定提出某种模型时，就要基当实际问题中不能根据一些假定提出某种模型时，就要基于数据建立经验模型。于数据建立经验模型。 二、基于数据建立经验模型二、基于数据建立经验模型例例2. Chesapeake 海湾的收成海湾的收成.1. 单项模型的选择单项模型的选择任务是预测蓝鱼的产量。任务是预测蓝鱼的产量。策略是变换数据，使得所产生的图形近似一条直线。策略是变换数据，使得所产生的图形近似一条直线。怎样确定这一变换呢？怎样确定这一变换呢？作出散点图：作出散点图：散点图中数据的倾向是增的、下凸的。散点图中数据的倾向是增的、下凸的。 对于一个向上凹的正值函数，对于一个向上凹的正

19、值函数，y = f ( x )，x 1.用阶梯中处于用阶梯中处于 z下方的某些变换，将下方的某些变换，将y值变为值变为 或或 log y 或更剧烈的变化，挤压右侧尾部向下，或更剧烈的变化，挤压右侧尾部向下，可能产生更接近直线的新函数。应该采用哪个可能产生更接近直线的新函数。应该采用哪个变换是反复试验、不断摸索的问题（或根据经变换是反复试验、不断摸索的问题（或根据经验）。另一种变换是改变验）。另一种变换是改变 x 的值为的值为 x2，x3等，等，拉伸右侧的尾部向右。拉伸右侧的尾部向右。y22 log1 1 1 zzzzzzz幂阶梯对于蓝鱼的产量数据，将对于蓝鱼的产量数据，将 x 的值改为阶梯向上

20、的几种值的值改为阶梯向上的几种值（x2，x3 等），不能产生有关线性图形；因此将等），不能产生有关线性图形；因此将y的值改的值改为阶梯向下的值为阶梯向下的值 或或 ln y 的值。经比较，选取的值。经比较，选取 ln y 对对 x 的模型，用最小二乘拟合下面模型：的模型，用最小二乘拟合下面模型：ylog yaxb 找出拟合的曲线找出拟合的曲线 其中其中x是基底年，是基底年，log y 是以是以10为底的对数，为底的对数，y 的单位是的单位是104磅。磅。log0.7231 0.1654yx 由由log y取指数得取指数得5.2857(1.4635)xy 数据拟合效果图：数据拟合效果图：例例3.

21、 估计蓝蟹估计蓝蟹.散点图中数据的倾向是增的、上凸的。散点图中数据的倾向是增的、上凸的。作出散点图：作出散点图： 对于一个向上凸的正值函数，对于一个向上凸的正值函数，y = f ( x )，x 1.这次尝试改变这次尝试改变 y 的值成为的值成为 y2、y3 等来拉伸右侧等来拉伸右侧尾部向上而线性化。另一种可能是尝试改变尾部向上而线性化。另一种可能是尝试改变x的的值成为值成为 或或log x ，或更剧烈的变换梯步（变，或更剧烈的变换梯步（变换阶梯表的下部），来使右侧尾部向左。换阶梯表的下部），来使右侧尾部向左。 x注意，阶梯表中注意，阶梯表中log 以下的变换有时会把一个增以下的变换有时会把一个

22、增函数逆转为减函数，这是不希望的。因此数据函数逆转为减函数，这是不希望的。因此数据分析时常加一个负号，以保持变换的数据与原分析时常加一个负号，以保持变换的数据与原数据有相同的增减性。数据有相同的增减性。322 log1 1 1 zzzzzzzz变换阶梯常用变换常用变换试图改变试图改变y的值为的值为 y2 或或 y3 和改变和改变x的值成为的值成为 或或log x 等等，最后选择一条过原点的直线最后选择一条过原点的直线ykx由最小二乘法得到由最小二乘法得到158.344yx其中其中 x 是基底年，是基底年，y的单位是的单位是104磅蓝蟹。磅蓝蟹。x曲线似乎是可取的。曲线似乎是可取的。曲线拟合效果

23、图：曲线拟合效果图：蓝鱼蓝鱼蓝蟹蓝蟹模型验证：对每一对数据计算残差和相对误差。模型验证：对每一对数据计算残差和相对误差。模型预测：预测模型预测：预测2010年海湾收成年海湾收成1445.2857(1.4635)1092.95(10)10.9.y 磅百万磅4158.344 14592.469(10)5.92.y 磅百万磅这些简单的单项模型应该用于内插，而在作外推时误差较大。这些简单的单项模型应该用于内插，而在作外推时误差较大。小结小结：构造一个预测模型时，细心分析收集到的数据，看：构造一个预测模型时，细心分析收集到的数据，看数据存在什么样的倾向？是否有明显处于倾向外的数据点？数据存在什么样的倾向

24、？是否有明显处于倾向外的数据点？如果这样的异常值存在，是否抛弃它？如果是实验观察到如果这样的异常值存在，是否抛弃它？如果是实验观察到的，重复该实验以检查数据。当某一种倾向确实清楚存在的，重复该实验以检查数据。当某一种倾向确实清楚存在时，找到一个将数据变换成一直线（近似地）的函数。时，找到一个将数据变换成一直线（近似地）的函数。如果选择了模型如果选择了模型 y = a xb, 应画一个应画一个lny对对 lnx 的图，看是否的图，看是否产生一近似直线。同样在研究模型产生一近似直线。同样在研究模型y = a ebx 的适当性时，的适当性时，应画应画ln y 对对x 的图是否产生一条直线。当认为一个

25、确定的模的图是否产生一条直线。当认为一个确定的模型类型比较符合数据的倾向时，可以用最小二乘及其他解型类型比较符合数据的倾向时，可以用最小二乘及其他解析方法估计模型参数。最后对原始模型作拟合优度分析。析方法估计模型参数。最后对原始模型作拟合优度分析。如果对这一拟合不满意，可以研究其他的单项模型。如果如果对这一拟合不满意，可以研究其他的单项模型。如果在单项模型框架内不能满意地拟合数据，则要使用其他的在单项模型框架内不能满意地拟合数据，则要使用其他的技术。这就是将要讨论的多项式拟合和样条插值。技术。这就是将要讨论的多项式拟合和样条插值。2. 高阶多项式高阶多项式(polynomials)模型模型单项

26、模型易于进行模型分析，然而由于其数学上的简单性，单项模型易于进行模型分析，然而由于其数学上的简单性，在追踪数据到的数据的趋势时，单项模型是有局限的。在某在追踪数据到的数据的趋势时，单项模型是有局限的。在某些情况下，必须考虑多于一项的模型，即多项式。当然多项些情况下，必须考虑多于一项的模型，即多项式。当然多项式也是有其缺点的。式也是有其缺点的。大家知道：有惟一的最高阶为大家知道：有惟一的最高阶为2的多项式的多项式 y = a0+ a1x + a2x2 能够通过三个不同的点。但存在无数个阶大于能够通过三个不同的点。但存在无数个阶大于2的多项式通的多项式通过这三个点。过这三个点。 例例4 带式录音机

27、的播放时间带式录音机的播放时间.收集一个特定的录音机的计数器的数据及相应的录音机的收集一个特定的录音机的计数器的数据及相应的录音机的播放时间。数据如下（播放时间。数据如下（ ci 表示计数器读数，表示计数器读数，ti 为对应的播为对应的播放总时间）放总时间）ci (读数读数) 100 200 300 400 500 600 700 800ti (秒秒) 205 430 677 945 1233 1542 1872 2224可以构造一个最高阶为可以构造一个最高阶为7的惟一多项式，通过这的惟一多项式，通过这8个数据个数据点。记此多项式为点。记此多项式为267701267( )P caa ca ca

28、 ca c因为要经过因为要经过8个点，所以有下面方程组（读数缩小个点，所以有下面方程组（读数缩小100倍）：倍）：解得：解得：对对P7(ci)的预测值舍入到小数第四位，给出了与观测到数据的预测值舍入到小数第四位，给出了与观测到数据完全一致的预测值，产生零绝对偏差。我们真的能够认为完全一致的预测值，产生零绝对偏差。我们真的能够认为这一模型比其他可以提供的模型更好吗？答案是否定的。这一模型比其他可以提供的模型更好吗？答案是否定的。后面将要谈到它的缺点。后面将要谈到它的缺点。在理论上可以证明在理论上可以证明：给定：给定(n+1)个不同的数据点，存在惟一个不同的数据点，存在惟一的一个最高阶为的一个最高

29、阶为n的多项式通过全部数据点。这就是多项的多项式通过全部数据点。这就是多项式的式的Lagrange形式。形式。 高阶多项式的优点与缺点：高阶多项式的优点与缺点：优点优点：因为多项式函数易于积分和微分。如果高阶多项式：因为多项式函数易于积分和微分。如果高阶多项式 模型对研究的问题是合适的，则后面的计算就比较容易。模型对研究的问题是合适的，则后面的计算就比较容易。缺点：缺点：虽然虽然 Lagrange 多项式能够通过给定的多项式能够通过给定的n个数据点，个数据点，但在区间的端点附近，多项式有严重的摆动。而在作预测但在区间的端点附近，多项式有严重的摆动。而在作预测和外推时，端点附近的估计往往是很需要

30、的。和外推时，端点附近的估计往往是很需要的。高阶多项式的另一个严重的缺点是它的系数对数据的微小高阶多项式的另一个严重的缺点是它的系数对数据的微小变化的敏感性。因此高阶多项式的这些缺点限制了它在建变化的敏感性。因此高阶多项式的这些缺点限制了它在建模中的应用。模中的应用。 例如，给定例如，给定17个数据点如下：个数据点如下：xi -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 yi 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0高阶多项式拟合图高阶多项式拟合图:x=-2:1:7; y=sin(x);p=polyfit(x,y,7);x1=-5

31、:0.05:10;f=polyval(p,x1);plot(x,y,o,x1,f,-)axis(-5 12 -5 5)又如，又如，设有数据集：设有数据集：xi= -2, -1, 0, 1, 2, 3, 4, 5, 6, 7. yi= sin xi共共8个数据点。用个数据点。用7阶多项式拟合，可以经过这阶多项式拟合，可以经过这8个点。个点。但在但在x 7 时，多项式的变动非常大。时，多项式的变动非常大。-4-2024681012-5-4-3-2-1012345 3. 低阶多项式对数据的光滑化低阶多项式对数据的光滑化 用低阶多项式拟合数据，不要求通过全部数据点。这用低阶多项式拟合数据，不要求通过全

32、部数据点。这将降低多项式的摆动倾向和它对数据中微小变化的敏感性。将降低多项式的摆动倾向和它对数据中微小变化的敏感性。低阶多项式拟合要求做出两个决定低阶多项式拟合要求做出两个决定: 1. 决定选择所用多项式的阶；决定选择所用多项式的阶； 2. 确定最佳拟合多项式的系数。确定最佳拟合多项式的系数。例例4 带式录音机的播放时间带式录音机的播放时间ci (读数读数) 100 200 300 400 500 600 700 800ti (秒秒) 205 430 677 945 1233 1542 1872 2224先用二阶多项式拟合带式录音机的播放时间，用最小二乘先用二阶多项式拟合带式录音机的播放时间，

33、用最小二乘法决定最佳系数法决定最佳系数.x=100 200 300 400 500 600 700 800;y=205 430 677 945 1233 1542 1872 2224;x1=ones(length(x),1),x,(x.2);regress(y,x1) % a=polyfit(x,y,2) 得到最佳拟合二次多项式为：得到最佳拟合二次多项式为： P2(c) = 0.14286 + 1.94226 c + 0.00105 c2 残差平方和为残差平方和为 0.9048.输入输入 Matlab 指令如下指令如下： 计算残差，得：计算残差，得：ci (读数读数) 100 200 300

34、400 500 600 700 800ti (秒秒) 205 430 677 945 1233 1542 1872 2224ti-P2(ci)0.167 -0.452 0.000 0.524 0.119 -0.214 -0.476 0.333这些残差与时间量的阶比较，相当小。这些残差与时间量的阶比较，相当小。 那么怎样选择低阶多项式的阶呢？几阶多项式为合适？那么怎样选择低阶多项式的阶呢？几阶多项式为合适？ 注意二次多项式函数的二阶导数是常数，三阶导数为零。注意二次多项式函数的二阶导数是常数，三阶导数为零。 对于离散的数据点集，可以计算对于离散的数据点集，可以计算差分或均差差分或均差。计算计算

35、yi+1 - yi , i=1,2,3,4, 构成一阶差分，记为构成一阶差分，记为D. D. 从从D D 列的相列的相继的一阶差分间算出二阶差分，记为继的一阶差分间算出二阶差分，记为D D2. 继续下去，直到继续下去，直到从从 n 个数据点算出个数据点算出D Dn-1. xi 0 2 4 6 8 yi 0 4 16 36 64 假设有数据集（虚拟）：假设有数据集（虚拟）：上述虚拟数据集中二阶差分是常数，三阶差分是零上述虚拟数据集中二阶差分是常数，三阶差分是零. 注意各注意各阶差分只是估计了各阶导数的分子。当阶差分只是估计了各阶导数的分子。当 xi 是不等跨距时，是不等跨距时，若要更细致的估计，

36、需要计算均差。若要更细致的估计，需要计算均差。 数数 据据 一阶均差一阶均差 二阶均差二阶均差 x1 y2 x2 y2 x3 y3 2121yyxx3232yyxx3221322131yyyyxxxxxx 计算上述虚拟数据集的均差：计算上述虚拟数据集的均差： 计算带式录音机数据的均差表：计算带式录音机数据的均差表：从上表看到，二阶均差基本上是常数，三阶均差为零（小从上表看到，二阶均差基本上是常数，三阶均差为零（小数点后第四位为零，不应期望准确地为零）数点后第四位为零，不应期望准确地为零）. 可以看出数据可以看出数据基本上是二次的，支持用二次多项式作为经验模型。基本上是二次的，支持用二次多项式作

37、为经验模型。例例5. 酵母培养物的增长。数据是随时间（小时）测量的酵母酵母培养物的增长。数据是随时间（小时）测量的酵母细胞的总数。数据与均差表如下：细胞的总数。数据与均差表如下：其均差表能够帮助决定低阶多项式是否能够提供一个满意其均差表能够帮助决定低阶多项式是否能够提供一个满意的经验模型。的经验模型。注意一阶均差其均差表注意一阶均差其均差表 到到 t =8小时一直增长，然后开始小时一直增长，然后开始下降。反映在下降。反映在 中表现为一连串的负号，指示了凹凸性的中表现为一连串的负号，指示了凹凸性的变化。这样，不能期望用仅有单个凹形的二次多项式追踪变化。这样，不能期望用仅有单个凹形的二次多项式追踪

38、这些数据的趋势。三阶均差这些数据的趋势。三阶均差 的正负号就没有一连串的的正负号就没有一连串的趋势。应该用三次多项式去拟合数据趋势。应该用三次多项式去拟合数据. 为了说明问题，我们也用二次的多项式拟合这组数据，最为了说明问题，我们也用二次的多项式拟合这组数据，最小二乘拟合结果如下：小二乘拟合结果如下： P = - 93.82 + 65.70 t - 1.12 t 2模型拟合很差，追踪数据趋势是失败的。模型拟合很差，追踪数据趋势是失败的。-505101520-400-2000200400600800二次多项式拟合二次多项式拟合(Matlab命令为命令为polyfit)的结果：的结果：-50510

39、15200100200300400500600700三次多项式拟合三次多项式拟合(Matlab命令为命令为polyfit)的结果：的结果：更好的模型应该是用逻辑斯蒂曲线拟合。更好的模型应该是用逻辑斯蒂曲线拟合。lnPrMtCMP00()00( )rM t tMPP tMPPe0.5470663.0220( )171.5763tP te0510152001002003004005006007004. 三阶样条模型三阶样条模型三阶样条插值是非常通用的建模技术。在三阶样条插值是非常通用的建模技术。在相邻的数据点间相邻的数据点间使用不同的三阶多项式使用不同的三阶多项式，追踪数据的趋势。既保证基本关，追

40、踪数据的趋势。既保证基本关系的特征，同时减少摆动的倾向和对数据变化的灵敏性。系的特征，同时减少摆动的倾向和对数据变化的灵敏性。模型进一步的讨论略。模型进一步的讨论略。 给出美国人口从给出美国人口从1790年到年到1990年间的人口如表年间的人口如表1（每（每10年为年为一个间隔），请估计出美国一个间隔），请估计出美国2010年的人口。年的人口。表表1 美国人口统计数据美国人口统计数据 三、模型举例三、模型举例年 份1790180018101820183018401850人口(106)3.95.37.29.612.917.123.2年 份1860187018801890190019101920人

41、口(106)31.438.650.262.976.092.0105.8年 份1930194019501960197019801990人口(106)122.8131.7150.7179.3203.2226.5248.7建模方法建模方法1： 作出散点图，由图可以发现美国人口的变化规律曲线近作出散点图，由图可以发现美国人口的变化规律曲线近似为一条指数函数曲线，因此我们假设美国的人口满足函数似为一条指数函数曲线，因此我们假设美国的人口满足函数关系关系 x = f(t), f(t) = ea+bt，a, b为待定常数为待定常数 .1750180018501900195020000501001502002

42、50300图图1：散点图：散点图输出输出 -35.4278 0.0208log( )35.42780.0208f tt x=1790:10:1990;y=3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9 76.0 92.0 105.8 122.8 131.7 150.7 179.3 203.2 226.5 248.7;plot(x,y,x,y,*)y1=log(y);x1=ones(length(x),1),x;b= inv(x1*x1)*(x1)*y1 %或者regress(y1,x1)两边取对数，得线性模型：两边取对数，得线性模型：log

43、( )f tabt输入输入 Matlab 指令如下指令如下：即有变换以后的最小二乘估计即有变换以后的最小二乘估计为为还原为还原为35.4278 0.0208( )tf te检验检验拟合效果拟合效果：很不理想很不理想.x2=1780:2000;y2=exp(b(1)+b(2).*x2);plot(x,y,+,x2,y2)y3=exp(b(1)+b(2).*x); sum(y-y3).2) %残差175018001850190019502000050100150200250300350400450500输出输出残差平方和残差平方和(很大很大)： Res1=24917图图2：线性化后拟合效果图：线性

44、化后拟合效果图 这个这个初值非常重要初值非常重要，因此前面变换数据以后最小二乘方，因此前面变换数据以后最小二乘方法得到的估计常常作为现在法得到的估计常常作为现在 lsqcurvefit 的初值。的初值。 为了改进拟合效果，直接用为了改进拟合效果，直接用matlab 软件的软件的lsqcurvefit命令命令, 进行曲线拟合（数值方法求进行曲线拟合（数值方法求 的的最小值点最小值点 ，其中，其中x = f (t), f (t) = ea+bt， ），不过计算时需），不过计算时需要提供被估计参数的初值。要提供被估计参数的初值。niiixtfbaE12)(),(x=1790:10:1990;y=3.

45、9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9 76.0 92.0 105.8 122.8 131.7 150.7 179.3 203.2 226.5 248.7;fun = inline(exp(b(1)+b(2).*x),b,x) b0=-35,0.02;b1 Res2=lsqcurvefit(fun,b0,x,y) b1 = -23.7274 0.0147Res2 = 1.8766e+003输出输出23.7274 0.0147( )tf te得到直接进行曲线拟合的最小二乘估计为得到直接进行曲线拟合的最小二乘估计为预测美国在预测美国在201

46、0年的人口数为年的人口数为356.254百万。百万。残差平方和残差平方和: Res2=1876.6，有所减少，仍比较大。，有所减少，仍比较大。输入输入 Matlab 指令如下指令如下：x2=1780:2000;y3=exp(b1(1)+b1(2).*x2); plot(x,y,+,x2,y3)175018001850190019502000050100150200250300350检验检验拟合效果拟合效果：有所改进有所改进.图图3 lsqcurvefit拟合效果图拟合效果图建模方法建模方法2： 上述模型对过去的统计数据吻合得较好，但也存在问上述模型对过去的统计数据吻合得较好，但也存在问题，即人

47、口是呈指数规律无止境地增长，此时人口的自然题，即人口是呈指数规律无止境地增长，此时人口的自然增长率随人口的增长而增长，这不可能。一般说来，当人增长率随人口的增长而增长，这不可能。一般说来，当人口较少时增长得越来越快，即增长率在变大；人口增长到口较少时增长得越来越快，即增长率在变大；人口增长到一定数量以后，增长就会慢下来，即增长率变小一定数量以后，增长就会慢下来，即增长率变小。这是因这是因为，自然资源、环境条件等因素不允许人口无限制地增长，为，自然资源、环境条件等因素不允许人口无限制地增长，它们对人口的增长起着阻滞作用，而且随着人口的增加，它们对人口的增长起着阻滞作用，而且随着人口的增加，阻滞作

48、用越来越大。而且阻滞作用越来越大。而且人口最终会饱和，趋于某一个常人口最终会饱和，趋于某一个常数数x ，我们假设人口的，我们假设人口的净增长率为净增长率为 r x(t)(1-x(t)/x )，即人，即人口的静增长率随着人口的增长而不断减小，当口的静增长率随着人口的增长而不断减小，当 t时，净时，净增长率趋于零。增长率趋于零。按照这个假设，得到按照这个假设，得到00(1)( )dxxrxdtxx tx这便是丹麦数学家这便是丹麦数学家 Verhulst（1804-1849）于）于19世纪中叶提出世纪中叶提出的阻滞增长模型（的阻滞增长模型（Logistic模型）。模型）。 满足满足 x(1790)=

49、3.9 的微分方程的微分方程(1)的解为的解为(1)rtexxx)1790(3910) 1(1(2)fun = inline(b(1)./(1+(10/39*b(1)-1)*exp(1790-x)*b(2), b, x) b0= 200 0.1;% (猜测)b1 Res3=lsqcurvefit(fun,b0,x,y) b1 = 308.5182 0.0280Res2 = 813.8362输出输出残差平方和残差平方和为为 Res3=813.8362，已经比指数模型时减少很多。，已经比指数模型时减少很多。参数估计为参数估计为 r =0.0280, x =308.5128，2010年美国的人口年美

50、国的人口预计为预计为265.1068百万人。这个结果还比较合理，当百万人。这个结果还比较合理，当t时，时，静增长率趋于零，人口数趋于静增长率趋于零，人口数趋于308.5182百万人。百万人。输入输入 Matlab 指令如下指令如下：图图4 b0=200,0.1时的拟合效果时的拟合效果x3=1780:2020;y4= b1(1)./(1+(10/39*b1(1)-1)*exp(1790-x3)*b1(2); plot(x,y,+,x3,y4)检验拟合效果：效果比前面两种情形都好检验拟合效果：效果比前面两种情形都好.17501800185019001950200020500501001502002

51、50300回顾回顾 logistic 曲线的来源，它是下述微分方程的解：曲线的来源，它是下述微分方程的解：00(1)( )dxxrxdtxx tx(1)其中第二个方程是初始条件：其中第二个方程是初始条件： x(1790)=3.9由此确定了由此确定了logistic曲线中有两个参数曲线中有两个参数 r 和和 x .rtexxx)1790(3910) 1(1(2)显然显然 x(1790)=3.9 这个条件对预测这个条件对预测2010年的美国人口并不是年的美国人口并不是非常重要的。因此考虑没有初始条件的下述微分方程非常重要的。因此考虑没有初始条件的下述微分方程(1)dxxrxdtx(3)(1)dxx

52、rxdtMMx记记则方程则方程 (3) 为为()dxrx MxdtM()Mdxrdtx Mx分离变量得分离变量得11dxrdtxMx11dxrdtxMxlnln()xMxrtClnxrtCMx(4)因此方程因此方程 (3) 的解可表示为的解可表示为 (4)，其中，其中C是积分常数，是积分常数，可以由初始条件确定。我们可以由初始条件确定。我们不急于确定它不急于确定它，而是看作而是看作一个参数一个参数。这样我们的曲线。这样我们的曲线 (4) 有三个参数有三个参数.lnxyrtCMx(4) 把把 (4) 的右端看成数据变换的结果的右端看成数据变换的结果 y：则则 y 与与 t 的关系是线性的。如果能

53、预先猜测美国的极的关系是线性的。如果能预先猜测美国的极限人口限人口 M，则可以用线性模型的最小二乘估计方法得，则可以用线性模型的最小二乘估计方法得到净增长率到净增长率 r 的估计和的估计和 C 的估计。的估计。顺便说明顺便说明 logistic 一词的来源于一词的来源于 log it. 这充分说明了在这充分说明了在19世纪世纪（计算机技术远远没有现在这么发展的时代），对数变换是多（计算机技术远远没有现在这么发展的时代），对数变换是多么的重要！么的重要！ln53.34870.0274310 xtx 如果猜测美国的极限人口为如果猜测美国的极限人口为310（百万），则对（百万），则对1790年至年至

54、1990年美国人口数据年美国人口数据变换以后的最小二乘估计为变换以后的最小二乘估计为（不需要（不需要初值）初值） 它的等价形式是它的等价形式是53.3487 0.02743101txelnxyrtCMx方程方程 (4) 解解 的等价形式是的等价形式是1b rtMxe(5)fun1=inline(b(1)./ (1 + exp(b(2) +b(3)*x), b, x)b0=310 53.3 -0.027;b2 Res4=lsqcurvefit(fun1,b0,x,y)b2 = 387.5591 44.6944 -0.0227Res4 = 358.1678输出为输出为残差平方和为残差平方和为 Re

55、s4 = 358.1678，比较小。，比较小。 参数估计为参数估计为 r =0.0227, x =387.5591，C=44.6944，预计，预计2010年美国的人口为年美国的人口为 283.8656 百万人。这个结果还比较合百万人。这个结果还比较合理，当理，当t时，静增长率趋于零，人口数趋于时，静增长率趋于零，人口数趋于387.5591百百万人。万人。输入输入 Matlab 指令如下指令如下：x3=1780:2020;y5=b2(1) ./ (1 + exp(b2(2) +b2(3)*x3); plot(x,y,+,x3,y5)检验拟合效果：效果比前面几种情形都好检验拟合效果：效果比前面几种情形都好.1750180018501900195020