解决数据拟合问题最重要方法是最小二乘法和回归分析。如，我们需要从一组测定的数据（例如N个点（xi，yi）（i=0,1，…，m））去求得自变量 x 和因变量 y 的一个近似解表达式 y=f（x），这就是由给定的 N 个点（xi，yi）（i=0,1，…，m）求数据拟合的问题。（注意数据拟合和数据插值是不同的，举个例子：因为测量数据往往不可避免地带有测试误差，而插值多项式又通过所有的点（xi，yi），这样就使插值多项式保留了这些误差，从而影响逼近精度，使得插值效果不理想） 所以使用最小二乘法曲线拟合法：即寻求已知函数的一个逼近函数y=f（x），使得逼近函数从总体与已知函数的偏差按某种方法度量能达到最小，而又不一定通过全部的点（xi，yi），这个时候就需要使用最小二乘法曲线拟合法。 数据拟合的具体做法是：对给定的数据（xi，yi）（i=0,1，…，m），在取定的函数类 ϕ \phi ϕ中使误差 r i = p ( x i ) − y i ( i = 0 , 1 , … , m ) r_{i}=p\left(x_{i}\right)-y_{i}(i=0,1, \ldots, m) ri​=p(xi​)−yi​(i=0,1,…,m)的平方和最小，即 [ ∑ i = 0 m r i 2 = ∑ i = 0 m [ p ( x i − y i ) ] 2 ] min ⁡ \left[\sum_{i=0}^{m} r_{i}^{2}=\sum_{i=0}^{m}\left[p\left(x_{i}-y_{i}\right)\right]^{2}\right]_{\min } [i=0∑m​ri2​=i=0∑m​[p(xi​−yi​)]2]min​ 从几何意义讲，即寻求与给定点 x i − y i ( i = 0 , 1 , … , m ) x_{i}-y_{i}(i=0,1, \ldots, m) xi​−yi​(i=0,1,…,m) 的距离平方和为最小的曲线y=p（x）。函数p（x）称为拟合函数或最小二乘解，求拟合函数p（x）的方法称为曲线拟合的最小二乘法。 在曲线拟合中，函数类 ϕ \phi ϕ 可有不同的选取方法。 MATLAB工具箱中提供了最小二乘拟合函数 polyfit() -->多项式曲线拟合 具体调用格式有三种：

（1）P = polyfit(X,Y,N) 返回次数为 n 的多项式 p(x) 的系数，该阶数是 y 中数据的最佳拟合（在最小二乘方式中）。p 中的系数按降幂排列，p 的长度为 n+1. 其中 X 为输入的向量x，Y 为得到的函数值，N 表示拟合的最高次数，返回的P值为拟合的多项式：

P ( 1 ) × X n + P ( 2 ) × X n − 1 + … + P ( N ) × X + P ( N + 1 ) P(1) \times X^{n}+P(2) \times X^{n-1}+\ldots+P(N) \times X+P(N+1) P(1)×Xn+P(2)×Xn−1+…+P(N)×X+P(N+1)

（2） [P,S] = polyfit(X,Y,N)还返回一个结构体 S，S可用作 polyval 的输入来获取误差估计值（S为由范德蒙矩阵的QR分解的R分量）。 其中 X 为输入的向量x，Y 为得到的函数值，N 表示拟合的最高次数，返回的P值为拟合的多项式：

P ( 1 ) × X n + P ( 2 ) × X n − 1 + … + P ( N ) × X + P ( N + 1 ) P(1) \times X^{n}+P(2) \times X^{n-1}+\ldots+P(N) \times X+P(N+1) P(1)×Xn+P(2)×Xn−1+…+P(N)×X+P(N+1)

（3）[P,S,MU] = polyfit(X,Y,N)还返回 mu，mu是一个二元素向量，包含中心化值和缩放值。mu(1) 是 mean(x)，mu(2) 是 std(x)。使用这些值时，polyfit 将 x 的中心置于零值处并缩放为具有单位标准差：

x ^ = x − x ˉ σ x \hat{x}=\frac{x-\bar{x}}{\sigma_{x}} x^=σx​x−xˉ​

这种中心化和缩放变换可同时改善多项式和拟合算法的数值属性。

下面来看一些具体的例子：（来源于帮助文档，改编）

使用 polyfit 将一个 7 次多项式与这些点拟合。

在更精细的网格上计算多项式并绘制结果图。

运行结果如下：

