假設對觀測數據進行擬合，得到的擬合曲線為。將觀測數據代入，得到，其和的偏差定義為

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?（1）

評價擬合結果好壞的函數稱為指標函數

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? （2）

擬合函數在觀測數據上總的偏差越小，說明擬合的越好，因此可以具體寫成

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?（3）

其中，是加權系數，表示每個觀測數據的重要程度。一般情況下，我們認為每個觀測數據是等重要的。因此，上式可以簡化為

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?（4）

其中偏差可以定義為

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? （5）

當p=0時，

當p=1時，

當p=2時，

。。。

評價函數一般寫成p范數的p次方

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? （6）

求最優擬合函數的過程是在和構成的空間上尋優的過程

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? （7）

對應的解是

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? （8）

指標函數是一個以函數為自變量的函數。至此，這個問題是變分問題（至于說如何使用變分來求解，暫且壓一壓，以后再細說）。

?

假設擬合函數為。，，。

如何求出呢？，從p=2開始推導。

評價函數對求導（因為在的取值范圍是連續的，且向量的是一個凸函數，至少是不凹的），

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? （9）

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? （10）

當時，，因此極值點為極小值點。

等價于

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? （11）

上式兩邊進行轉置，得

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?（12）

X為列滿秩矩陣，的逆矩陣存在。因此

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? （13）

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?（14）

當p=0時，。該式表示，擬合直線穿過觀測數據點越多越好。該方法對數據的分布以及其中摻雜的噪聲比較敏感，解不穩定。該方法可能不能擬合出真實數據的曲線，而是擬合了噪聲數據。

當p=1時，。評價函數導數不連續，求解不像二范數求解那么方便，因此很少使用，但在有些情況下，1范數擬合的曲線更好。

下面用一段代碼，簡單說明一下p范數對擬合結果的影響

?

lens=50; b=50;x=1:lens; y=2*x+b*randn(1,lens);x=[1 2 3 4 5]; y=[1 2 2 3 5];for k=1:200for b=1:200yk=(k-100)/10*x+(b-100)/10;d=y-yk;sign=ones(1,lens);sign(find(abs(d)<0.5))=0;err0(k,b)=sum(sign);err1(k,b)=sum(abs(d));err2(k,b)=sum(d.*d);end end figure(1) mesh(err0) figure(2) mesh(err1) figure(3) mesh(err2)figure(4) plot(x,y,'ro') hold on[a1 a2]=min(err0); [b1 b2]=min(a1); ka0=(a2(b2)-100)/10; ba0=(b2-100)/10; yk0=ka0*x+ba0; plot(x,yk0,'md')[a1 a2]=min(err1); [b1 b2]=min(a1); ka1=(a2(b2)-100)/10; ba1=(b2-100)/10; yk1=ka1*x+ba1; plot(x,yk1,'g+')[a1 a2]=min(err2); [b1 b2]=min(a1); ka2=(a2(b2)-100)/10; ba2=(b2-100)/10; yk2=ka2*x+ba2; plot(x,yk2,'bs')e10=sum(abs(y-yk0)) e11=sum(abs(y-yk1)) e12=sum(abs(y-yk2))e20=sum((y-yk0).*(y-yk0)) e21=sum((y-yk1).*(y-yk1)) e22=sum((y-yk2).*(y-yk2))legend('data','L0','L1','L2') hold off

這段代碼首先生成一組數據（x，y），然后分別使用0、1和2范數進行擬合求解。搜索范圍k=[-10,10]，b=[-10,10]。在搜索空間中找到擬合誤差最小的最小p乘解。

?

?

?

?

從評價函數的2范數出發，的解的形式中包含著復雜的矩陣運算關系，這其中應該蘊含著什么。讓我們首先從一個簡單的例子入手吧。

有a，b兩個向量，b在a上的投影p可以寫成投影長度x和a方向上單位向量的乘積

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? （15）

往a方向上做投影的投影變換矩陣為

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?（16）

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?（17）

向量b到向量a的距離等于e=b-p的模，垂直方向上的投影矩陣為

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? （18）

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?（19）

在數據擬合中，數據的個數遠遠多于未知數（待求解參數）的個數，因此這個方程不能得到精確解。我們需要找到距離Y最近的一個空間，將Y投影到該空間中。X的列向量構成的空間叫做列空間。該空間內的任意向量Xv都是X的列向量的線性組合得到，v就是組合系數。Y的近似解形如

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?（20）

近似解與Y之間的向量為

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? （21）

該向量垂直于X的列空間中的任意向量Xv，有

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?（22）

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?（23）

因為v為任意向量，因此

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?（24）

（20）式兩邊同時左乘矩陣X的轉置，將（24）式代入得到

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?（25）

解出

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? （26）

近似解

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?（27）

Y在投影變換矩陣的作用下，投影到X的列空間上，投影后得到的向量為，誤差向量，誤差的大小為。

投影矩陣從另外一個角度解釋了最小二乘法，同時也是最小p（p>2）乘法的解釋。我們發現，最小二乘法給出的解是近似解。誤差的來源方方面面，比如系統偏差，觀測誤差，記錄誤差等等。這些因素之間是一個什么樣的關系，線性的還是非線性的，一時半會兒說不清楚，我們偷個懶，將擬合問題重新形式化如下

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?（28）

我們要求的問題等價于

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? （29）

