数据拟合与模型选择

1、数据拟合与模型选择 Curve Fit and Model Selecting,Jie Fei MPTC June 20, 2012,引例 美国人口预测,给出美国人口从1790年到1990年间的人口如表1（每10年为一个间隔），请估计出美国2010年的人口。,表1 美国人口统计数据,一、数据拟合的方法,平面上绘出已知数据的分布图 (散点图，Scatter Plot). 通过直观观察或经验公式猜测人口随时间的变化规律 (函数关系). 3. 利用函数拟合的方法确定拟合函数中的未知参数. 4. 利用拟合函数估计出2010年的美国人口.,曲线拟合问题的提法 什么是曲线拟合,已知一组二元数据， 这组数据

2、形成 平面上的一组散点. 在某一类函数中寻找一个函数 使得函数曲线在某种准则下与所有数据点最为接近， 这种构造近似函数的方法称为曲线拟合， 称为拟合函数.,曲线拟合的图示,图1 数据散点图,设有 m 个数据点（xi, yi）i=1, 2, , m. 作散点图如下：,Scatter Plot,要想对图1所示的数据拟合模型 f (x) = a x + b. 应如何选择 a 和 b ，使直线最好地拟合数据？即最佳拟合的原则与方法是什么？,最佳拟合准则,从图上看，存在两个以上点时，不能期望它们精确地处于一直线上。数据点和直线间总存在一些纵向差异。称这些纵向差异为绝对偏差 。,最佳拟合的准则： （1）极

3、小化这些绝对偏差的和. 即极小化,图2 极小化绝对偏差的和,解这一问题必须用数值最优化方法 解出模型参数a, b的估计.,（2）极小化所有点的最大偏差.,图3 极小化最大绝对偏差,（3） 极小化这些绝对偏差的平方和（最小二乘原则）,即极小化,Max| yi f (xi) |， i =1,2, ,m,解这个问题可能需要高级的数学方法，或者需要计算机的数值算法。简单情形，会转化成一个线性规划问题。,（Least-Squares Criterion）,即极小化,最小二乘准则,函数中的参数作为自变量，绝对偏差的平方和作为目标函数，利用多元函数的极值理论就可以解决。,给定某一函数类型 y = f (x)

4、 , 以及m个数据点（xi, yi）的集合，极小化绝对偏差 | yi - f(xi)| 的平方和，即确定函数 y = f (x)中的参数，极小化,用最小二乘准则来估计各种类型曲线参数的数学过程如下：,一、拟合直线,设预期模型的形式为 y =Ax+B，用a、b记A、B的最小二乘估计，这时要求极小化：,等价于下述线性方程组的最小二乘解：,根据多元函数的极值理论，最优的必要条件是两个偏导数等于0.,整理得正规方程：,（正规方程）,正规方程组解的矩阵表示：,记,若 可逆，正规方程的解为,则正规方程组为,解得,可以编写计算机程序解 a 和 b. 用Matlab编程，只要输入数据，再输入 regress命

5、令即可得到系数的最小二乘估计。,统计上利用最小二乘准则估计拟合直线的参数称为线性回归。,例1：弹簧的弹性系数的测定.,下表给出了悬挂不同重量的物体时弹簧的长度，试由这组数据推测弹簧的弹性系数.,已知弹簧的弹力与弹簧长度有线性关系： y=ax+b 其中常数 a与b需要根据数据测定.,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,拟合效果图：,二、拟合多项式,最小二乘估计即要求极小化（以

6、n=2为为例）,S对参数求偏导，令其等于零得方程：,设预期模型的形式为 n是固定的.,正规方程组解的矩阵表示：,记,若 可逆，正规方程的解为,则正规方程组为,多项式拟合的Matlab指令为：polyfit 其调用格式为 a=polyfit(adata, ydata, n),其中 n拟合多项式的次数. xdata, ydate拟合的数据，以数组方式输入. a输出参数，拟合多项式的系数.,多项式在 x处的取值y可用如下命令格式计算： y=polyval(a, x),三、经变换的最小二乘拟合,这个方程组是非线性的，手工难于求解。,理论上最小二乘准则可用于其它模型，方法都是对参数求导，令其等于零，解得

