曲線擬合與回歸分析

《曲線擬合與回歸分析》由會員分享，可在線閱讀，更多相關(guān)《曲線擬合與回歸分析（48頁珍藏版）》請在裝配圖網(wǎng)上搜索。

1、MATLAB 程式設(shè)計進(jìn)階篇曲線擬合與迴歸分析,張智星 (Roger Jang) jangmirlab.org http://mirlab.org/jang 清大資工系 多媒體檢索實驗室,資料擬和簡介,資料擬合（Data Fitting） 給定一組資料（含輸入及輸出），建立一個數(shù)學(xué)模型，來逼近此資料的輸入輸出特性 如果此資料包含一維輸入及輸出，則此數(shù)學(xué)模型可以表示成一條曲線，在此情況下又稱為曲線擬合（Curve Fitting） 迴歸分析（Regression Analysis） 使用統(tǒng)計的方法來進(jìn)行資料擬和，並分析每一個變數(shù)的統(tǒng)計特性，此過程稱為迴歸分析,曲線擬合簡介,曲線擬合（Curve

2、Fitting） 欲建立的數(shù)學(xué)模型是單輸入、單輸出（Single-input Single-output，簡稱SISO），其特性可用一條曲線來表示 迴歸分析之分類 若模型是線性模型，則此類問題稱為線性迴歸（Linear Regression） 若模型是非線性模型，則稱為非線性迴歸（Nonlinear Regression）。,觀察資料是美國自 1790 至 1990 年（以 10 年為一單位）的總?cè)丝?，此資料可由載入檔案 census.mat 得到： 範(fàn)例10-1： censusPlot01.m,曲線擬合：美國人口範(fàn)例,load census.mat% 載入人口資料 plot(cdate, p

3、op, o);% cdate 代表年度，pop 代表人口總數(shù) xlabel(年度); ylabel(美國人口總數(shù));,模型選取 由上圖資料分佈，可猜測這適合的曲線可能是二次拋物線 其中y為輸出，x為輸入， 則為此模型的參數(shù)。由於參數(shù)相對於y呈線性關(guān)係，所以此模型為稱線性模型 目標(biāo) 找出最好的參數(shù)值，使得模型輸出與實際資料越接近越好，此過程即稱為線性迴歸,曲線擬合之模型選取,,,,,,,,,,,曲線擬和的平方誤差 假設(shè)觀察資料可寫成 ,i= 121。當(dāng)輸入為 時，實際輸出為 。 模型的預(yù)測值為 平方誤差： 總平方誤差 是參數(shù) 的函數(shù)，這也是我們要最小化的目標(biāo)函數(shù)，可表示如下：,曲線擬合之

4、目標(biāo)函數(shù),,,,,,求得參數(shù) 、 、 的最佳值 求出 對 、 、 的導(dǎo)式，令其為零，即可解出 、 、 的最佳值。 總平方誤差 為 的二次式 導(dǎo)式 、 及 為 的一次式 令上述導(dǎo)式為零之後，我們可以得到一組三元一次線性聯(lián)立方程式，就可以解出參數(shù) 、 、 的最佳值。,目標(biāo)函數(shù)之求解,假設(shè) 21 個觀察點均通過此拋物線，將這 21 個點帶入拋物線方程式，得到下列21個等式： 亦可寫成 其中 、 為已知， 為未知向量。,矩陣表示法,,,,,,,觀察 上述21個方程式，只有 3 個未知數(shù) ，所以通常不存在一組解來滿足這 21 個方程式。 在一般情況下，只能找到一組 ，使得等號兩邊的差異為最小

