【科学计算】数据拟合

数据拟合方法

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

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

目录

一、线性数据拟合

线性数据拟合指的是用线性函数拟合数据点,一般给定的数据点是形如\((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 \]

上一篇:【洛谷P4781】【模板】拉格朗日插值


下一篇:线性系统粗浅认识——第三次作业