7、到的方程，求出模型的参数。但在实践上可能有困难。 例如要拟合模型 ,作最小二乘估计时需要极小化,处理方法：在模型 取对数以后，变成线性模型：,可以利用最小二乘求解。,非线性关系变成线性关系，常见形式有一下几种：,一般地，变换以后再进行最小二乘拟合与直接进行最小二乘拟合的结果是不一样的。,在原始问题中，寻找曲线时，是极小化原始数据的偏差平方和，而在变换后的问题中，极小化变换后的变量的偏差平方和。,非线性拟合的Matlab指令为：lsqcurvefit 其调用格式为 p= lsqcurvefit(Fun, p0 , xdata, ydata),其中 Fun表示拟合函数的M文件或内联函数. xdat

8、a, ydate拟合的数据，以数组方式输入. p0 拟合参数的初值. p拟合的参数.,若要求在 x处的取值y可用如下命令格式计算： y=Fun(p, x),2. 按照一个或一些选出的模型类型对数据进行拟合 (确定模型中的最佳参数).,3. 对拟合模型的评价 (不同模型之间的优劣).,4. 根据收集的数据做出预报, 评价预报的效果.,在分析一个数据集合时，有四个可能需要解决的任务：,1. 模型选择 , 线性模型, 多项式模型, 指数模型或是其它模型.,怎样评价一个好模型？,评估模型是否很好地拟合了数据的途径是：计算残差。即计算模型点与实际点间的偏差。 要回答哪个模型最好时要以考虑模型的目的、实际

9、情况要求的精度、数据的准确性以及使用模型时独立变量的值的范围等因素。,1. 单项模型的选择； 2. 通过数据点的高阶多项式； 3. 低阶多项式对数据光滑化； 4. 三次样条插值。,当实际问题中不能根据一些假定提出某种模型时，就要基于数据建立经验模型。,二、基于数据建立经验模型,例2. Chesapeake 海湾的收成.,1. 单项模型的选择,任务是预测蓝鱼的产量。,策略是变换数据，使得所产生的图形近似一条直线。,怎样确定这一变换呢？,作出散点图：,散点图中数据的倾向是增的、下凸的。,对于一个向上凹的正值函数，y = f ( x )，x 1.,用阶梯中处于 z下方的某些变换，将y值变为 或 lo

10、g y 或更剧烈的变化，挤压右侧尾部向下，可能产生更接近直线的新函数。应该采用哪个变换是反复试验、不断摸索的问题（或根据经验）。另一种变换是改变 x 的值为 x2，x3等，拉伸右侧的尾部向右。,对于蓝鱼的产量数据，将 x 的值改为阶梯向上的几种值（x2，x3 等），不能产生有关线性图形；因此将y的值改为阶梯向下的值 或 ln y 的值。经比较，选取 ln y 对 x 的模型，用最小二乘拟合下面模型：,找出拟合的曲线,其中x是基底年，log y 是以10为底的对数，y 的单位是104磅。,由log y取指数得,数据拟合效果图：,例3. 估计蓝蟹.,散点图中数据的倾向是增的、上凸的。,作出散点图：

11、,对于一个向上凸的正值函数，y = f ( x )，x 1.,这次尝试改变 y 的值成为 y2、y3 等来拉伸右侧尾部向上而线性化。另一种可能是尝试改变x的值成为 或log x ，或更剧烈的变换梯步（变换阶梯表的下部），来使右侧尾部向左。,注意，阶梯表中log 以下的变换有时会把一个增函数逆转为减函数，这是不希望的。因此数据分析时常加一个负号，以保持变换的数据与原数据有相同的增减性。,试图改变y的值为 y2 或 y3 和改变x的值成为 或log x 等，最后选择一条过原点的直线,由最小二乘法得到,其中 x 是基底年，y的单位是104磅蓝蟹。,曲线似乎是可取的。,曲线拟合效果图：,蓝鱼,蓝蟹,模