5、，此差異可寫成 此即為前述的總平方誤差 MATLAB 提供一個簡單方便的左除（）指令，來解出最佳的 ，使得總平方誤差為最小。,MATLAB的最小平方解,,,利用左除來算出最佳的參數(shù)值，並同時畫出具有最小平方誤差的二次曲線 範(fàn)例10-2： census01.m,曲線擬合運(yùn)算範(fàn)例,load census.mat% 載入人口資料 plot(cdate, pop, o);% cdate 代表年度，pop 代表人口總數(shù) A = ones(size(cdate)), cdate, cdate.2; y = pop; theta = Ay;% 利用左除，找出最佳的 theta 值 plot(cdate,

6、 pop, o, cdate, A*theta, -); legend(實際人口數(shù), 預(yù)測人口數(shù)); xlabel(年度); ylabel(美國人口總數(shù));,曲線擬合結(jié)果,由上述範(fàn)例，我們可以找出最佳的 因此具有最小平方誤差的拋物線可以寫成：,,,,提示：左除及右除,左除的概念，可記憶如下：原先的方程式是 A*theta = y，我們可將 A移項至等號右邊，而得到 theta = Ay。必須小心的是：原先 A 在乘式的第一項，所以移到等號右邊後，A 仍然必須是除式的第一項。 若我們要解的方程式是 theta*A = y，則同樣的概念可得到最小平方解 theta = y/A。,根據(jù)上拋物線數(shù)學(xué)模

7、型，我們可以預(yù)測美國在 2000 年的人口總數(shù)為： 範(fàn)例10-3： census02.m,以模型預(yù)測人口總數(shù),load census.mat% 載入人口資料 A = ones(size(cdate)), cdate, cdate.2; theta = Apop;% 利用左除，找出最佳的 theta 值 t=2000; pop2000 = 1, t, t2*theta; % 在 2000 年美國人口線數(shù)預(yù)測值 t=2010; pop2010 = 1, t, t2*theta; % 在 2010 年美國人口線數(shù)預(yù)測值 fprintf(美國人口在2000年的預(yù)測值 = %g （百萬人）n, pop2

8、000); fprintf(美國人口在2010年的預(yù)測值 = %g （百萬人）n, pop2010);, 美國人口在2000年的預(yù)測值 = 274.622 （百萬人） 美國人口在2010年的預(yù)測值 = 301.824 （百萬人）,上述例子推廣，得到一個 n 次多項式： 利用多項式的數(shù)學(xué)模型來進(jìn)行曲線擬合，通稱為多項式擬合（Polynomial Fitting） 由於多項式擬和的應(yīng)用面很廣，MATLAB 提供 polyfit 指令來找出多項式最佳參數(shù) Polyval 指令來進(jìn)行多項式的求值,多項式擬和,,使用 polyfit % 進(jìn)行二次多項式擬合，找出 theta 值 fprintf(2000

9、年的預(yù)測值 = %g （百萬人）n, polyval(theta, 2000)); fprintf(2010年的預(yù)測值 = %g （百萬人）n, polyval(theta, 2010));, 在2000年的預(yù)測值 = 274.622 （百萬人） 在2010年的預(yù)測值 = 301.824 （百萬人）,使用polyfit xx, yy, zz = peaks(pointNum); zz = zz + randn(size(zz));% 加入雜訊 surf(xx, yy, zz); axis tight,現(xiàn)在我們要用已知的基底函數(shù)，來找出最佳的 、 和 範(fàn)例10-6： peaks02.m 由此找

10、出的 值和最佳值 相當(dāng)接近。,曲面擬合範(fàn)例（4/6）,pointNum = 10; xx, yy, zz = peaks(pointNum); zz = zz + randn(size(zz))/10;% 加入雜訊 x = xx(:);% 轉(zhuǎn)為行向量 y = yy(:);% 轉(zhuǎn)為行向量 z = zz(:);% 轉(zhuǎn)為行向量 A = (1-x).2.*exp(-(x.2)-(y+1).2), (x/5-x.3-y.5).*exp(-x.2-y.2), exp(-(x+1).2-y.2); theta = Az% 最佳的 theta 值,theta = 3.0088 -10.0148 -0.2

