曲線擬合與回歸分析ppt課件

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

2、輸入、單輸出Single-input Single-output，簡稱SISO，其特性可用一條曲線來表示 n迴歸分析之分類n假設(shè)模型是線性模型，則此類問題稱為線性迴歸Linear Regression n假設(shè)模型是非線性模型，則稱為非線性迴歸Nonlinear Regression。n觀察資料是美國自 1790 至 1990 年以 10 年為一單位的總?cè)丝冢速Y料可由載入檔案 census.mat 得到：n範(fàn)例10-1： censusPlot01.m 曲線擬合：美國人口範(fàn)例 load census.mat% 載入人口資料plot(cdate, pop, o);% cdate 代表年度，pop

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

4、總平方誤差 是參數(shù) 的函數(shù)，這也是我們要最小化的目標(biāo)函數(shù)，可表示如下： 曲線擬合之目標(biāo)函數(shù) ),(iiyxixiyiy2210210),;(iiixaxaaaaaxf2)(iixfy 211222102112210)(),(iiiiiiixaxaayxfyaaaE21, 0,aaaEn求得參數(shù) 、 、 的最正確值n求出 對 、 、 的導(dǎo)式，令其為零，即可解出 、 、 的最正確值。n總平方誤差 為 的二次式n導(dǎo)式 、 及 為 的一次式n令上述導(dǎo)式為零之後，我們可以得到一組三元一次線性聯(lián)立方程式，就可以解出參數(shù) 、 、 的最正確值。目標(biāo)函數(shù)之求解 0a1a2a0aE1aE2aEE0a1a2aE0a

5、1a2a0a1a2a21, 0,aaa21, 0,aaan假設(shè) 21 個觀察點均通過此拋物線，將這 21 個點帶入拋物線方程式，得到以下21個等式：n亦可寫成n其中 、 為知， 為未知向量。矩陣表示法 212212211022222101212110yxaxaayxaxaayxaxaa yAyyyxxxxxx212132122121222211111Ayn觀察n上述21個方程式，只需 3 個未知數(shù) ，所以通常不存在一組解來滿足這 21 個方程式。n在普通情況下，只能找到一組 ，使得等號兩邊的差異為最小，此差異可寫成 n n 此即為前述的總平方誤差nMATLAB 提供一個簡單方便的左除指令，來解

6、出最正確的 ，使得總平方誤差為最小。MATLAB的最小平方解T321,)()(2AyAyAyTEn利用左除來算出最正確的參數(shù)值，並同時畫出具有最小平方誤差的二次曲線n範(fàn)例10-2： census01.m曲線擬合運算範(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, pop, o, cdate, A*theta, -);legend(實

7、際人口數(shù), 預(yù)測人口數(shù));xlabel(年度);ylabel(美國人口總數(shù));曲線擬合結(jié)果 n由上述範(fàn)例，我們可以找出最正確的 n因此具有最小平方誤差的拋物線可以寫成： 00654. 0,51.23,21130,210aaa2221000654. 051.2321130)(xxxaxaaxfy175018001850190019502000050100150200250年 度美國人口總數(shù) 實 際 人 口 數(shù)預(yù) 測 人 口 數(shù)提示：左除及右除n左除的概念，可記憶如下：原先的方程式是 A*theta = y，我們可將 A移項至等號右邊，而得到 theta = Ay。必須小心的是：原先 A 在乘式的

8、第一項，所以移到等號右邊後，A 依然必須是除式的第一項。n假設(shè)我們要解的方程式是 theta*A = y，則同樣的概念可得到最小平方解 theta = y/A。n根據(jù)上拋物線數(shù)學(xué)模型，我們可以預(yù)測美國在 2000 年的人口總數(shù)為：n範(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=2021; pop

9、2021 = 1, t, t2*theta;% 在 2021 年美國人口線數(shù)預(yù)測值fprintf(美國人口在2000年的預(yù)測值 = %g 百萬人n, pop2000);fprintf(美國人口在2021年的預(yù)測值 = %g 百萬人n, pop2021); 美國人口在2000年的預(yù)測值 = 274.622 百萬人 美國人口在2021年的預(yù)測值 = 301.824 百萬人n上述例子推廣，得到一個 n 次多項式：n利用多項式的數(shù)學(xué)模型來進行曲線擬合，通稱為多項式擬合Polynomial Fittingn由於多項式擬和的應(yīng)用面很廣，MATLAB 提供npolyfit 指令來找出多項式最正確參數(shù)nPol

10、yval 指令來進行多項式的求值多項式擬和 nnxaxaaxfy10)(n運用 polyfit & polyval 使程式碼更加簡潔n範(fàn)例10-4： census03.mnpolyfit(cdate, pop, 2)中的 2 代表用到的模型是 2 次多項式load census.mat% 載入人口資料theta = polyfit(cdate, pop, 2);% 進行二次多項式擬合，找出 theta 值fprintf(2000年的預(yù)測值 = %g 百萬人n, polyval(theta, 2000);fprintf(2021年的預(yù)測值 = %g 百萬人n, polyval(theta, 20

11、21); 在2000年的預(yù)測值 = 274.622 百萬人 在2021年的預(yù)測值 = 301.824 百萬人運用polyfit & polyval的範(fàn)例 nMATLAB 下輸入census，可對 census 資料進行曲線擬合的結(jié)果，如下：n上述圖形可以看出，當(dāng)多項式的次數(shù)越來越高時，外插常會出現(xiàn)不可信的結(jié)果。n這表示選用的模型參數(shù)太高，雖然誤差的平方和變小了，但是預(yù)測的可靠度也下降了。模型複雜度對預(yù)測的影響 n線性迴歸的勝利與否，與所選取的模型息息相關(guān)n模型所含的參數(shù)越多，平方誤差會越小n假設(shè)參數(shù)個數(shù)等於資料點個數(shù)，平方誤差會等於零，但這並不表示預(yù)測會最準(zhǔn)，因為資料點含有雜訊n完全吻合資料的

12、模型亦代表此模型受雜訊的影響最大，預(yù)測之準(zhǔn)確度也會較差n如何選取模型，是線性迴歸的一個艱苦課題！nLeave-one-out method線性迴歸的模型選取n多輸入、單輸出的線性迴歸數(shù)學(xué)模型寫成n其中 為輸入，y 為輸出， 、 、 為此模型的參數(shù)， ，則是知的函數(shù)，稱為基底函數(shù)Basis Functionsn所給的資料點為 ，稱為取樣資料Sample Data或訓(xùn)練資料Training Data。多輸入之線性迴歸1a2ananifi1),(xmiyii1),(x)()()()(2211xxxxnnfafafafyxyn將上述資料點帶入模型後可得：n或可表示成矩陣格式：矩陣表示法 )()()()

13、()()()()(221122111mmmm1111xxxxxxxxnnmnnfafafafyfafafafyymnAmnmnyyaaffff111111)()()()(x xx xx xx xn由於 (即資料點個數(shù)遠大於可變參數(shù)個數(shù))，欲使上式成立，須加上一誤差向量 e：n平方誤差則可寫成n求 的最正確值n直接取 對 的偏微分，並令其等於零，即可得到一組 n 元一次的線性聯(lián)立方程式n用矩陣運算來表示， 的最正確值可表示成 n也可以運用 MATLAB 的左除來算出 的最正確值，即 。平方誤差和的最小化 nmyeA)()()(2AyAyeeeETT)(EyAAATT1yA n理論上，最正確的 值

14、為 ，但是 容易呵斥電腦內(nèi)部計算的誤差，MATLAB 實際在計算左除時，會按照矩陣 A 的特性而選用最正確的方法，因此可以得到較穩(wěn)定且正確的數(shù)值解。n欲知如何推導(dǎo)最正確的 值，可參考筆者另一著作：NeuroFuzzy and Soft Computing，Prentice Hall，1997。提示yAAATT11AATn在 MATLAB 下輸入 peaks，可以畫出一個凹凸有致的曲面，如下：n此函數(shù)的方程式如下：曲面擬合範(fàn)例1/6 222222)1(53)1(231510)1 (3yxyxyxeeyxxexzn在以下說明中，假設(shè)：n數(shù)學(xué)模型的基底函數(shù)知n訓(xùn)練資料包含正規(guī)分佈的雜訊n上述函數(shù)可寫

15、成：n其中我們假設(shè) 、 和 是未知參數(shù)，n 則是平均為零、變異為 1 的正規(guī)分佈雜訊。曲面擬合範(fàn)例2/6nyxfyxfyxfneyxfyxfneeyxxexzyxyxyxyx),(),(),(31),(10),(331510)1 (3332211)1(21)1(53)1(222222222123n假設(shè)要獲得 100 筆訓(xùn)練資料n範(fàn)例10-5： peaks01.mnrandn 指令的運用即在參與正規(guī)分佈雜訊。上圖為我們搜集到的訓(xùn)練資料，由於雜訊很大，所以和原先未帶雜訊的圖形差異很大。曲面擬合範(fàn)例3/6 pointNum = 10;xx, yy, zz = peaks(pointNum);zz =

16、 zz + randn(size(zz); % 參與雜訊 surf(xx, yy, zz);axis tight-3-2-10123-202-4-20246n現(xiàn)在我們要用知的基底函數(shù)，來找出最正確的 、 和 n範(fàn)例10-6： peaks02.mn由此找出的 值和最正確值 相當(dāng)接近。曲面擬合範(fàn)例4/6 123pointNum = 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(-

17、(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.2924 31,10, 3n根據(jù)上求得之參數(shù)，可以輸入較密的點，得到迴歸後的曲面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)為行向量

18、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ù)測的曲面pointNum = 31; xx, yy = meshgrid(linspace(-3, 3, pointNum), linspace(-3, 3, pointNum); n在上圖中，我們猜對了基底函數(shù)，因此得到非常好的曲面擬合。n只需基底函數(shù)正確，而且雜訊是正規(guī)分佈，那麼當(dāng)資料點越來越多，上述的最小平方法就可以逼近參數(shù)的真正數(shù)值。曲面擬合範(fàn)例6/6 x = x

19、x(:); 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 = reshape(A*theta, pointNum, pointNum);surf(xx, yy, zz);axis tight -3-2-10123-202-4-20246n非線性迴歸Nonlinear Regression是一個比較困難的問題，緣由如下：n無法一次找到最正確解。n無法保證能夠找到最正確解。n須援用各種非線性最正確化的方法。n各種相關(guān)數(shù)學(xué)性質(zhì)並不明顯。n以數(shù)

20、學(xué)來描畫，假設(shè)所用的數(shù)學(xué)模型是n其中 是輸入向量， 是可變非線性函數(shù)，y 是輸出變數(shù)。n總平方誤差為非線性迴歸 ),(xfy x21),()(miiixfyEn用普通最正確化Optimization的方法，來找出 的最小值，例如n梯度下降法Gradient DescentnSimplex下坡式搜尋Simplex Downhill search：此方法即為fminsearch指令所採用的方法n假設(shè)所用的數(shù)學(xué)模型為n其中 、 為線性參數(shù)，但1、2 為非線性參數(shù)n總平方誤差可表示：n欲找出使 為最小的 、 、1及2，需將 E 寫成一函式，並由其它最正確化的方法來求出此函式的最小值。非線性迴歸：誤差值

21、的最小化)(Exxeaeay21211a2amixxieaeayaaEi12212121221),(E1a2an我們運用errorMeasure1.m來計算誤差值n範(fàn)例10-8： errorMeasure01.mn其中 theta 是參數(shù)向量，包含了 、 、1 及 2，data 則是觀察到的資料點，傳回的值 則是總平方誤差。運用fminsearch的範(fàn)例1/3 function squaredError = errorMeasure1(theta, data)x = data(:,1);y = data(:,2);y2 = theta(1)*exp(theta(3)*x)+theta(2)*e

22、xp(theta(4)*x);squaredError = sum(y-y2).2); 1a2aEn欲求出 的最小值，我們可運用 fminsearch 指令n範(fàn)例10-9： nonlinearFit01.m運用fminsearch的範(fàn)例2/3 Eload data.txttheta0 = 0 0 0 0;tictheta = fminsearch(errorMeasure1, theta0, , data);fprintf(計算時間 = %gn, toc);x = data(:, 1);y = data(:, 2);y2 = theta(1)*exp(theta(3)*x)+theta(2)*

23、exp(theta(4)*x); plot(x, y, ro, x, y2, b-); legend(Sample data, Regression curve);fprintf(誤差平方和 = %dn, sum(y-y2).2);n上圖的曲線為 fminsearch 指令產(chǎn)生的迴歸曲線。nfminsearch 指令是一個運用 Simplex 下坡式搜尋法Downhill Simplex Search的最正確化方法，用來找出 errorMeasure1 的極小值，並傳回 theta 的最正確值。運用fminsearch的範(fàn)例3/3計算時間 = 0.03誤差平方和 = 5.337871e-001

24、00.81.820123456 Sample dataRegression curven改良方向n上述方法把一切參數(shù)全部視為非線性參數(shù)?；斐煞▽⒕€性與非線性參數(shù)分開，各用不同的方法來處理。n上例數(shù)學(xué)模型為n 、 線性參數(shù)：最小平方法，即左除或n1、2 非線性參數(shù)： Simplex 下坡式搜尋即 fminsearchn混成法的優(yōu)點n最小平方法能夠在非線性參數(shù)固定的情況下，一次找到最好的線性參數(shù)的值，因為搜尋空間的維度由4降為2，最正確化會更有效率上述範(fàn)例的改良 xxeaeay21211a2an運用上述混成Hybrid的方法，函式 errorMeasure1

25、須改寫成 errorMeasure2n範(fàn)例10-10： errorMeasure2.mnlambda 是非線性參數(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(lambda(2)*x);a = Ay;y2 = a(1)*exp(lambda(1)*x)+a(2)*exp(lambda(2)

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

27、);legend(Sample data, Regression curve);fprintf(誤差平方和 = %dn, sum(y-y2).2); n此種混成法可以產(chǎn)生較低的誤差平方和，同時所需的計算時間也比較短?；斐煞ü?fàn)例3/3 計算時間 = 0.02誤差平方和 = 1.477226e-000.81.820123456 Sample dataRegression curven我們亦可利用變形法Transformation，將一數(shù)學(xué)模型轉(zhuǎn)換成只包含線性參數(shù)的模型。n假設(shè)一模型為： n取自然對數(shù)，可得：n ln(a) 及 b 成為線性參數(shù)，我們即可用最

28、小平方法找出其最正確值n範(fàn)例10-12： transformFit01.m變形法 bxaey bxay lnlnload data2.txtx = data2(:, 1);% 知資料點的 x 座標(biāo)y = data2(:, 2);% 知資料點的 y 座標(biāo)A = ones(size(x) x;a = 4.3282b =-1.8235誤差平方和 = 8.744185e-001變形法範(fàn)例1/2theta = Alog(y);subplot(2,1,1)plot(x, log(y), o, x, A*theta); xlabel(x); ylabel(ln(y);title(ln(y) vs. x);l

29、egend(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);n第一個小圖是ln(y)對x的作圖，n第二個小則是 y對 x的作圖。n經(jīng)由變形法之後，此最小平方法所得到的最小總

30、平方誤差是n而不是原模型的總平方誤差：n通常E為最小值時，E不一定是最小值，但亦離最小值不遠矣！變形法範(fàn)例2/0.81.82-3-2-1012xln(y)ln(y) vs. x Actual valuePredicted value00.81.82012345xyy vs. x Actual valuePredicted value21lnlnmiiibxayEmibxiiaeyE12n假設(shè)要求取原始E的最小值，可再用 fminsearch ，並以變形法得到的 a 及 b 為搜尋的起點n範(fàn)例10-13： trans

31、formFit02.m變形法之改進範(fàn)例1/2 load data2.txtx = data2(:, 1);% 知資料點的 x 座標(biāo)y = data2(:, 2);% 知資料點的 y 座標(biāo)A = ones(size(x) x;theta = Alog(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(

32、1)*exp(theta(2)*x);誤差平方和 = 1.680455e-001變形法之改進範(fàn)例2/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);00.81.8200.511.522.533.544.55xyy vs. x Actual valuePredicted valuen由上述範(fàn)例可以看出，我們可以先運用變形法，先找出大約的參數(shù)值，再用 fminsearch 來對誤差平方和進行最小化，因此得到的誤差平方和，比只用變形法還要小。n以下是一些變