【科學計算】數據擬合

數據擬合方法

數據擬合相對於插值放寬了限制，盡管都是找一個特定函數來盡可能地估計待估函數，但數據擬合不要求函數一定要經過數據節點，而是希望函數與數據節點的差異盡可能小。數據擬合被廣泛運用於發掘數據之中隱含的模式，盡管有多種多樣的數據擬合方法，但最常用的擬合還是基於最小二乘法的擬合。

本文將從最簡單的一元線性數據擬合開始，向多元數據擬合拓展，並介紹多項式擬合。

目錄

• 數據擬合方法
  • 一、線性數據擬合
    • 二元線性擬合
    • 多元線性擬合
    • 解超定方程組
  • 二、非線性二元數據的擬合
    • 非線性數據轉化為線性數據擬合
    • 帶參數函數擬合
    • 正交多項式擬合
    • 正交多項式舉例

一、線性數據擬合

線性數據擬合指的是用線性函數擬合數據點，一般給定的數據點是形如\((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 \]