11、924,,根據(jù)上求得之參數(shù)，可以輸入較密的點，得到迴歸後的曲面 範(fàn)例10-7： peaks03.m,曲面擬合範(fàn)例（5/6）,pointNum = 10; xx, yy, zz = peaks(pointNum); zz = zz + randn(size(zz))/10;% 加入雜訊 x = xx(:); y = yy(:); z = zz(:);% 轉(zhuǎn)為行向量 A = (1-x).2.*exp(-(x.2)-(y+1).2), (x/5-x.3-y.5).*exp(-x.2-y.2), exp(-(x+1).2-y.2); theta = Az;% 最佳的 theta 值 % 畫出預(yù)測的曲面

12、 pointNum = 31; xx, yy = meshgrid(linspace(-3, 3, pointNum), linspace(-3, 3, pointNum));,在上圖中，我們猜對了基底函數(shù)，因此得到非常好的曲面擬合。 只要基底函數(shù)正確，而且雜訊是正規(guī)分佈，那麼當(dāng)資料點越來越多，上述的最小平方法就可以逼近參數(shù)的真正數(shù)值。,曲面擬合範(fàn)例（6/6）,x = xx(:); y = yy(:);% 轉(zhuǎn)為行向量 A = (1-x).2.*exp(-(x.2)-(y+1).2), (x/5-x.3-y.5).*exp(-x.2-y.2), exp(-(x+1).2-y.2); zz = r

13、eshape(A*theta, pointNum, pointNum); surf(xx, yy, zz); axis tight,非線性迴歸（Nonlinear Regression）是一個比較困難的問題，原因如下： 無法一次找到最佳解。 無法保證能夠找到最佳解。 須引用各種非線性最佳化的方法。 各種相關(guān)數(shù)學(xué)性質(zhì)並不明顯。 以數(shù)學(xué)來描述，假設(shè)所用的數(shù)學(xué)模型是 其中 是輸入向量， 是可變非線性函數(shù)，y 是輸出變數(shù)。 總平方誤差為,非線性迴歸,,,,,用一般最佳化（Optimization）的方法，來找出 的最小值，例如 梯度下降法（Gradient Descent） Simplex下坡式搜尋

14、（Simplex Downhill search）：此方法即為fminsearch指令所採用的方法 假設(shè)所用的數(shù)學(xué)模型為 其中 、 為線性參數(shù)，但1、2 為非線性參數(shù) 總平方誤差可表示： 欲找出使 為最小的 、 、1及2，需將 E 寫成一函式，並由其它最佳化的方法來求出此函式的最小值。,非線性迴歸：誤差值的最小化,,,,,,,我們使用errorMeasure1.m來計算誤差值 範(fàn)例10-8： errorMeasure01.m 其中 theta 是參數(shù)向量，包含了 、 、1 及 2，data 則是觀察到的資料點，傳回的值 則是總平方誤差。,使用fminsearch的範(fàn)例（1/3）,functi

15、on squaredError = errorMeasure1(theta, data) x = data(:,1); y = data(:,2); y2 = theta(1)*exp(theta(3)*x)+theta(2)*exp(theta(4)*x); squaredError = sum((y-y2).2);,欲求出 的最小值，我們可使用 fminsearch 指令 範(fàn)例10-9： nonlinearFit01.m,使用fminsearch的範(fàn)例（2/3）,load data.txt theta0 = 0 0 0 0; tic theta = fminsearch(errorMeas

16、ure1, theta0, , data); fprintf(計算時間 = %gn, toc); x = data(:, 1); y = data(:, 2); y2 = theta(1)*exp(theta(3)*x)+theta(2)*exp(theta(4)*x); plot(x, y, ro, x, y2, b-); legend(Sample data, Regression curve); fprintf(誤差平方和 = %dn, sum((y-y2).2));,上圖的曲線為 fminsearch 指令產(chǎn)生的迴歸曲線。 fminsearch 指令是一個使用 Simplex 下坡式搜

