最小二乘法擬合直線c++程式碼
最近公司的一個專案需要計算TVDI(Temperature Vegetation Dryness Index ,溫度植被幹旱指數) ,TVDI的計算公式如下(具體原理自行百度):
其中,為任意像元的地表溫度;為某一NDVI對應的最小地表溫度,對應的是溼邊;為某一NDVI對應的最大地表溫度,對應的是幹邊;a,b為溼邊的擬合方程係數,c,d為幹邊的擬合方程係數。
在擬合幹邊和溼邊的過程中,需要利用最小二乘方法來對有效的NDVI和Lst資料來進行線性擬合。因此,本文記錄在工作中用C++實現的最小二乘擬合直線,關鍵是理解最小二乘擬合直線的基本原理,實現起來比較簡單。具體的最小二乘原理再此不做過多的闡述,網上有大量的介紹資料,這裡只給出形如
- /*************************************************************************
- 最小二乘法擬合直線,y = a*x + b; n組資料; r-相關係數[-1,1],fabs(r)->1,說明x,y之間線性關係好,fabs(r)->0,x,y之間無線性關係,擬合無意義
- a = (n*C - B*D) / (n*A - B*B)
- b = (A*D - B*C) / (n*A - B*B)
- r = E / F
- 其中:
- A = sum(Xi * Xi)
- B = sum(Xi)
- C = sum(Xi * Yi)
- D = sum(Yi)
- E = sum((Xi - Xmean)*(Yi - Ymean))
- F = sqrt(sum((Xi - Xmean)*(Xi - Xmean))) * sqrt(sum((Yi - Ymean)*(Yi - Ymean)))
- **************************************************************************/
- void LineFitLeastSquares(float *data_x, float *data_y, int data_n, vector<float> &vResult)
- {
- float A = 0.0;
- float B = 0.0;
- float C = 0.0;
- float D = 0.0;
- float E = 0.0;
- float F = 0.0;
- for (int i=0; i<data_n; i++)
- {
- A += data_x[i] * data_x[i];
- B += data_x[i];
- C += data_x[i] * data_y[i];
- D += data_y[i];
- }
- // 計算斜率a和截距b
- float a, b, temp = 0;
- if( temp = (data_n*A - B*B) )// 判斷分母不為0
- {
- a = (data_n*C - B*D) / temp;
- b = (A*D - B*C) / temp;
- }
- else
- {
- a = 1;
- b = 0;
- }
- // 計算相關係數r
- float Xmean, Ymean;
- Xmean = B / data_n;
- Ymean = D / data_n;
- float tempSumXX = 0.0, tempSumYY = 0.0;
- for (int i=0; i<data_n; i++)
- {
- tempSumXX += (data_x[i] - Xmean) * (data_x[i] - Xmean);
- tempSumYY += (data_y[i] - Ymean) * (data_y[i] - Ymean);
- E += (data_x[i] - Xmean) * (data_y[i] - Ymean);
- }
- F = sqrt(tempSumXX) * sqrt(tempSumYY);
- float r;
- r = E / F;
- vResult.push_back(a);
- vResult.push_back(b);
- vResult.push_back(r*r);
- }
為了驗證該演算法的有效性,給出如下測試資料,資料來源為某論文的實驗資料:
[cpp] view plain copy- float pY[25] = { 10.98, 11.13, 12.51, 8.40, 9.27,
- 8.73, 6.36, 8.50, 7.82, 9.14,
- 8.24, 12.19, 11.88, 9.57, 10.94,
- 9.58, 10.09, 8.11, 6.83, 8.88,
- 7.68, 8.47, 8.86, 10.38, 11.08 };
- float pX[25] = { 35.3, 29.7, 30.8, 58.8, 61.4,
- 71.3, 74.4, 76.6, 70.7, 57.5,
- 46.4, 28.9, 28.1, 39.1, 46.8,
- 48.5, 59.3, 70.0, 70.0, 74.5,
- 72.1, 58.1, 44.6, 33.4, 28.6 };
該資料在Excel的擬合結果為,其中。
轉載地址 http://blog.csdn.net/pl20140910/article/details/51926886
在工程實踐中,經常遇到類似的問題:
我們做了n次實驗,獲得了一組資料
然後,我們希望知道x和y之間的函式關係。所以我們將其描繪在XOY直角座標系下,得到下面這麼一張點雲圖:
然後,我們發現,x和y「可能」是線性的關係,因為我們可以用一條直線大致的將所有的樣本點串連起來,如下圖:
所以,我們可以「猜測」。接下來的問題,就是求出a和b的值。
這看起來是一個很簡單的問題,a和b是兩個未知數,我們只需要隨意找出兩個樣本點,列出方程組:
兩個未知數,兩個方程,就可以求解出a和b的值。
然而,在這裡是不對的,或者說是不準確的。為什麼呢?因為 這個函式關係,是我們「猜測」的,並不一定是客觀正確的(雖然也許是正確的)。所以我們不能這麼簡單粗暴的方程組求解。
那怎麼辦呢?既然是「猜測」的,那麼就存在誤差。那麼我們將這個函式關係稍加修正為:
這裡, 分別是第i次實驗的因變數、自變數、誤差。
既然是「猜測」,那我們當然希望猜得準一點。那怎麼衡量準確呢?自然和e有關係。
上式變型後可得:
在這裡,a和b才是自變數,e是函式值。
這裡是最容易搞糊塗的地方,為什麼a,b是自變數,而不是x,y?
這就要提及「曲線擬合」的概念。所謂「擬合」就是說我們要找到一個函式,來「匹配」我們在實驗中獲得的樣本值。放到上面的例子,就是我們要調整a和b的值,來使得這個函式和實驗中獲得的資料更加「匹配」。所以,a和b才是「曲線擬合」過程中的自變數。
接下來,繼續如何讓誤差更小的問題。
「最小二乘法」的思想核心,就是定義一個損失函式:
顯然,如果我們調整a和b,使得Q達到最小,那麼「曲線擬合」的誤差也會最小。
這裡,Q是a,b的函式。根據高等數學的只是,Q的最小值點必然是其導數為0的點。
所以,我們令:
求解上述方程組以後得到關於a,b的一個二元二次方程組,因此可以解得a和b的值。這就是最小二乘法的整個過程。
最後說明:
(1)最小二乘法英文名Least Squares,其實翻譯成「最小平方法」,更容易讓人理解。其核心就是定義了損失函式;
(2)評價誤差的方法不止一個,還有諸如 等(當然這就不是最小二乘法了);
(3)最小二乘法不僅可以用於一次函式的擬合,還可以用於更高次函式的擬合;
(4)最小二乘法既是一種曲線擬合的方法,也可用於最優化。