12、型验证：对每一对数据计算残差和相对误差。,模型预测：预测2010年海湾收成,这些简单的单项模型应该用于内插，而在作外推时误差较大。,小结：构造一个预测模型时，细心分析收集到的数据，看数据存在什么样的倾向？是否有明显处于倾向外的数据点？如果这样的异常值存在，是否抛弃它？如果是实验观察到的，重复该实验以检查数据。当某一种倾向确实清楚存在时，找到一个将数据变换成一直线（近似地）的函数。,如果选择了模型 y = a xb, 应画一个lny对 lnx 的图，看是否产生一近似直线。同样在研究模型y = a ebx 的适当性时，应画ln y 对x 的图是否产生一条直线。当认为一个确定的模型类型比较符合数据的

13、倾向时，可以用最小二乘及其他解析方法估计模型参数。最后对原始模型作拟合优度分析。如果对这一拟合不满意，可以研究其他的单项模型。如果在单项模型框架内不能满意地拟合数据，则要使用其他的技术。这就是将要讨论的多项式拟合和样条插值。,2. 高阶多项式(polynomials)模型,单项模型易于进行模型分析，然而由于其数学上的简单性，在追踪数据到的数据的趋势时，单项模型是有局限的。在某些情况下，必须考虑多于一项的模型，即多项式。当然多项式也是有其缺点的。,大家知道：有惟一的最高阶为2的多项式 y = a0+ a1x + a2x2 能够通过三个不同的点。但存在无数个阶大于2的多项式通过这三个点。,例4 带

14、式录音机的播放时间.,收集一个特定的录音机的计数器的数据及相应的录音机的播放时间。数据如下（ ci 表示计数器读数，ti 为对应的播放总时间）,可以构造一个最高阶为7的惟一多项式，通过这8个数据点。记此多项式为,因为要经过8个点，所以有下面方程组（读数缩小100倍）：,解得：,对P7(ci)的预测值舍入到小数第四位，给出了与观测到数据完全一致的预测值，产生零绝对偏差。我们真的能够认为这一模型比其他可以提供的模型更好吗？答案是否定的。后面将要谈到它的缺点。,在理论上可以证明：给定(n+1)个不同的数据点，存在惟一的一个最高阶为n的多项式通过全部数据点。这就是多项式的Lagrange形式。,高阶多

