最小二乘線性及平面擬合原理及C++實現

一、線性最小二乘擬合

使用一個簡單函數在整體上逼近已知函數，使其在整體上盡可能與原始數據曲線近似。記為：

 

 稱之為擬合曲線，若該函數為插值多項式，則所有偏差為零。

但實際情況中，我們不可能要求近似曲線

y =

 嚴格通過這么多數據點。但為了使其盡可能反映所給數據的變化趨勢，我們可以要求偏差的絕對值盡可能小，甚至是所有偏差中的最大值盡可能小。我們可以通過使選取的近似曲線在節點x_i 處的偏差的平方和達到最小來實現這一目標，這一原則就是 最小二乘原則。

按最小原則選擇的擬合曲線就稱為最小二乘擬合曲線，此方法稱為最小二乘法。

實用公式推導：

 

 

 假設我們此處有這樣一組數據點，這些點的分布接近於在一條直線上，因此選一條直線（一條曲線則加入對應方程推導）來擬合這組數據，令：

根據最小二乘原則，有：

 

 

 

 令a0,a1為未知數，則此處轉換為求二元函數S(a0,a1)的極小點問題：

 

 由此可得：

 

 

 

 

 聯立解得：

 

 

 

 

 

即得到了待求的擬合直線段。

C++實現：

bool gFittingLine(double *xArray, double *yArray, int firstIndex, int lastIndex,
				  double &a, double &b)
{
　　 int count = lastIndex-firstIndex+1;
    if(count < 2) return false;
    double s0 = (double)count, s1 = 0, s2 = 0, t0 = 0, t1 = 0;
    for(int i=firstIndex;i<=lastIndex;i++) 
	{
		s1 += xArray[i];
		s2 += (xArray[i]*xArray[i]);
		t0 += yArray[i];
		t1 += (xArray[i]*yArray[i]);
    }
    double d = s0*s2-s1*s1;
    b = (s2*t0-s1*t1)/d;
    a = (s0*t1-s1*t0)/d;
    return true;
}

　　實現對二維平面離散點的曲線擬合

二、最小二乘面擬合

對空間中的一系列散點，尋求一個近似平面，與線性最小二乘一樣，只是變換了擬合方程：

使用平面的一般方程：

Ax + By + CZ + D = 0

可以令平面方程為：

 

 

 由最小二乘法知：

 

 

 同樣分別取 a0,a1,a2的偏導數：

 

 

 即是：

 

 

 換算為矩陣形式有：

 可以直接通過克拉默法則求出a0,a1,a2的行列式表達式，有：

 

 

 

 

 c++實現(gDaterm3() 為自定義的三階行列式計算函數)：

bool gFittingPlane(double *x, double *y, double *z, int n, double &a, double &b, double &c)
{
	int i;
	double x1, x2, y1, y2, z1, xz, yz, xy, r;

	x1 = x2 = y1 = y2 = z1 = xz = yz = xy = 0;
	for(i=0; i<n; i++)
	{
		x1 += x[i];
		x2 += x[i]*x[i];
		xz += x[i]*z[i];

		y1 += y[i];
		y2 += y[i]*y[i];
		yz += y[i]*z[i];

		z1 += z[i]; 
		xy += x[i]*y[i];
	}

	r = gDeterm3(x2, xy, x1, xy, y2, y1, x1, y1, n);
	if(IS_ZERO(r)) return false;

	a = gDeterm3(xz, xy, x1, yz, y2, y1, z1, y1, n) / r;
	b = gDeterm3(x2, xz, x1, xy, yz, y1, x1, z1, n) / r;
	c = gDeterm3(x2, xy, xz, xy, y2, yz, x1, y1, z1) / r;

	return true;
}