1. 程式人生 > >最小二乘法擬合直線c++程式碼

最小二乘法擬合直線c++程式碼

最近公司的一個專案需要計算TVDI(Temperature Vegetation Dryness Index ,溫度植被幹旱指數) ,TVDI的計算公式如下(具體原理自行百度):

其中,為任意像元的地表溫度;為某一NDVI對應的最小地表溫度,對應的是溼邊為某一NDVI對應的最大地表溫度,對應的是幹邊;a,b為溼邊的擬合方程係數,c,d為幹邊的擬合方程係數。

在擬合幹邊和溼邊的過程中,需要利用最小二乘方法來對有效的NDVI和Lst資料來進行線性擬合。因此,本文記錄在工作中用C++實現的最小二乘擬合直線,關鍵是理解最小二乘擬合直線的基本原理,實現起來比較簡單。具體的最小二乘原理再此不做過多的闡述,網上有大量的介紹資料,這裡只給出形如

的線性迴歸計算a,b係數以及r^2的最終計算公式,相關程式碼如下:

[cpp]  view plain  copy
  1. /************************************************************************* 
  2.  最小二乘法擬合直線,y = a*x + b; n組資料; r-相關係數[-1,1],fabs(r)->1,說明x,y之間線性關係好,fabs(r)->0,x,y之間無線性關係,擬合無意義
     
  3.  a = (n*C - B*D) / (n*A - B*B) 
  4.  b = (A*D - B*C) / (n*A - B*B) 
  5.  r = E / F 
  6.  其中: 
  7.  A = sum(Xi * Xi) 
  8.  B = sum(Xi)
     
  9.  C = sum(Xi * Yi) 
  10.  D = sum(Yi) 
  11.  E = sum((Xi - Xmean)*(Yi - Ymean)) 
  12.  F = sqrt(sum((Xi - Xmean)*(Xi - Xmean))) * sqrt(sum((Yi - Ymean)*(Yi - Ymean))) 
  13.  
  14. **************************************************************************/  
  15. void LineFitLeastSquares(float *data_x, float *data_y, int data_n, vector<float> &vResult)  
  16. {  
  17.     float A = 0.0;  
  18.     float B = 0.0;  
  19.     float C = 0.0;  
  20.     float D = 0.0;  
  21.     float E = 0.0;  
  22.     float F = 0.0;  
  23.   
  24.     for (int i=0; i<data_n; i++)  
  25.     {  
  26.         A += data_x[i] * data_x[i];  
  27.         B += data_x[i];  
  28.         C += data_x[i] * data_y[i];  
  29.         D += data_y[i];  
  30.     }  
  31.   
  32.     // 計算斜率a和截距b  
  33.     float a, b, temp = 0;  
  34.     if( temp = (data_n*A - B*B) )// 判斷分母不為0  
  35.     {  
  36.         a = (data_n*C - B*D) / temp;  
  37.         b = (A*D - B*C) / temp;  
  38.     }  
  39.     else  
  40.     {  
  41.         a = 1;  
  42.         b = 0;  
  43.     }  
  44.   
  45.     // 計算相關係數r  
  46.     float Xmean, Ymean;  
  47.     Xmean = B / data_n;  
  48.     Ymean = D / data_n;  
  49.   
  50.     float tempSumXX = 0.0, tempSumYY = 0.0;  
  51.     for (int i=0; i<data_n; i++)  
  52.     {  
  53.         tempSumXX += (data_x[i] - Xmean) * (data_x[i] - Xmean);  
  54.         tempSumYY += (data_y[i] - Ymean) * (data_y[i] - Ymean);  
  55.         E += (data_x[i] - Xmean) * (data_y[i] - Ymean);  
  56.     }  
  57.     F = sqrt(tempSumXX) * sqrt(tempSumYY);  
  58.   
  59.     float r;  
  60.     r = E / F;  
  61.   
  62.     vResult.push_back(a);  
  63.     vResult.push_back(b);  
  64.     vResult.push_back(r*r);  
  65. }  

為了驗證該演算法的有效性,給出如下測試資料,資料來源為某論文的實驗資料:

[cpp]  view plain  copy
  1. float pY[25] = { 10.98, 11.13, 12.51, 8.40, 9.27,  
  2.                       8.73, 6.36, 8.50, 7.82, 9.14,  
  3.                       8.24, 12.19, 11.88, 9.57, 10.94,  
  4.                       9.58, 10.09, 8.11, 6.83, 8.88,  
  5.                       7.68, 8.47, 8.86, 10.38, 11.08 };  
  6.   
  7. float pX[25] = { 35.3, 29.7, 30.8, 58.8, 61.4,  
  8.                       71.3, 74.4, 76.6, 70.7, 57.5,  
  9.                       46.4, 28.9, 28.1, 39.1, 46.8,  
  10.                       48.5, 59.3, 70.0, 70.0, 74.5,  
  11.                       72.1, 58.1, 44.6, 33.4, 28.6 };  

該資料在Excel的擬合結果為,其中


轉載地址   http://blog.csdn.net/pl20140910/article/details/51926886


在工程實踐中,經常遇到類似的問題:

我們做了n次實驗,獲得了一組資料

然後,我們希望知道xy之間的函式關係。所以我們將其描繪在XOY直角座標系下,得到下面這麼一張點雲圖:


然後,我們發現,xy「可能」是線性的關係,因為我們可以用一條直線大致的將所有的樣本點串連起來,如下圖:

所以,我們可以「猜測」。接下來的問題,就是求出ab的值。


這看起來是一個很簡單的問題,ab是兩個未知數,我們只需要隨意找出兩個樣本點,列出方程組:


兩個未知數,兩個方程,就可以求解出ab的值。


然而,在這裡是不對的,或者說是不準確的。為什麼呢?因為  這個函式關係,是我們「猜測」的,並不一定是客觀正確的(雖然也許是正確的)。所以我們不能這麼簡單粗暴的方程組求解。

那怎麼辦呢?既然是「猜測」的,那麼就存在誤差。那麼我們將這個函式關係稍加修正為:

這裡,  分別是第i次實驗的因變數、自變數、誤差。

既然是「猜測」,那我們當然希望猜得準一點。那怎麼衡量準確呢?自然和e有關係。

上式變型後可得:

在這裡,ab才是自變數,e是函式值。

這裡是最容易搞糊塗的地方,為什麼a,b是自變數,而不是x,y

這就要提及「曲線擬合」的概念。所謂「擬合」就是說我們要找到一個函式,來「匹配」我們在實驗中獲得的樣本值。放到上面的例子,就是我們要調整ab的值,來使得這個函式和實驗中獲得的資料更加「匹配」。所以,ab才是「曲線擬合」過程中的自變數。


接下來,繼續如何讓誤差更小的問題。

「最小二乘法」的思想核心,就是定義一個損失函式:


顯然,如果我們調整ab,使得Q達到最小,那麼「曲線擬合」的誤差也會最小。

這裡,Qa,b的函式。根據高等數學的只是,Q的最小值點必然是其導數為0的點。

所以,我們令:



求解上述方程組以後得到關於a,b的一個二元二次方程組,因此可以解得ab的值。這就是最小二乘法的整個過程。

 

最後說明:

(1)最小二乘法英文名Least Squares,其實翻譯成「最小平方法」,更容易讓人理解。其核心就是定義了損失函式

(2)評價誤差的方法不止一個,還有諸如  等(當然這就不是最小二乘法了);

(3)最小二乘法不僅可以用於一次函式的擬合,還可以用於更高次函式的擬合;

(4)最小二乘法既是一種曲線擬合的方法,也可用於最優化。