15、项式的优点与缺点：,优点：因为多项式函数易于积分和微分。如果高阶多项式 模型对研究的问题是合适的，则后面的计算就比较容易。,缺点：虽然 Lagrange 多项式能够通过给定的n个数据点，但在区间的端点附近，多项式有严重的摆动。而在作预测和外推时，端点附近的估计往往是很需要的。,高阶多项式的另一个严重的缺点是它的系数对数据的微小变化的敏感性。因此高阶多项式的这些缺点限制了它在建模中的应用。,例如，给定17个数据点如下：,高阶多项式拟合图:,x=-2:1:7; y=sin(x); p=polyfit(x,y,7); x1=-5:0.05:10; f=polyval(p,x1); plot(x,y,

16、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 时，多项式的变动非常大。,3. 低阶多项式对数据的光滑化,用低阶多项式拟合数据，不要求通过全部数据点。这将降低多项式的摆动倾向和它对数据中微小变化的敏感性。,低阶多项式拟合要求做出两个决定: 1. 决定选择所用多项式的阶； 2. 确定最佳拟合多项式的系数。,例4 带式录音机的播放时间,先用二阶多项式拟合带式录音机的播放时间，用最小二乘法决定最佳系数.,x=100 20

17、0 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 指令如下：,计算残差，得：,这些残差与时间量的阶比较，相当小。,那么怎样选择低阶多项式的阶呢？几阶多项式为合适？,注意二次多项式函数的二阶导数是常数，三阶导数为零。,对于离散的数据点

18、集，可以计算差分或均差。,计算 yi+1 - yi , i=1,2,3,4, 构成一阶差分，记为D. 从D 列的相继的一阶差分间算出二阶差分，记为D2. 继续下去，直到从 n 个数据点算出Dn-1.,假设有数据集（虚拟）：,上述虚拟数据集中二阶差分是常数，三阶差分是零. 注意各阶差分只是估计了各阶导数的分子。当 xi 是不等跨距时，若要更细致的估计，需要计算均差。,计算上述虚拟数据集的均差：,计算带式录音机数据的均差表：,从上表看到，二阶均差基本上是常数，三阶均差为零（小数点后第四位为零，不应期望准确地为零）. 可以看出数据基本上是二次的，支持用二次多项式作为经验模型。,例5. 酵母培养物的增

19、长。数据是随时间（小时）测量的酵母细胞的总数。数据与均差表如下：,其均差表能够帮助决定低阶多项式是否能够提供一个满意的经验模型。,注意一阶均差其均差表 到 t =8小时一直增长，然后开始下降。反映在 中表现为一连串的负号，指示了凹凸性的变化。这样，不能期望用仅有单个凹形的二次多项式追踪这些数据的趋势。三阶均差 的正负号就没有一连串的趋势。应该用三次多项式去拟合数据.,为了说明问题，我们也用二次的多项式拟合这组数据，最小二乘拟合结果如下： P = - 93.82 + 65.70 t - 1.12 t 2,模型拟合很差，追踪数据趋势是失败的。,二次多项式拟合(Matlab命令为polyfit)的结

20、果：,三次多项式拟合(Matlab命令为polyfit)的结果：,更好的模型应该是用逻辑斯蒂曲线拟合。,4. 三阶样条模型,三阶样条插值是非常通用的建模技术。在相邻的数据点间使用不同的三阶多项式，追踪数据的趋势。既保证基本关系的特征，同时减少摆动的倾向和对数据变化的灵敏性。,模型进一步的讨论略。,给出美国人口从1790年到1990年间的人口如表1（每10年为一个间隔），请估计出美国2010年的人口。,表1 美国人口统计数据,三、模型举例,建模方法1：,作出散点图，由图可以发现美国人口的变化规律曲线近似为一条指数函数曲线，因此我们假设美国的人口满足函数关系 x = f(t), f(t) = ea

21、+bt，a, b为待定常数 .,图1：散点图,输出 -35.4278 0.0208,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),两边取对数，得线性模型：,输入 Matlab 指令如下：,即有变换以后的最小二

22、乘估计为,还原为,检验拟合效果：很不理想.,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) %残差,输出残差平方和(很大)： Res1=24917,图2：线性化后拟合效果图,这个初值非常重要，因此前面变换数据以后最小二乘方法得到的估计常常作为现在 lsqcurvefit 的初值。,为了改进拟合效果，直接用matlab 软件的lsqcurvefit命令, 进行曲线拟合（数值方法求 的最小值点 ，其中x = f (t), f (t) = ea+bt， ），不过计算时需

23、要提供被估计参数的初值。,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; 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.0147 Res2 = 1.8766e+003,输出,得到直接进行曲线拟合的最小二乘估计为,预测美国在2010年的人口数为3

24、56.254百万。,残差平方和: Res2=1876.6，有所减少，仍比较大。,输入 Matlab 指令如下：,x2=1780:2000; y3=exp(b1(1)+b1(2).*x2); plot(x,y,+,x2,y3),检验拟合效果：有所改进.,图3 lsqcurvefit拟合效果图,建模方法2：,上述模型对过去的统计数据吻合得较好，但也存在问题，即人口是呈指数规律无止境地增长，此时人口的自然增长率随人口的增长而增长，这不可能。一般说来，当人口较少时增长得越来越快，即增长率在变大；人口增长到一定数量以后，增长就会慢下来，即增长率变小。这是因为，自然资源、环境条件等因素不允许人口无限制地增