将 4 次多项式与 5 个点拟合。通常，对于 n 个点，可以拟合 n-1 次多项式以便完全通过这些点。

在由 0 和 2 之间的点组成的更精细网格上计算原始函数和多项式拟合。

在更大的区间 [0,2] 中绘制函数值和多项式拟合，其中包含用于获取以圆形突出显示的多项式拟合的点。多项式拟合在原始 [0,1] 区间中的效果较好，但在该区间外部很快与拟合函数出现差异。

运行结果如下：

确定6次逼近多项式的系数。

p =

为了查看拟合情况如何，在各数据点处计算多项式，并生成说明数据、拟合和误差的一个表。

T =

48×4 table

在该区间中，插值与实际值非常符合。创建一个绘图，以显示在该区间以外，外插值与实际数据值如何快速偏离。

T =

11×2 table

使用带三个输入的 polyfit 拟合一个使用中心化和缩放的 5 次多项式，这将改善问题的数值属性。polyfit 将 year 中的数据以 0 为进行中心化，并缩放为具有标准差 1，这可避免在拟合计算中出现病态的 Vandermonde （范德蒙）矩阵。

使用带四个输入的 polyval，根据缩放后的年份 (year-mu(1))/mu(2) 计算 p。绘制结果对原始年份的图。

计算在 x 中的点处拟合的多项式 p。用这些数据绘制得到的线性回归模型。

计算以 p 为系数的一次多项式在 x 中各点处的拟合值。将误差估计结构体指定为第三个输入，以便 polyval 计算标准误差的估计值。标准误差估计值在 delta 中返回。

绘制原始数据、线性拟合和 95% 预测区间 y ± 2 Δ y \pm 2 \Delta y±2Δ。

最小二乘拟合多项式系数，以向量的形式返回。p 的长度为 n+1，包含按降幂排列的多项式系数，最高幂为 n。如果 x 或 y 包含 NaN 值且 n < length(x)，则 p 的所有元素均为 NaN。

使用 polyval 计算 p 在查询点处的解。

如果 y 中的数据是随机的，则 p 的估计协方差矩阵是 (Rinv*Rinv’)*normr^2/df，其中 Rinv 是 R 的逆矩阵。

如果 y 中数据的误差呈独立正态分布，并具有常量方差，则 [y,delta] = polyval(…) 可生成至少包含 50% 的预测值的误差边界。即 y ± delta 至少包含 50% 对 x 处的未来观测值的预测值。

中心化值和缩放值，以二元素向量形式返回。mu(1) 为 mean(x)，mu(2) 为 std(x)。这些值以单位标准差将 x 中的查询点的中心置于零值处。

使用 mu 作为 polyval 的第四个输入以计算 p 在缩放点 (x - mu(1))/mu(2) 处的解。

局限性

算法： polyfit 使用 x 构造具有 n+1 列和 m = length(x) 行的 Vandermonde 矩阵 V 并生成线性方程组 ( x 1 n x 1 n − 1 … 1 n 2 x 2 − 1 … 1 ⋮ ⋮ ⋱ ⋮ n m x m − 1 … 1 ) ( p 1 p 2 ⋮ p n + 1 ) = ( y 1 y 2 ⋮ y m ) \left(\begin{array}{cccc} x_{1}^{n} & x_{1}^{n-1} & \ldots & 1 \\ n_{2} & x_{2}-1 & \ldots & 1 \\ \vdots & \vdots & \ddots & \vdots \\ n_{m} & x_{m}-1 & \ldots & 1 \end{array}\right)\left(\begin{array}{c} p_{1} \\ p_{2} \\ \vdots \\ p_{n+1} \end{array}\right)=\left(\begin{array}{c} y_{1} \\ y_{2} \\ \vdots \\ y_{m} \end{array}\right) ⎝⎜⎜⎜⎛​x1n​n2​⋮nm​​x1n−1​x2​−1⋮xm​−1​……⋱…​11⋮1​⎠⎟⎟⎟⎞​⎝⎜⎜⎜⎛​p1​p2​⋮pn+1​​⎠⎟⎟⎟⎞​=⎝⎜⎜⎜⎛​y1​y2​⋮ym​​⎠⎟⎟⎟⎞​

其中 polyfit 使用 p = V\y 求解。由于 Vandermonde 矩阵中的列是向量 x 的幂，因此条件数 V 对于高阶拟合来说通常较大，生成一个奇异系数矩阵。在这些情况下，中心化和缩放可改善系统的数值属性以产生更可靠的拟合。