17、尋法（Downhill Simplex Search）的最佳化方法，用來找出 errorMeasure1 的極小值，並傳回 theta 的最佳值。,使用fminsearch的範(fàn)例（3/3）,計算時間 = 0.03 誤差平方和 = 5.337871e-001,改良方向 上述方法把所有參數(shù)全部視為非線性參數(shù)。混成法將線性與非線性參數(shù)分開，各用不同的方法來處理。 上例數(shù)學(xué)模型為 、 線性參數(shù)：最小平方法，即左除或 1、2 非線性參數(shù)： Simplex 下坡式搜尋（即 fminsearch） 混成法的優(yōu)點 最小平方法能夠在非線性參數(shù)固定的情況下，一次找到最好的線性參數(shù)的值，因為搜尋空間的維度由4降為2

18、，最佳化會更有效率,上述範(fàn)例的改良,,使用上述混成（Hybrid）的方法，函式 errorMeasure1 須改寫成 errorMeasure2 範(fàn)例10-10： errorMeasure2.m lambda 是非線性參數(shù)向量， data 仍是觀察到的資料點，a 是利用最小平方法算出的最佳線性參數(shù)向量，傳回的 squareError 仍是總平方誤差,混成法範(fàn)例（1/3）,function squaredError = errorMeasure2(lambda, data) x = data(:,1); y = data(:,2); A = exp(lambda(1)*x) exp(lambd

19、a(2)*x); a = Ay; y2 = a(1)*exp(lambda(1)*x)+a(2)*exp(lambda(2)*x); squaredError = sum((y-y2).2);,欲用此混成法求出誤差平方和的最小值 範(fàn)例10-11： nonlinearFit02.m,混成法範(fàn)例（2/3）,load data.txt lambda0 = 0 0; tic lambda = fminsearch(errorMeasure2, lambda0, , data); fprintf(計算時間 = %gn, toc); x = data(:, 1); y = data(:, 2); A =

20、exp(lambda(1)*x) exp(lambda(2)*x); a = Ay; y2 = A*a; plot(x, y, ro, x, y2, b-); legend(Sample data, Regression curve); fprintf(誤差平方和 = %dn, sum((y-y2).2));,此種混成法可以產(chǎn)生較低的誤差平方和，同時所需的計算時間也比較短。,混成法範(fàn)例（3/3）,計算時間 = 0.02 誤差平方和 = 1.477226e-001,我們亦可利用變形法（Transformation），將一數(shù)學(xué)模型轉(zhuǎn)換成只包含線性參數(shù)的模型。 假設(shè)一模型為： 取自然對數(shù)，可得：

21、ln(a) 及 b 成為線性參數(shù)，我們即可用最小平方法找出其最佳值 範(fàn)例10-12： transformFit01.m,變形法,,,load data2.txt x = data2(:, 1);% 已知資料點的 x 座標(biāo) y = data2(:, 2);% 已知資料點的 y 座標(biāo) A = ones(size(x)) x;,a = 4.3282 b =-1.8235 誤差平方和 = 8.744185e-001,變形法範(fàn)例（1/2）,theta = Alog(y); subplot(2,1,1) plot(x, log(y), o, x, A*theta); xlabel(x); ylabel(l

22、n(y)); title(ln(y) vs. x); legend(Actual value, Predicted value); a = exp(theta(1))% 辨識得到之參數(shù) b = theta(2)% 辨識得到之參數(shù) y2 = a*exp(b*x); subplot(2,1,2); plot(x, y, o, x, y2); xlabel(x); ylabel(y); legend(Actual value, Predicted value); title(y vs. x); fprintf(誤差平方和 = %dn, sum((y-y2).2));,第一個小圖是ln(y)對x的作圖

