1. 程式人生 > >數字訊號處理:曲線擬合演算法-----最小二乘法

數字訊號處理:曲線擬合演算法-----最小二乘法

	在迴歸分析中,一般任意的資料都可以用一條曲線來表示,這個曲線可以用某一個高次方的代數多項式 y= a + b*x + c*(x)2 + ...來描述,其中 a , b , ...是常數。但是這樣通過每個點的曲線是沒意義的,也不能表示y和x的真實的相關關係。趨勢變化才是兩者相關關係的合理解釋。在已知y=f(x)的形式時,運用最小二乘法可以計算出它的引數。
	現在我們通過求直線迴歸方程來了解最小二乘法。
	假定 y1 為預測值,它和 x之間的關係是 y1 = a + b*x,y2 為實測值。根據誤差理論,算術平均值是最佳值,由此算出的均方誤差最小。我們定義 y3 = y2 - y1為實測值相對於迴歸直線的誤差,那麼滿足均方誤差最小的迴歸直線方程,就是最佳配置直線。
	由上文中的滿足均方誤差最小的迴歸直線方程,可以轉化為:
			∑(y2-b*x-a)2 =  ∑ (y3)2 = min(表示最小值)
	要使等式成立,根據微分學的知識,必須使等式左邊對a和b的偏微分分別為零,然後微分得到:
		   ∑(y2-b*x-a)*x = 0
		   ∑(y2-b*x-a) = 0
	分解開來:
			∑ (x*y2) - b*∑(x)2 -a*∑ x = 0
		  ∑ y2 - b*∑ x - Na =0
	其中N為資料個數。
	解上面的方程組為:
	     b = (N*∑(x*y2)-(∑ y2)*(∑ x))/(N*∑(x)2 - (∑ x)2)
	     a = ((∑ (x)2)*(∑ y2)-(∑ (x*y2))*(∑ x))/(N*∑ (x)2-(∑ x)2)
   我們定義 y4 和 x4 為所有 y2 和x的平均值,那麼就有
        a = y4 - b*x4
        b = ∑((y2 - y4)*(x - x4))/∑(x - x4)2
  
  舉例:採集一段時間靜止下的加速度資料,然後求這個擬合曲線。

y = load(‘F:\matlab\impact\s.txt’); %匯入加速度資料,這y是一個一列多行的矩陣 x = 1:length(y); %獲取x的資料序列,x是一個一行多列的矩陣 n = length(x); %獲取資料長度 x = x’; %因為後面要用到矩陣運算,所以x 和y的矩陣維度必須一致,而這個程式碼就是把x的列轉成行 x_ave = mean(x); %求x的平均值 y_ave = mean(y); %求y的平均值 b = sum((y-y_ave).(x-x_ave))/sum((x-x_ave).^2); %b的計算 a = y_ave-bx_ave; %a的計算 r = sum((y-y_ave).(x-x_ave))/sqrt(sum((x-x_ave).2).*sum((y-y_ave).

2)); %相關係數計算 yhat = a+bx; %y的估計值 lyy = sum((y-y_ave).^2); %總平方和 U = sum((yhat-y).^2); %迴歸平方和 Q = lyy-U; %剩餘平方和 Sy = sqrt(lyy/(n-1)); %總均方和 Sq = sqrt(Q); %迴歸均方差 Syhat = sqrt(U/(n-2)); %剩餘均方差 Sa = Syhatsqrt(1/n+x_ave2/sum((x-x_ave).2)); %a的均方差 Sb = Syhat/sqrt(sum((x-x_ave).^2)); %b的均方差 xx = min(x):1:max(x); % yy = a+b
xx; %迴歸直線,擬合直線 deltayy = Syhat*sqrt(1+1/n+(xx-x_ave).2./sum((x-x_ave).2)); %迴歸直線誤差 y1 = yy+deltayy; %誤差上限 y2 = yy-deltayy; %誤差下限 plot(xx,yy,’-’,x,y,‘s’); xlabel(‘X’);ylabel(‘Y’); 程式執行得到a = -400.9312 b = 0.0022 那麼擬合直線為 y = 0.0022 *x -400.9312 a的值很小,可以認為y的值與x的值沒有太大關係。這也符合加速度計靜止的事實。 圖中藍色直線就是擬合直線