第十章
曲线拟合
MATLAB可以用来求与一组数据最吻合的函数。为了达到目的，本章我们将学习一些简单的技术。
186 第十章曲线拟合
线性函数拟合
我们想到的最简单的情况是由线性函数很好描述的一个数集，也就是说如果我们所考虑的数据是以y = f(x)的形式给出，并且其中f(x)满足：
y = mx + b
要求得m和b的值，我们可以使用一个称为polyfit(x, y, n)的MATLAB函数，其中n是我们要MATLAB求出的多项式的次数，对于y = mx + b形式的方程，我们把n设为等于1，因此调用的语句将是polyfit(x, y, 1)。polyfit函数采用最小二乘法计算。让我们使用一些简单的例子看看如何使用它。
例10-1
解10-1
解这个问题并不需要高尔夫的相关知识，你所需要知道的是我们设想差点与平均成绩之间存在线性关系，而我们要求得一个方程来描述它。一旦你有了方程y = mx + b，你就可以预知你没有的x所对应的y值。首先我们把这两组数据输进MATLAB：
接着我们调用polyfit让MATLAB计算拟合数据的多项式的系数。要让MATLAB产生y = mx + b形式的一次多项式，首先我们需要确定x（独立变量）和y（因变量）分别是什么。在本例中充当x角色的独立变量是差点(Handicap)，而充当y角色的独立变量是平均成绩(Average)。既然我们想要产生一个一次多项式，我们用下面的方式调用polyfit：
下一步我们需要提取出MATLAB求得的系数，一般来说，MATLAB会以下面的方式产生多项式的系数：
p(x) =
p1xn + p2xn-1 + ... + pn-1x2 + pnx + pn+1
如果我们使用p = polyfit(x, y, n)调用polyfit，p(j)即引用第j个系数。在我们的例子中，polyfit是以p(1)*x + p(2)的形式返回方程的。因此我们提取系数的方法是：
现在我们产生一个函数绘出y = mx + b的直线。第一步我们需要建立x轴：
第十章曲线拟合 187
然后创建绘制直线的函数：
我们也沿着直线把实际数据标出来，它们将各自以独立的点显示。使用下面的命令，我们把这些数据用小圆圈表示：
图10-1 例10-1中产生的最小二乘法拟合线及数据的图象
结果如图10-1所示。不过拟合并不完美，很多数据点分布在我们产生的直线较远的地方。让我们看看本例中由拟合直线对特定的差点所预计的成绩值的情况，这可以通过执行下面的命令计算：
回头看刚才的原始数据表格，我们看到这些数据还是相当接近的。由于这些数据是使用平滑直线y = mx + b产生的，因此w中保存的数据通常称为平滑数据。
在很多情况下测量得到的数据都带有误差，而方程y = mx + b可以比原始数据更好地表达出变量之间的实际关系。
现在让我们了解一下拟合达到了什么程度。我们可以通过计算误差来评估拟合程度。第一项可以用来了解拟合程度的是残差。假设我们要把函数f(x)与数集ti进行拟合，残差平方的总和计算公式是：
N
A = ∑[f(x=1i) - ti]2
i现在用-t表示数集的均值。数集与均值偏差的平方和计算公式是：
N
S =
∑(xi - -t)2 i=1
那么r平方值就是：
r2 = 1 - A
如果r2 = 1，那么函数将与数据完美拟合，因此r2越接近1拟合就越好。在MATLAB，我们可以相当简单地实现这些公式。首先，我们计算这些数据的均值，一共有10个数据点，
188 第十章曲线拟合
因此：
我们可以使用MATLAB的内置函数计算数组中的数据的均值：
然后： (不用对Ave的使用感到不解，Ave是一个包含平均成绩的数组)。现在我们计算A，此时需要相关数据点的拟合值
——我们前面创建的数组w就是了，因此：
现在我们计算r2： 结果r平方值并不怎么接近1，因此并不是很完美的拟合，不过也不太差，因为与0相比它更接近1。
让我们尝试另一个更加有难度的例子。 例10-2
下面的表格列出了缅因州欢乐谷房子的平方英尺数与其对应的平均售价。求拟合这些数
首先我们把数据输入两个数组： 现在我们绘制数据的图象： 图象如图10-2所示。
虽然4000平方英尺的房子的价格看起来有点偏离正常，其它的数据还是基本上在一条
第十章曲线拟合 189
直线的周围的，让我们找出最拟合这些数据的一条直线。在我们尝试求出y = mx + b的过程中，房子的SQFT（平方英尺数）充当x的角色而平均售价充当y的角色。使用polyfit找出我们需要的系数，只需把数据传递给它并告知它我们在求一次的多项式。调用的语句是：
图10-2 例10-2中使用的数据
在本例中polyfit函数返回两个元素。它们是p(1) = m和p(2) = b。我们取出这些值：
下一步我们就可以绘制MATLAB得到的拟合直线，并把它和实际数据点比较。首先我们产生x轴所用的数据：
现在我们使用上面求得的系数产生y：
最后，我们调用subplot，把两条曲线在同一个图上显示出来：
190 第十章曲线拟合
与直线拟合的数据同时也在图10-3中显示出来。很多价格数据在很大程序上与直线相靠近，但最小和最大建筑面积不是，即实际上这并不是那么很好的拟合。polyfit程序确定了能最好描述这些数据的函数是：
图10-3 MATLAB产生的拟合例10-2中的数据的图象
y = (0.1032)x - 28.4909
其中x是房子的建筑面积而y是以千为单位的房子价格。这个模型将预计一个1600平方英尺的房子的售价是：
y = 0.1032×1600 - 28.4909 = 136.63
或者说大约$136,000。参考表格中的数据，可以看出这还是有点差距的。另一种我们可以用来体现拟合程度的方法是计算一个称为均方根（Root-Mean-Square，RMS）误差的数值，要用过这种方法来评估多项式与实际数据的拟合程度，我们可以使用下面的命令：
来说，预计值与实际售价的差的百分比是：
|$230 - $220|
≈ 4% 而对于一个1500平方英尺的房子：
|$142 - $126|
≈ 11% 房地产经纪人可能并不要求很高的精度，因此4%的误差是可以接受的，但11%的误差可能就没法让人接受了。现在让我们计算RMS误差，以便对这个模型有更好的理解。要计算RMS误差，我们应用下面的公式：
N1/2(Ti-Ai)2
RMS error = i=0
其中Ti是准确值，Ai是模型的近似值或预计值，而N是数据点的总个数。首先我们计算差额：
>> d = price - w; 模型的数据点总个数为N = 10：
∑
第十章曲线拟合 191
为了求出RMS误差，需要把差额进行平方。我们使用下面的命令把差额数组中的每个元素进行平方：
现在我们可以求RMS误差了：
确实是一个非常大的误差！RMS误差本来应该是小于1的。如何改善这种状况呢？我们可以尝试拟合更高阶的多项式。让我们使用一个二次多项式看看。使用下面的步骤来做：
图象如图10-4所示。从中可以看出，这个模型能更好地拟合数据，虽然它还未能达到完美。
现在我们求RMS
：
这次RMS减少了一半，有了大改善。注意检查w数组中的元素，它的预计值与真实值更加接近了。让我们计算r平方值：
192
第十章曲线拟合
图10-4 使用二次多项式之后提高了拟合程度
在二次多项式拟合这种情况下r平方的数值已经相当接近1了。所以我们可以认为这种模型是非常恰当的。
例10-3
一金属块加热到300°F
解10-3
第十章曲线拟合 193
我们在MATLAB输入数据：
让我们绘制数据的图象，如图10-5所示。
图10-5 冷却块的温度图象
现在我们调用polyfit产生三次多项式的系数： 既然多项式是三次的，我们的拟合函数具有下面的形式：
y = p1x3 + p
2x2 + p3x + p4
现在我们提取系数：
让我们为时间轴创建数据：
接着定义拟合函数：
我们还需要一组确切时间点上的温度数据：
让我们在同一个图中绘制数据和拟合曲线的图象： 结果如图10-6所示。看起来这个三次多项式能够很好地拟合。 让我们计算r2。首先我们计算数据的均值和S：
194
第十章曲线拟合
图10-6 三次多项式的拟合图象
现在我们求A：
我们求得的r2是： 看起来与数据完美拟合。然而这是一个误导，事实上描述这些数据的真实函数是：
T = 100 + 202.59e-0.23t - 2.59e-18.0t
让MATLAB以
long格式输出r平方的值：
是的，正如你从图中所看到的，拟合非常精确。让我们计算RMS。
RMS值比1还小，因此我们可以认为它是一个很好的近似。既然我们有函数y(t)，我们就能够估算出还没取得数据的任意时刻的温度。例如，在上面7小时之外的温度值。让我们计算经过更长时间的函数值。首先我们延伸时间线：
第十章曲线拟合 195 重新产生拟合多项式：
我们可以使用find命令提问与数据有关的问题。例如，温度小于80度是在什么时候，因为我们有可能会冒险拿起这个金属块：
find命令返回了满足这个要求的数组。我们可以使用这些位置引用时间数组t和温度数组y中的数据。首先我们提取这些时刻，并在行末添加单引号把它们转置成列向量：
现在我们取相对应的温度：
我们可以把数据放进一个两列的表格中，左边一列是时间，右边一列是温度： 举个例子，我们可以看到在14.9小时这个物块的预计值约为75度。现在让我们产生从10小时到15小时之间的数据的图象。首先我们求10小时这个时刻在数组中的索引。我们已经以0.1步长产生时间数据，因此我们可以搜索t>9.9的情况：
196
接着我们创建一个新的数组来包含这些时间值：
现在让我们攫取这个范围内的温度值并储存它们：
这时我们就可以绘制它们的图象了：
第十章曲线拟合
结果如图10-7所示。这个模型预言在5小时内金属块的温度将下降40多度。
图10-7 在所收集的数据之外的一段时间内的温度预计值的图象 指数函数的拟合 除了多项式之外还存在着其它类型函数拟合的可能，这里我们简要地提出一种。在某些情况下需要把数据拟合成指数函数。这种情况下的拟合是： y = b(10)mx 其中x是独立变量而y是因变量。现在我们定义关系： w = log10y z = x 然后用下面的形式进行数据拟合： w = p1z + p2 其中的系数p1和p2通过调用polyfit
产生。这种拟合可以在MATLAB使用下面的m = p1 , b = 102 p习题 考虑下面的数据，一位奥林匹克举重教练已经收集了每个举重运动员的年龄和最大举重磅数，他相信这两者之间存在着函数关系。
第十章曲线拟合 197 1. 使用MATLAB
2. 创建一个年龄数组，以便估算当前队伍中队员（要求每一岁都有）的最大举重。
3. 创建一个函数y实现一次拟合。
4. 这个模型预计17岁的队员最大举重为多少？
5. 实际的举重均值为多少？
6. 这个模型预计的举重均值为多少？
7. 为了计算r2，请先求出S和A。
8. 本模型的r2是多少？
9. 尝试求二次拟合，这个函数是什么？
10. 二次拟合的r2值是多少？
【参考答案在第242页】