23、， 第二個小則是 y對 x的作圖。 經(jīng)由變形法之後，此最小平方法所得到的最小總平方誤差是 而不是原模型的總平方誤差： 通常E為最小值時，E不一定是最小值，但亦離最小值不遠(yuǎn)矣！,變形法範(fàn)例（2/2）,,,若要求取原始E的最小值，可再用 fminsearch ，並以變形法得到的 a 及 b 為搜尋的起點 範(fàn)例10-13： transformFit02.m,變形法之改進(jìn)範(fàn)例（1/2）,load data2.txt x = data2(:, 1);% 已知資料點的 x 座標(biāo) y = data2(:, 2);% 已知資料點的 y 座標(biāo) A = ones(size(x)) x; theta = Alog(

24、y); a = exp(theta(1))% 辨識得到之參數(shù) b = theta(2)% 辨識得到之參數(shù) theta0 = a, b;% fminsearch 的啟始參數(shù) theta = fminsearch(errorMeasure3, theta0, , data2); x = data2(:, 1); y = data2(:, 2); y2 = theta(1)*exp(theta(2)*x);,誤差平方和 = 1.680455e-001,變形法之改進(jìn)範(fàn)例（2/2）,plot(x, y, o, x, y2); xlabel(x); ylabel(y); legend(Actual val

25、ue, Predicted value); title(y vs. x); fprintf(誤差平方和 = %dn, sum((y-y2).2));,由上述範(fàn)例可以看出，我們可以先使用變形法，先找出大略的參數(shù)值，再用 fminsearch 來對誤差平方和進(jìn)行最小化，因此得到的誤差平方和，比只用變形法還要小。,以下是一些變形法可用的非線性模型，以及相關(guān)的轉(zhuǎn)換方法：,可使用變形法的數(shù)學(xué)模型（1/3）,,,,,,,,,,可使用變形法的數(shù)學(xué)模型（2/3）,,,,,,,,,,,,,,,可使用變形法的數(shù)學(xué)模型（3/3）,,,,,,,,,,曲線擬合包含下列幾個步驟： 觀察資料，並剔除明顯不合理的資料（這些資

26、料在統(tǒng)計學(xué)上稱為 Outliers）。 根據(jù)資料，選定數(shù)學(xué)模型及相關(guān)參數(shù)。 根據(jù)線性或非線性迴歸的各種方法，以及一組給定的訓(xùn)練資料（Training Data），算出參數(shù)的最佳值。 觀察模型的誤差，以及模型對於其它測試資料（Test Data）的效能，以決定此模型的適用性。若適用，則停止此演算法。反之，若不適用，則根據(jù)模型的對於訓(xùn)練及測試資料的誤差程度，重新修正模型，並回到步驟 3。,曲線擬合工具箱的使用（1/4）,以上步驟，需要經(jīng)驗，且必須反覆進(jìn)行，可能耗費(fèi)大量時間，有鑑於此，MathWorks 公司在 MATLAB 6.x 後，推出了曲線擬合工具箱（Curve Fitting Toolbo

27、x），讓使用者能以 GUI （圖形使用者介面）的方式，來進(jìn)行曲線擬合，並能快速地檢視擬合的結(jié)果和成效。,曲線擬合工具箱的使用（2/4）,說明曲線擬合工具箱的使用 首先，先載入 enso.mat，裡面包含兩個變數(shù)： month：每個資料點發(fā)生的相對月份 pressure：在復(fù)活島（Easter Island）和澳洲的達(dá)爾文（Darwin）兩地的大氣壓力差值，取其整個月的平均值。此差值會導(dǎo)引整個南半球的貿(mào)易風(fēng)（Trade Winds）流向 根據(jù)這兩個變數(shù)，就可以呼叫曲線擬合工具箱來進(jìn)行曲線的分析與擬合 範(fàn)例10-14： cftool01.m,曲線擬合工具箱的使用（3/4）,load enso.mat% 載入 month 和 pressure 變數(shù) cftool(month, pressure);% 呼叫曲線擬合工具箱,此時此工具箱會將資料點畫出來，並將相關(guān)的操作介面顯示如上圖。,曲線擬合工具箱的使用（4/4）,