曲線(xiàn)擬合與回歸分析

《曲線(xiàn)擬合與回歸分析》由會(huì)員分享，可在線(xiàn)閱讀，更多相關(guān)《曲線(xiàn)擬合與回歸分析（48頁(yè)珍藏版）》請(qǐng)?jiān)谘b配圖網(wǎng)上搜索。

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

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

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

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

5、，此差異可寫(xiě)成 此即為前述的總平方誤差 MATLAB 提供一個(gè)簡(jiǎn)單方便的左除（）指令，來(lái)解出最佳的 ，使得總平方誤差為最小。,MATLAB的最小平方解,,,利用左除來(lái)算出最佳的參數(shù)值，並同時(shí)畫(huà)出具有最小平方誤差的二次曲線(xiàn) 範(fàn)例10-2： census01.m,曲線(xiàn)擬合運(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í)際人口數(shù), 預(yù)測(cè)人口數(shù)); xlabel(年度); ylabel(美國(guó)人口總數(shù));,曲線(xiàn)擬合結(jié)果,由上述範(fàn)例，我們可以找出最佳的 因此具有最小平方誤差的拋物線(xiàn)可以寫(xiě)成：,,,,提示：左除及右除,左除的概念，可記憶如下：原先的方程式是 A*theta = y，我們可將 A移項(xiàng)至等號(hào)右邊，而得到 theta = Ay。必須小心的是：原先 A 在乘式的第一項(xiàng)，所以移到等號(hào)右邊後，A 仍然必須是除式的第一項(xiàng)。 若我們要解的方程式是 theta*A = y，則同樣的概念可得到最小平方解 theta = y/A。,根據(jù)上拋物線(xiàn)數(shù)學(xué)模

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

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

9、年的預(yù)測(cè)值 = %g （百萬(wàn)人）n, polyval(theta, 2000)); fprintf(2010年的預(yù)測(cè)值 = %g （百萬(wàn)人）n, polyval(theta, 2010));, 在2000年的預(yù)測(cè)值 = 274.622 （百萬(wàn)人） 在2010年的預(yù)測(cè)值 = 301.824 （百萬(wàn)人）,使用polyfit xx, yy, zz = peaks(pointNum); zz = zz + randn(size(zz));% 加入雜訊 surf(xx, yy, zz); axis tight,現(xiàn)在我們要用已知的基底函數(shù)，來(lái)找出最佳的 、 和 範(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ù)，可以輸入較密的點(diǎn)，得到迴歸後的曲面 範(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 值 % 畫(huà)出預(yù)測(cè)的曲面

12、 pointNum = 31; xx, yy = meshgrid(linspace(-3, 3, pointNum), linspace(-3, 3, pointNum));,在上圖中，我們猜對(duì)了基底函數(shù)，因此得到非常好的曲面擬合。 只要基底函數(shù)正確，而且雜訊是正規(guī)分佈，那麼當(dāng)資料點(diǎn)越來(lái)越多，上述的最小平方法就可以逼近參數(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,非線(xiàn)性迴歸（Nonlinear Regression）是一個(gè)比較困難的問(wèn)題，原因如下： 無(wú)法一次找到最佳解。 無(wú)法保證能夠找到最佳解。 須引用各種非線(xiàn)性最佳化的方法。 各種相關(guān)數(shù)學(xué)性質(zhì)並不明顯。 以數(shù)學(xué)來(lái)描述，假設(shè)所用的數(shù)學(xué)模型是 其中 是輸入向量， 是可變非線(xiàn)性函數(shù)，y 是輸出變數(shù)。 總平方誤差為,非線(xiàn)性迴歸,,,,,用一般最佳化（Optimization）的方法，來(lái)找出 的最小值，例如 梯度下降法（Gradient Descent） Simplex下坡式搜尋

14、（Simplex Downhill search）：此方法即為fminsearch指令所採(cǎi)用的方法 假設(shè)所用的數(shù)學(xué)模型為 其中 、 為線(xiàn)性參數(shù)，但1、2 為非線(xiàn)性參數(shù) 總平方誤差可表示： 欲找出使 為最小的 、 、1及2，需將 E 寫(xiě)成一函式，並由其它最佳化的方法來(lái)求出此函式的最小值。,非線(xiàn)性迴歸：誤差值的最小化,,,,,,,我們使用errorMeasure1.m來(lái)計(jì)算誤差值 範(fàn)例10-8： errorMeasure01.m 其中 theta 是參數(shù)向量，包含了 、 、1 及 2，data 則是觀察到的資料點(diǎn)，傳回的值 則是總平方誤差。,使用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(計(jì)算時(shí)間 = %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));,上圖的曲線(xiàn)為 fminsearch 指令產(chǎn)生的迴歸曲線(xiàn)。 fminsearch 指令是一個(gè)使用 Simplex 下坡式搜

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

18、，最佳化會(huì)更有效率,上述範(fàn)例的改良,,使用上述混成（Hybrid）的方法，函式 errorMeasure1 須改寫(xiě)成 errorMeasure2 範(fàn)例10-10： errorMeasure2.m lambda 是非線(xiàn)性參數(shù)向量， data 仍是觀察到的資料點(diǎn)，a 是利用最小平方法算出的最佳線(xiàn)性參數(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(計(jì)算時(shí)間 = %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)生較低的誤差平方和，同時(shí)所需的計(jì)算時(shí)間也比較短。,混成法範(fàn)例（3/3）,計(jì)算時(shí)間 = 0.02 誤差平方和 = 1.477226e-001,我們亦可利用變形法（Transformation），將一數(shù)學(xué)模型轉(zhuǎn)換成只包含線(xiàn)性參數(shù)的模型。 假設(shè)一模型為： 取自然對(duì)數(shù)，可得：

21、ln(a) 及 b 成為線(xiàn)性參數(shù)，我們即可用最小平方法找出其最佳值 範(fàn)例10-12： transformFit01.m,變形法,,,load data2.txt x = data2(:, 1);% 已知資料點(diǎn)的 x 座標(biāo) y = data2(:, 2);% 已知資料點(diǎn)的 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í)得到之參數(shù) b = theta(2)% 辨識(shí)得到之參數(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));,第一個(gè)小圖是ln(y)對(duì)x的作圖

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

24、y); a = exp(theta(1))% 辨識(shí)得到之參數(shù) b = theta(2)% 辨識(shí)得到之參數(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 來(lái)對(duì)誤差平方和進(jìn)行最小化，因此得到的誤差平方和，比只用變形法還要小。,以下是一些變形法可用的非線(xiàn)性模型，以及相關(guān)的轉(zhuǎn)換方法：,可使用變形法的數(shù)學(xué)模型（1/3）,,,,,,,,,,可使用變形法的數(shù)學(xué)模型（2/3）,,,,,,,,,,,,,,,可使用變形法的數(shù)學(xué)模型（3/3）,,,,,,,,,,曲線(xiàn)擬合包含下列幾個(gè)步驟： 觀察資料，並剔除明顯不合理的資料（這些資

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

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