【科学计算】数据拟合

数据拟合方法

数据拟合相对于插值放宽了限制，尽管都是找一个特定函数来尽可能地估计待估函数，但数据拟合不要求函数一定要经过数据节点，而是希望函数与数据节点的差异尽可能小。数据拟合被广泛运用于发掘数据之中隐含的模式，尽管有多种多样的数据拟合方法，但最常用的拟合还是基于最小二乘法的拟合。

本文将从最简单的一元线性数据拟合开始，向多元数据拟合拓展，并介绍多项式拟合。

目录

• 数据拟合方法
  • 一、线性数据拟合
    • 二元线性拟合
    • 多元线性拟合
    • 解超定方程组
  • 二、非线性二元数据的拟合
    • 非线性数据转化为线性数据拟合
    • 带参数函数拟合
    • 正交多项式拟合
    • 正交多项式举例

一、线性数据拟合

线性数据拟合指的是用线性函数拟合数据点，一般给定的数据点是形如\((x_{i1},x_{i2},\cdots,y_i)\)的，这里下标\(i\)代表数据点的序号，后面的序号代表观测变量的取值，而\(y\)则是因变量的取值。往往我们会先获得一组这样的数据点，再选择一个函数来拟合，最简单的拟合函数就是线性函数。

二元线性拟合

如再将自变量限制成只有一个，则观测数据点形如\((x_i,y_i)\)，用于拟合的线性函数就是\(y=a+bx\)。然而，同样的\(a,b\)不会让所有数据点都落在一条直线上，因此，我们使用残差来度量数据点和目标函数的偏离程度：对于给定的\(a,b\)和第\(i\)组数据点，残差是

\[r_i=y_i-(a+bx_i). \]

最小二乘法的目标是最小化残差平方和\(Q\)，也就是找到一组\(a,b\)使得\(Q\)最小。

\[\min Q=\sum_{i=1}^n[y_i-(a+bx_i)]^2. \]

为了最小化\(Q\)，应对\(a,b\)求偏导，从而得到一个可以解出\(a,b\)的方程组，这个方程组就称为正规方程组。

\[\frac{\partial Q}{\partial a}=2\sum_{i=1}^{n}[y_i-(a+bx_i)]=0,\\ \frac{\partial Q}{\partial b}=2\sum_{i=1}^nx_i[y_i-(a+bx_i)]=0. \]

从中可以解出\(a,b\)的值，从而得到线性拟合函数的形式（以下\(\bar x\)和\(\bar y\)分别代表\(\frac{1}{n}\sum_{i=1}^{n}x_i\)和\(\frac{1}{n}\sum_{i=1}^{n}y_i\)）。

\[b=\frac{\sum_{i=1}^{n}(x_i-\bar x)(y_i-\bar y)}{\sum_{i=1}^{n}(x_i-\bar x)^2},\quad a=\bar y-b\bar x. \]

事实上平方和分解式应该很熟悉，在处理方差、协方差时常用，即

\[\sum_{i=1}^{n}(x_i-\bar x)^2=\sum_{i=1}^{n}x_i^2-n\bar x^2,\quad \sum_{i=1}^{n}(x_i-\bar x)(y_i-\bar y)=\sum_{i=1}^{n}x_iy_i-n\bar x\bar y. \]

多元线性拟合

二元线性拟合由于足够简单，所以可以直接显式写出\(a,b\)的表达式，只要直接代入计算即可。不过，往往数据点会是多元的，形如\((x_{i1},x_{i2},\cdots,x_{im},y_i)\)，线性拟合需要找到一个这样的表达式：

\[y_i=a_0+a_1x_{i1}+a_2x_{i2}+\cdots+a_mx_{im}. \]

为方便讨论，先引入一些符号。以下总记\(n\)为观测数，\(m\)为自变量数，称\(A\)为系数矩阵，\(Y\)为右端向量，\(X\)为解向量，\(\varepsilon\)为残差向量。