搞工程的搞來搞去不經意的發現，當隨機不可知的因素很多，獨立隨機試驗的次數很大時，由這些隨機因素造成的隨機誤差服從高斯分布。學術界的人按耐不住了，不能讓搞工程的壓下去，整出了一個中心極限定理，各種分布都會漸進服從高斯分布。

我們俗氣一把，還是從高斯分布入手，假設誤差 服從零均值高斯分布，即。為什么是零均值，因為我們不想估計出個沒用的有偏的分布出來。這種零均值誤差還有個很帶感的名字，白噪聲。對其進行傅里葉變換后，各個頻率都有響應，也就是說這種噪聲是由不同頻率的噪聲合成的。我們常見的由各個頻率合成的事物就是白色的光。所以，按照國際慣例，這種噪聲叫做白噪聲。

另外一個假設是獨立同分布，也就是每次實驗都是獨立的，但是服從的分布相同。

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?（30）

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? （31）

令上式對的導數為零，又見。也就是說，最小二乘法的概率解釋是，擬合誤差服從零均值高斯分布，擬合直線通過均值點。

假設噪聲服從形如，那么其解對應的就是最小p范數解。這種分布貌似就是指數族分布，當p=1時叫做拉普拉斯分布，p=2時就是大名鼎鼎的高斯分布。

通過上面的分析，我們可以察覺到，噪聲的模型影響了最終的擬合結果。如果噪聲是形如1范數的，那么用2范數的最小二乘法擬合出來的直線就存在偏差。如果噪聲中有粗大誤差，那么如果不能事先去除，擬合的結果很有可能擬合了噪聲而沒能擬合真實數據。從這個角度看，其實我們能夠解決的問題還比較有限，因為我們對于噪聲的認識還不夠。

從p范數的角度出發，均值是誤差形如2范數的解，中位數是誤差形如1范數的解，眾數是誤差形如0范數的解。猜想，p階矩就是誤差形如p范數的解。

到這里，我喪心病狂的把凸優化，線性代數和概率論貌似完美的在數據擬合的框架下聯系在了一起。從其反方向看，這些理論本來就是為了解決數據擬合而被提出來的。只是數學教學時，人為的把本來應該在一起的拆開了。

圖中的數據在y=2x上加入了標準差為10的零均值高斯白噪聲。使用上面推導的公式解出擬合曲線為

?

擬合誤差的樣本均值為0.abc%#%$^#e-13，樣本標準差為9.44。同時，該擬合曲線通過樣本的均值點 。

?

?

a=4; b=10; lens=100;x=1:lens; y=zeros(1,lens); y=2*x+b*randn(1,lens);X=[x;y]';Y0=zeros(lens,2); Y1=zeros(lens,2); Y2=zeros(lens,2); Y3=zeros(lens,2); pY=zeros(lens+10,2);A=zeros(lens,4); V=zeros(lens,4); D=zeros(lens,4); pY1=zeros(lens+10,2);figure for i=1:lensY0(i,:)=X(i,:);if i>3 && a>=4Y3(i,:)=Y0(i,:)-3*Y0(i-1,:)+3*Y0(i-2,:)-Y0(i-3,:);endif i>2 && a>=3Y2(i,:)=Y0(i,:)-2*Y0(i-1,:)+Y0(i-2,:);endif i>1 && a>=2Y1(i,:)=Y0(i,:)-Y0(i-1,:);endpY(i+1,:)=Y0(i,:)+Y1(i,:)+0.5*Y2(i,:)+1/6*Y3(i,:);if i<lenst=X(i+1,:)'*pinv(X(i,:)');pY1(i+2,:)=(t*X(i+1,:)')';[v d]=eig(t(1:2,1:2));A(i,:)=reshape(t,1,4);V(i,:)=reshape(v,1,4);D(i,:)=reshape(d,1,4);end endX1=[ones(1,lens);x]'; beta=inv(X1'*X1)*X1'*y'; pY2=X1*beta;plot(x,y,'go') hold on plot(pY(a:lens,1),pY(a:lens,2),'r+'); plot(pY1(a:lens,1),pY1(a:lens,2),'bs'); plot(x,pY2,'mv');

上面的代碼中pY是基于泰勒級數展開的近似估計，pY1是基于局部最小二乘估計，pY2是全局最小二乘估計。

說了這么多，我們原地踏步在一階線性估計上，高階怎么求解。改改X和就可以了。求解方法還是原來的配方，還是原來的味道。這里就不在多啰嗦了。

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? （32）

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? （33）

給大家再拜上一記大殺器

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? （34）

小伙伴們不要害怕，看一副賣萌的樣子就知道這貨是傳說中的核函數。加權函數一塊來的時候，更加兇殘的公式也就來了

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? （35）

從數據出發，我們總能對它們進行合適的解釋，發現合適的模型，擬合出合適的曲線，并給出這種解釋的好壞程度。到目前為止，我們對于數據擬合能夠給出的最簡潔的表達是

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?（36）

貌似我們可以滿足的洗洗睡了，但是一個關鍵的問題是X的具體形式能不能由數據自己說出來，也就是說從無到有可不可能。此外，還有兩個問題沒有解決，如何在保證上式成立的同時將參數的空間最小化，也就說能夠用線性擬合的不用二階多項式擬合。第二個問題就是迭代擬合。現實情況中，我們接收到的數據往往是以時間序列的形式呈現的，那t時刻和t-1時刻的模型出現不一致時，怎么處理。

未完待續

?

總結

以上是生活随笔為你收集整理的透过表象看本质！？之二数据拟合的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。