MATLAB程式设计进阶篇曲线拟合与回归分析(精)

MATLAB程式設計進階篇曲線擬合與迴歸分析張智星(RogerJang)[e-mail]清大資工系多媒體檢索實驗室資料擬和簡介資料擬合（DataFitting）給定一組資料（含輸入及輸出），建立一個數學模型，來逼近此資料的輸入輸出特性如果此資料包含一維輸入及輸出，則此數學模型可以表示成一條曲線，在此情況下又稱為曲線擬合（CurveFitting）迴歸分析（RegressionAnalysis）使用統計的方法來進行資料擬和，並分析每一個變數的統計特性，此過程稱為迴歸分析曲線擬合簡介曲線擬合（CurveFitting）欲建立的數學模型是「單輸入、單輸出」（Single-inputSingle-output，簡稱SISO），其特性可用一條曲線來表示迴歸分析之分類若模型是線性模型，則此類問題稱為線性迴歸（LinearRegression）若模型是非線性模型，則稱為非線性迴歸（NonlinearRegression）。觀察資料是美國自1790至1990年（以10年為一單位）的總人口，此資料可由載入檔案census.mat得到：範例10-1：censusPlot01.m曲線擬合：美國人口範例loadcensus.mat%載入人口資料plot(cdate,pop,'o');%cdate代表年度，pop代表人口總數xlabel('年度');ylabel('美國人口總數');175018001850190019502000050100150200250年度美國人口總數模型選取由上圖資料分佈，可猜測這適合的曲線可能是二次拋物線其中y為輸出，x為輸入，則為此模型的參數。由於參數相對於y呈線性關係，所以此模型為稱線性模型目標找出最好的參數值，使得模型輸出與實際資料越接近越好，此過程即稱為線性迴歸曲線擬合之模型選取2210210),,;(xaxaaaaaxfy21,0,aaa曲線擬和的平方誤差假設觀察資料可寫成,i=1~21。當輸入為時，實際輸出為。模型的預測值為平方誤差：總平方誤差是參數的函數，這也是我們要最小化的目標函數，可表示如下：曲線擬合之目標函數),(iiyxixiyiy2210210),,;(iiixaxaaaaaxf2)(iixfy211222102112210)(),,(iiiiiiixaxaayxfyaaaE21,0,aaaE求得參數、、的最佳值求出對、、的導式，令其為零，即可解出、、的最佳值。總平方誤差為的二次式導式、及為的一次式令上述導式為零之後，我們可以得到一組三元一次線性聯立方程式，就可以解出參數、、的最佳值。目標函數之求解0a1a2a0aE1aE2aEE0a1a2aE0a1a2a0a1a2a21,0,aaa21,0,aaa假設21個觀察點均通過此拋物線，將這21個點帶入拋物線方程式，得到下列21個等式：亦可寫成其中、為已知，為未知向量。矩陣表示法212212211022222101212110yxaxaayxaxaayxaxaayAyyyxxxxxx212132122121222211111Ay觀察上述21個方程式，只有3個未知數，所以通常不存在一組解來滿足這21個方程式。在一般情況下，只能找到一組，使得等號兩邊的差異為最小，此差異可寫成此即為前述的總平方誤差MATLAB提供一個簡單方便的「左除」（\）指令，來解出最佳的，使得總平方誤差為最小。MATLAB的最小平方解T321,,)()(2AyAyAyTE利用「左除」來算出最佳的參數值，並同時畫出具有最小平方誤差的二次曲線範例10-2：census01.m曲線擬合運算範例loadcensus.mat%載入人口資料plot(cdate,pop,'o');%cdate代表年度，pop代表人口總數A=[ones(size(cdate)),cdate,cdate.^2];y=pop;theta=A\y;%利用「左除」，找出最佳的theta值plot(cdate,pop,'o',cdate,A*theta,'-');legend('實際人口數','預測人口數');xlabel('年度');ylabel('美國人口總數');曲線擬合結果由上述範例，我們可以找出最佳的因此具有最小平方誤差的拋物線可以寫成：00654.0,51.23,21130,,210aaa2221000654.051.2321130)(xxxaxaaxfy175018001850190019502000050100150200250年度美國人口總數實際人口數預測人口數提示：左除及右除左除的概念，可記憶如下：原先的方程式是A*theta=y，我們可將A移項至等號右邊，而得到theta=A\y。必須小心的是：原先A在乘式的第一項，所以移到等號右邊後，A仍然必須是除式的第一項。若我們要解的方程式是theta*A=y，則同樣的概念可得到最小平方解theta=y/A。根據上拋物線數學模型，我們可以預測美國在2000年的人口總數為：範例10-3：census02.m以模型預測人口總數loadcensus.mat%載入人口資料A=[ones(size(cdate)),cdate,cdate.^2];theta=A\pop;%利用「左除」，找出最佳的theta值t=2000;pop2000=[1,t,t^2]*theta;%在2000年美國人口線數預測值t=2010;pop2010=[1,t,t^2]*theta;%在2010年美國人口線數預測值fprintf('美國人口在2000年的預測值=%g（百萬人）\n',pop2000);fprintf('美國人口在2010年的預測值=%g（百萬人）\n',pop2010);美國人口在2000年的預測值=274.622（百萬人）美國人口在2010年的預測值=301.824（百萬人）上述例子推廣，得到一個n次多項式：利用多項式的數學模型來進行曲線擬合，通稱為「多項式擬合（PolynomialFitting）」由於多項式擬和的應用面很廣，MATLAB提供polyfit指令來找出多項式最佳參數Polyval指令來進行多項式的求值多項式擬和nnxaxaaxfy10)(使用polyfit&polyval使程式碼更加簡潔範例10-4：census03.mpolyfit(cdate,pop,2)」中的2代表用到的模型是2次多項式loadcensus.mat%載入人口資料theta=polyfit(cdate,pop,2);%進行二次多項式擬合，找出theta值fprintf('2000年的預測值=%g（百萬人）\n',polyval(theta,2000));fprintf('2010年的預測值=%g（百萬人）\n',polyval(theta,2010));在2000年的預測值=274.622（百萬人）在2010年的預測值=301.824（百萬人）使用polyfit&polyval的範例MATLAB下輸入「census」，可對census資料進行曲線擬合的結果，如下：上述圖形可以看出，當多項式的次數越來越高時，「外插」常會出現不可信的結果。這表示選用的模型參數太高，雖然誤差的平方和變小了，但是預測的可靠度也下降了。模型複雜度對預測的影響線性迴歸的成功與否，與所選取的模型息息相關模型所含的參數越多，平方誤差會越小若參數個數等於資料點個數，平方誤差會等於零，但這並不表示預測會最準，因為資料點含有雜訊完全吻合資料的模型亦代表此模型受雜訊的影響最大，預測之準確度也會較差如何選取模型，是線性迴歸的一個重大課題！Leave-one-outmethod線性迴歸的模型選取「多輸入、單輸出」的線性迴歸數學模型寫成其中為輸入，y為輸出，、、…、為此模型的參數，，則是已知的函數，稱為基底函數（BasisFunctions）所給的資料點為，稱為取樣資料（SampleData）或訓練資料（TrainingData）。多輸入之線性迴歸1a2ananifi1),(xmiyii1),,(x)()()()(2211xxxxnnfafafafyxy將上述資料點帶入模型後可得：或可表示成矩陣格式：矩陣表示法)()()()()()()()(221122111mmmm1111xxxxxxxxnnmnnfafafafyfafafafyymnAmnmnyyaaffff111111)()()()(xxxx由於(即資料點個數遠大於可變參數個數)，欲使上式成立，須加上一誤差向量e：平方誤差則可寫成求的最佳值直接取對的偏微分，並令其等於零，即可得到一組n元一次的線性聯立方程式用矩陣運算來表示，的最佳值可表示成也可以使用MATLAB的「左除」來算出的最佳值，即。平方誤差和的最小化nmyeθA)()()(2AyAyeeeETT)(EyAAATT1yA\理論上，最佳的值為，但是容易造成電腦內部計算的誤差，MATLAB實際在計算「左除」時，會依照矩陣A的特性而選用最佳的方法，因此可以得到較穩定且正確的數值解。欲知如何推導最佳的值，可參考筆者另一著作：Neuro–FuzzyandSoftComputing，PrenticeHall，1997。提示yAAATT11AAT在MATLAB下輸入peaks，可以畫出一個凹凸有致的曲面，如下：此函數的方程式如下：曲面擬合範例（1/6）222222)1(53)1(231510)1(3yxyxyxeeyxxexz在下列說明中，假設：數學模型的基底函數已知訓練資料包含正規分佈的雜訊上述函數可寫成：其中我們假設、和是未知參數，n則是平均為零、變異為1的正規分佈雜訊。曲面擬合範例（2/6）nyxfyxfyxfneyxfyxfneeyxxexzyxyxyxyx),(),(),(31),(10),(331510)1(3332211)1(21)1(53)1(222222222123若要取得100筆訓練資料範例10-5：peaks01.mrandn指令的使用即在加入正規分佈雜訊。上圖為我們收集到的訓練資料，由於雜訊很大，所以和原先未帶雜訊的圖形差異很大。曲面擬合範例（3/6）pointNum=10;[xx,yy,zz]=peaks(pointNum);zz=zz+randn(size(zz));%加入雜訊surf(xx,yy,zz);axistight-3-2-10123-202-4-20246現在我們要用已知的基底函數，來找出最佳的、和範例10-6：peaks02.m由此找出的值和最佳值相當接近。曲面擬合範例（4/6）123pointNum=10;[xx,yy,zz]=peaks(pointNum);zz=zz+randn(size(zz))/10;%加入雜訊x=xx(:);%轉為行向量y=yy(:);%轉為行向量z=zz(:);%轉為行向量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=A\z%最佳的theta值theta=3.0088-10.0148-0.