MATLAB實現最小二乘法
阿新 • • 發佈:2019-02-13
最小二乘法
最小二乘法(又稱最小平方法)是一種數學優化技術。它通過最小化誤差的平方和尋找資料的最佳函式匹配。
利用最小二乘法可以簡便地求得未知的資料,並使得這些求得的資料與實際資料之間誤差的平方和為最小。最小二乘法還可用於曲線擬合。其他一些優化問題也可通過最小化能量或最大化熵用最小二乘法來表達。
線性函式模型
典型的一類函式模型是線性函式模型。最簡單的線性式是寫成矩陣形式為:
直接給出該式的引數解:
其中,為x的算術平均值,也可解得如下形式:
一般線性情況
若含有更多不相關模型變數,可如組成線性函式的形式:
即線性方程組:
通常人們將xij記作資料矩陣 A,引數bj記做引數向量b,觀測值yi記作Y,則線性方程組又可寫成:
上述方程運用最小二乘法匯出為線性平方差計算的形式為:。
最後的最優解為:。
示例
實驗得到4個數據(x, y):(1, 6)、(2, 5)、(3, 7)、(4, 10)。希望找出一條和這四個點最匹配的直線:。
最小二乘法採用的手段是儘量使得等號兩邊的方差最小:
求通過對β1,β2求偏導:
得到β1=3.5,β2=1.4。所以:最佳。
MATLAB實現
例一:小車時間與位移關係
% 小車時間(xi)和位移關係(yi)關係 x = [0 1 2 3 4 5 6 7 8 9]; y = [0 2 4 7 8 9 12 14 15 18]; %{ subplot(m,n,p) 其中前兩個引數 m,n是指將你的圖分成 m*n個柵格, 每個柵格用 p 來編號,而編號是按行(橫著)編號的,所以,當 m = 2,n = 2時編號規則為 1 | 2 ------ 3 | 4 所以subplot(2,2,[1 3]),就說明你這一個子圖佔據的是 1, 3兩個柵格, 而subplot(2,2,2)說明子圖僅佔據第2個柵格. %} subplot(1,2,1); plot(x,y,'o'); % 圖形的一些設定 xlabel('時間(秒)'); ylabel('位移(米)'); title('原始資料離散點') %{ grid on:是開啟網格 grid off:是關閉網格 而grid是切換兩種狀態,如果在grid off的狀態下,輸入grid,相當於grid on 相反,如果在grid on狀態下輸入grid 等價於grid off %} grid on %{ polyfit函式是matlab中用於進行曲線擬合的一個函式。其數學基礎是最小二乘法曲線擬合原理。 曲線擬合:已知離散點上的資料集,即已知在點集上的函式值,構造一個解析函式(其圖形為一曲線)使在原離散點上儘可能接近給定的值。 呼叫方法:polyfit(x,y,n)。用多項式求過已知點的表示式, 其中x為源資料點對應的橫座標,可為行向量、矩陣; y為源資料點對應的縱座標,可為行向量、矩陣; n為你要擬合的階數,一階直線擬合,二階拋物線擬合,並非階次越高越好,看擬合情況而定。 多項式在x處的值y可用下面程式計算:y=polyval(a,x,m) %} p = polyfit(x,y,1) % 0:0.01:9 起始為0,終點為9,步長0.01 x1 = 0:0.01:9; y1 = polyval(p,x1); x2 = 0:0.01:9; %{ MATLAB中的插值函式為interp1,其呼叫格式為: yi= interp1(x,y,xi,'method') 其中x,y為插值點,yi為在被插值點xi處的插值結果;x,y為向量, 'method'表示採用的插值方法,MATLAB提供的插值方法有幾種: 'nearest'是最鄰近插值, 'linear'線性插值; 'spline'三次樣條插值; 'pchip'立方插值.預設時表示線性插值 注意:所有的插值方法都要求x是單調的,並且xi不能夠超過x的範圍。 %} y2 = interp1(x,y,x2,'spline'); subplot(1,2,2); plot(x1,y1,'k',x2,y2,'r') xlabel('時間(秒)'); ylabel('位移(米)'); title('黑線為最小二乘法擬合,紅色為插值法擬合') grid on
例二:溫度和時間關係
%{
例如:對某日隔兩小時測一次氣溫。設時間為ti,氣溫為Ci,i = 0,2 ,4,…,24。資料如下:
表2 溫度(Ci)隨時間(ti)變化關係
-----------------------------------------------------------
ti 0 2 4 6 8 10 12 14 16 18 20 22 24
-----------------------------------------------------------
ci 15 14 14 16 20 23 28 27 26 25 22 18 16
-----------------------------------------------------------
%}
x = [0 2 4 6 8 10 12 14 16 18 20 22 24]
y = [15 14 14 16 20 23 26 27 26 25 22 18 16]
plot(x,y,'o')
grid on
%{
hold on 和hold off,是相對使用的
前者的意思是,你在當前圖的軸(座標系)中畫了一幅圖,再畫另一幅圖時,原來的圖還在,與新圖共存,都看得到
後者表達的是,你在當前圖的軸(座標系)中畫了一幅圖,此時,狀態是hold off,則再畫另一幅圖時,
原來的圖就看不到了,在軸上繪製的是新圖,原圖被替換了
%}
hold on
% 三階擬合 得到的 p = -0.0061 0.1474 -0.0246 13.7390是個多項式的係數
% 即擬合的曲線y = -0.0061*x3 + 0.1474*x2 - 0.0246*x + 13.7390 (其中x3表示x的3次方,x2同理)
p = polyfit(x,y,3)
x1 = 0:0.01:24
y1 = polyval(p,x1)
plot(x1,y1,'r')
% axis座標軸範圍設定
axis([0 24 12 28])
xlabel('溫度(度)');
ylabel('時間(點)');
title('溫度變化圖','position', [18,10])