25、长，它们对人口的增长起着阻滞作用，而且随着人口的增加，阻滞作用越来越大。而且人口最终会饱和，趋于某一个常数x，我们假设人口的净增长率为 r x(t)(1-x(t)/x)，即人口的静增长率随着人口的增长而不断减小，当 t时，净增长率趋于零。,按照这个假设，得到,这便是丹麦数学家 Verhulst（1804-1849）于19世纪中叶提出的阻滞增长模型（Logistic模型）。,满足 x(1790)=3.9 的微分方程(1)的解为,(1),(2),fun = inline(b(1)./(1+(10/39*b(1)-1)*exp(1790-x)*b(2), b, x) b0= 200 0.1;% (猜

26、测) b1 Res3=lsqcurvefit(fun,b0,x,y),b1 = 308.5182 0.0280 Res2 = 813.8362,输出,残差平方和为 Res3=813.8362，已经比指数模型时减少很多。,参数估计为 r =0.0280, x=308.5128，2010年美国的人口预计为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

27、); plot(x,y,+,x3,y4),检验拟合效果：效果比前面两种情形都好.,回顾 logistic 曲线的来源，它是下述微分方程的解：,(1),其中第二个方程是初始条件： x(1790)=3.9 由此确定了logistic曲线中有两个参数 r 和 x .,(2),显然 x(1790)=3.9 这个条件对预测2010年的美国人口并不是非常重要的。因此考虑没有初始条件的下述微分方程,(3),记,则方程 (3) 为,分离变量得,(4),因此方程 (3) 的解可表示为 (4)，其中C是积分常数，可以由初始条件确定。我们不急于确定它，而是看作一个参数。这样我们的曲线 (4) 有三个参数.,(4),

28、把 (4) 的右端看成数据变换的结果 y：,则 y 与 t 的关系是线性的。如果能预先猜测美国的极限人口 M，则可以用线性模型的最小二乘估计方法得到净增长率 r 的估计和 C 的估计。,顺便说明 logistic 一词的来源于 log it. 这充分说明了在19世纪（计算机技术远远没有现在这么发展的时代），对数变换是多么的重要！,如果猜测美国的极限人口为310（百万），则对1790年至1990年美国人口数据变换以后的最小二乘估计为（不需要初值）,它的等价形式是,方程 (4) 解 的等价形式是,(5),fun1=inline(b(1)./ (1 + exp(b(2) +b(3)*x), b, x

29、) b0=310 53.3 -0.027; b2 Res4=lsqcurvefit(fun1,b0,x,y),b2 = 387.5591 44.6944 -0.0227 Res4 = 358.1678,输出为,残差平方和为 Res4 = 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

30、)*x3); plot(x,y,+,x3,y5),检验拟合效果：效果比前面几种情形都好.,图5 b0=310, 53.3, -0.027时的拟合效果,作业,一、职业篮球队员的薪金,所有行业的薪金都在一次一次地上升，职业篮球也不例外。下面的表格展示出全国篮球协会 (即NBA)的运动员的平均薪金，时间从19851986赛季到20012002赛季，省略了有关19861987和19891990赛季的数据。,考虑所有数据，求一个线性函数与该组数据拟合。 作出线性函数的图形。 用所求函数计算关于19851986赛季到20012002赛季NBA的平均 薪金。这些计算结果是否合理？ 用所求函数预测200520

31、06和20492050赛季NBA的薪金。 用所求函数预测什么时候NBA的平均薪金达到100百万美元。,1.,2.,考虑所有数据，求一个二次多项式函数与该组数据拟合。 作出二次多项式函数的图形。 用所求函数计算关于19851986赛季到20012002赛季NBA的平均 薪金。这些计算结果是否合理？ 用所求函数预测20052006和20492050赛季NBA的薪金。 用所求函数预测什么时候NBA的平均薪金达到100百万美元。,3.,考虑所有数据，求一个三次多项式函数与该组数据拟合。 作出三次多项式函数的图形。 用所求函数计算关于19851986赛季到20012002赛季NBA的平均 薪金。这些计算结果是否合理？ 用所求函数预测20052006和20492050赛季NBA的薪金。 用所求函数预测什么时候NBA的平均薪金达到100百万美元。,考虑所有数据，求一个四次多项式函数与该组数据拟合。 作出四次多项式函数的图形。