\[A=\begin{bmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1m} \\ 1 & x_{21} & x_{22} & \cdots & x_{2m} \\ \vdots & \vdots & \vdots & & \vdots \\ 1 & x_{n1} & x_{n2} & \cdots & x_{nm} \end{bmatrix},Y=\begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix},X=\begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_m \end{bmatrix},\varepsilon=\begin{bmatrix} r_1 \\ r_2 \\ \vdots \\ r_n \end{bmatrix}, \]

于是有\(AX+\varepsilon=Y\)，我们的目标是解出一组\(X\)使得残差平方和最小，注意到这里残差平方和就是\(\varepsilon'\varepsilon\)。接下来的推导过程需要运用一些矩阵求导的公式，事实上矩阵求导不过是多重偏导的简化写法。

\[\begin{aligned} \frac{\partial \varepsilon'\varepsilon}{\partial X}&= \frac{\partial}{\partial X}(Y-AX)'(Y-AX)\\ &=\frac{\partial}{\partial X}(Y'Y-X'A'Y-Y'AX+X'A'AX)\\ &=-2A'Y+2A'AX\\ &=0 \end{aligned} \]

从中我们知道\(A'AX=A'Y\)。如果\(A'A\)可逆，就得到方程的解为

\[X=(A'A)^{-1}A'Y. \]

解超定方程组

最小二乘法的一个重要应用是解超定方程组，即约束条件个数多于未知量个数的方程组。以下是一个具体的例子：

\[2x_1-x_2=1,\\ 8x_1+4x_2=0,\\ 2x_1+x_2=1,\\ 7x_1-x_2=8,\\ 4x_1=3. \]

不论我们怎么选取\(x_1,x_2\)的值，都只能让以上方程组中的至多两个同时成立，因此我们的目标不再是给出方程的准确解，而是一个最小误差解。一般地，将其转化为下面的规划问题：

\[\min \sum r_i^2\\ \text{s.t. }\left\{\begin{array}l 2x_1 - x_2 + r_1 = 1,\\ 8x_1 + 4x_2 + r_2 = 0,\\ 2x_1 + x_2 + r_3 = 1,\\ 7x_1 - x_2 + r_4 = 8,\\ 4x_1 + r_5 = 3. \end{array} \right. \]

这个方程组形式与上面多元线性拟合的形式类似，因此解也相同的，只需类似地引入符号即可：

\[A=\begin{bmatrix} 2 & -1 \\ 8 & 4 \\ 2 & 1 \\ 7 & -1 \\ 4 & 0 \end{bmatrix},b=\begin{bmatrix} 1 \\ 0 \\ 1 \\ 8 \\ 3 \end{bmatrix},\\ X=(A'A)^{-1}A'b=\begin{bmatrix} 0.7927 \\ -1.4641 \end{bmatrix} \]

注意，原式中的方程不可再约，否则将影响残差前面的权重，即\(8x_1+4x_2=0\)是不能约化成\(2x_1+x_2=0\)求解的（否则答案将不一致）。

二、非线性二元数据的拟合

接下来我们将对二元数据拟合作更详细的讨论，因为二元数据的可视化最为容易，且呈现出的关系往往不是线性函数这么简单。

非线性数据转化为线性数据拟合

在对二元数据进行拟合之前，我们往往会先绘制散点图观察数据的形态，如果数据呈现的不是线性形态，那么用线性函数来拟合肯定是效果不佳的。

如果可以通过一定的变换，使得数据点经过变换后变得近似线性，再使用线性拟合效果就会更好。即找到一个变换\(y'=g(y),x'=f(x)\)，使得\((x',y')=(f(x),g(y))\)在坐标点内近似线性，再对以下方程进行拟合：

\[y'=a+bx'. \]

这就将非线性的拟合问题转化为了线性拟合问题。注意，虽然我们在选择变换函数\(f,g\)时往往选择非线性的变换函数（如指数、对数、倒数变换），但是当获得数据观测值后，诸\(x'\)和\(y'\)都是可计算的，故这仍然是一个线性拟合，我们一般称之为“模型对参数是线性的”。

如果不知道如何选择合适的变换函数，选择多项式无疑是最简单的做法，此时往往设\(x_i=x^i\)，这样就可以衍生出若干个变量，将一元线性拟合转化为多元线性拟合的问题，即我们假设

\[y=b_0+b_1x+b_2x^2+\cdots \]

然后用最小二乘法计算一组最合适的回归系数\((b_0,b_1,\cdots)\)。

带参数函数拟合

现在我们对最小二乘法进行推广，使其不局限于线性函数\(y=AX\)。考虑一个具有\(k\)个参数的多元函数

\[\varphi(x)=F(x;a_1,\cdots,a_k). \]

这个意思是，我们知道这个函数的形式，知道他含有多少个参数，但是不知道这些参数的具体值，也就不知道这个函数的具体表达式，即我们假定总体符合一个函数族，并根据样本对参数进行估计。依然使用最小二乘法，使残差平方和最小，也就是如下的最优化问题：

\[\min Q=\sum_{i=1}^{n}[y_i-F(x_i;a_1,\cdots,a_k)]^2. \]

分别对每一个参数求偏导数，并让偏导数为\(0\)，这个方程组称为法方程组。

\[\frac{\partial Q}{\partial a_i}=0,\quad i=1,2,\cdots,k. \]

当然，这个方程组不一定容易求解，更不一定是线性的。

正交多项式拟合

正交多项式拟合是对普通多项式拟合的优化，前述的多项式拟合当多项式系数过高时会出现病态性，即当观测数据出现一点微小扰动时，会对拟合系数造成很大的波动，这样的拟合曲线显然是不稳定的。因此，我们往往选择一系列多项式\(P_1(x),\cdots,P_k(x)\)，用这些多项式的线性组合来表示一个最终的拟合多项式，即

\[P(x)=\sum_{j=1}^{k}a_jP_j(x). \]

如果\(P_1(x),\cdots,P_k(x)\)已知，则\(P(x)\)是一个含参方程，此时再用最小二乘法，法方程组是一个线性方程组：

\[Q=\sum_{i=1}^{n}\left[y_i-\sum_{j=1}^{k}a_jP_j(x_i) \right]^2,\\ \frac{\partial Q}{\partial a_l}=\sum_{i=1}^{n}P_l(x_i)\left[y_i-\sum_{j=1}^{k}a_jP_j(x_i) \right]=0,\\ \sum_{j=1}^{k}\left(\sum_{i=1}^{n}P_j(x_i)P_l(x_i) \right)a_j=\sum_{i=1}^{n}y_iP_l(x_i). \]

注意到当\(P_j\)、\(x_i\)、\(y_i\)全部已知的情况下，只有\(a_j\)是我们需要求解的参数，因此这是一个线性方程组，可以用矩阵形式表示为

\[\begin{bmatrix} \displaystyle\sum_{i=1}^{n}P_1(x_i)P_1(x_i) & \displaystyle\sum_{i=1}^{n}P_2(x_i)P_1(x_i) & \cdots & \displaystyle\sum_{i=1}^{n}P_k(x_i)P_1(x_i) \\ \displaystyle\sum_{i=1}^{n}P_1(x_i)P_2(x_i) & \displaystyle\sum_{i=1}^{n}P_2(x_i)P_2(x_i) & \cdots & \displaystyle\sum_{i=1}^{n}P_k(x_i)P_2(x_i) \\ \vdots & \vdots & & \vdots \\ \displaystyle\sum_{i=1}^{n}P_1(x_i)P_k(x_i) & \displaystyle\sum_{i=1}^{n}P_2(x_i)P_k(x_i) & \cdots & \displaystyle\sum_{i=1}^{n}P_k(x_i)P_k(x_i) \\ \end{bmatrix}\begin{bmatrix} a_1 \\ a_2 \\ \vdots \\ a_k \end{bmatrix}=\begin{bmatrix} \displaystyle\sum_{i=1}^{n}y_iP_1(x_i) \\ \displaystyle\sum_{i=1}^{n}y_iP_2(x_i) \\ \vdots \\ \displaystyle\sum_{i=1}^{n}y_iP_k(x_i) \end{bmatrix} \]

由于符号过于繁琐，我们现在来简化这个式子。注意到这样的求和满足内积的对称性、正定性、线性性，所以我们定义存在一个函数\(f\)使得\(y_i=f(x_i)\)，并且定义两个函数之间的内积为

\[\langle f,g\rangle =\sum_{i=1}^{n}f(x_i)g(x_i), \]

这样就将函数空间转化为一个内积空间，并且将以上的矩阵方程转化为

\[\begin{bmatrix} \langle P_1,P_1\rangle & \langle P_1,P_2\rangle & \cdots & \langle P_1,P_k\rangle \\ \langle P_2,P_1\rangle & \langle P_2,P_2\rangle & \cdots & \langle P_2,P_k\rangle \\ \vdots & \vdots & & \vdots \\ \langle P_k,P_1\rangle & \langle P_k,P_2\rangle & \cdots & \langle P_k,P_k\rangle \end{bmatrix}=\begin{bmatrix} a_1 \\ a_2 \\ \vdots \\ a_k \end{bmatrix}=\begin{bmatrix} \langle f,P_1 \rangle \\ \langle f,P_2 \rangle \\ \vdots \\ \langle f,P_k \rangle \end{bmatrix}. \]

如何求解这个矩阵方程而得到每一个系数使得最小二乘法成立？我们自然能想到，如果系数矩阵是非退化的，那么Cramer法则能保证方程组有唯一的解，从而可以用一些数值计算方法求解线性方程组。但最简单的方法，无疑是将系数矩阵转化为对角矩阵，从而可以直接写出方程组的解：

\[a_j=\frac{\langle f, P_j\rangle}{\langle P_j,P_j\rangle}. \]

这相当于要求我们预先选定的多项式组\((P_1,P_2,\cdots,P_k)\)是一个正交组，即\(\forall i\ne j\)，有\(\langle P_i,P_j\rangle=0\)。

现在回过头来注意我们的推导过程，在我们的定义中，内积使用样本点的函数值求和来定义的，但在事先选择多项式的时候，我们还没有获得样本点，从而自然无法计算这些和式，也就无法求出这里精确定义的“内积”。但如果我们假设样本具有“变异性”，即样本可以在一个预定的范围\([a,b]\)内均匀取值，那么在这一区间上求和，就近似于在这一区间上的区间，所以我们将更改内积的定义为

\[\langle f,g\rangle =\int_a^b f(x)g(x)\mathrm{d}x, \]

从而使用这样的内积定义，来确定我们所需要的正交多项式。同时，我们会希望我们的拟合函数可以是任意次数的，所以这一个正交组，还应当是一个正交基，即我们要求（为了美观对称，增加了一个\(P_0\)对应常数项）\(\mathrm{span}(P_0,\cdots,P_k)=\mathcal P_k(\mathbb{R})\)。只要我们能找到这样的正交多项式，就可以将系数矩阵近似为对角矩阵，从而直接计算各个\(a_j\)的取值。

注意，这里的近似指的是，我们在定义时使用积分作为内积的定义，但是实际应用时，还应当选择求和来求内积（否则由于\(f\)的未知性，无法计算右端\(\langle f,P_j\rangle\)的值），这就使得系数矩阵并不严格是一个对角阵。不过，我们依然将其近似视为一个对角阵，然后用简单的方法计算各个系数。同时也要注意，定义时给样本\(x\)的取值拟定了一个上下界\([a,b]\)，这在实际生活时也不一定存在。

正交多项式举例

现在，我们知道了正交多项式的定义以及引入缘由，也将给出一些正交多项式的例子。

数学分析中提到的Legendre多项式是定义在\([-1,1]\)上的正交多项式组，它的定义是

\[P_n(x)=\frac{1}{2^nn!}\frac{\mathrm{d}^n}{\mathrm{d}x^n}(x^2-1)^n,\quad n=0,1,2,\cdots \]

勒让德多项式具有如下的性质：

1. 递推关系：\(P_0(x)=1\)，\(P_1(x)=x\)，

   \[(n+1)P_{n+1}(x)=(2n+1)xP_n(x)-nP_{n-1}(x). \]

2. 正交性：

   \[\int_{-1}^{1}P_n(x)P_m(x)\mathrm{d}x=\langle P_n,P_m \rangle=0,\\ \int_{-1}^{1}P_n^2(x)\mathrm{d}x=\langle P_n,P_n\rangle=\frac{2}{2n+1}. \]

3. 可变换性。Legendre多项式是在\([-1,1]\)上的正交多项式，如果要将其变为\([a,b]\)上的正交多项式，只需引入变换

   \[x=\frac{a+b}{2}+\frac{b-a}{2}t, \]

   此时如果\(t\)在\([-1,1]\)上变换，则\(x\)在\([a,b]\)上变换，故

   \[\tilde P_n(x)=P_n(t)=P_n\left(\frac{2x-(b+a)}{b-a} \right). \]

在介绍其他的正交多项式之前，先介绍权函数。所谓权函数，实际上指的是在观测数据时，不同的残差随着自变量的不同可能对问题具有不同的影响，可以直观地想想一下，如果真实数据点是\((1,2)\)和\((1000,2000)\)，那么量级为\(1\)的残差会使数据点变为\((1,3)\)、\((1000,2001)\)，显然残差在\(x=1\)时起到的影响更大，在\(x=1000\)时起到的影响更小，但是最小二乘法作用时可不管自变量的影响，同样的残差在计算残差平方和时起到相同的作用，这有时并不是我们需要的效果。

因此，更加良好的做法是根据自变量的不同，为残差赋予一定的权重，不妨设\(w_i\)是残差\(e_i\)的权重，则此时残差平方和变为\(Q=\sum w_i[y_i-P(x_i)]^2\)。确定权重的过程，一般认为权重与自变量的取值有一定关系，所以将残差的权重视为自变量的某个函数，即

\[Q=\sum_{i=1}^{n}w(x_i)[y_i-P(x_i)]^2. \]

这里，函数\(w(x)\)就被成为权函数，它依赖于自变量的取值，给出残差的一个权重。这时对内积的定义变为

\[\langle f, g\rangle=\int_{a}^b w(x)f(x)g(x)\mathrm{d}x, \]

在这样的内积定义下，正交多项式也类似定义，而Legendre多项式则可以视为权函数\(w(x)\equiv 1\)的正交多项式。下面给出一些带非常数权函数的正交多项式例。

Laguerre多项式是\([0,\infty)\)上，关于权函数\(w(x)=e^{-x}\)的正交多项式，即\(x\)越大权重越小，其定义为

\[L_n(x)=e^x\frac{\mathrm{d}^n}{\mathrm{d}x^n}(x^ne^{-x}),\quad n=0,1,2,\cdots \]

满足

\[\int_0^{\infty}e^{-x}L_n^2(x)\mathrm{d}x=(n!)^2,\\ \int_0^{\infty}e^{-x}L_n(x)L_m(x)\mathrm{d}x=0. \]

Chebyshev多项式是\([-1,1]\)上，关于权函数\(w(x)=\frac{1}{\sqrt{1-x^2}}\)的正交多项式，即\(x\)越接近\(1\)权重越大，其定义为

\[T_n(x)=\cos(n\cdot\arccos x). \]

递推式为\(T_0(x)=1\)，\(T_1(x)=x\)，

\[T_{n+1}(x)=2xT_n(x)-T_{n-1}(x),\quad n=1,2,\cdots \]