曲线拟合与回归分析.ppt

MATLAB 程式設計進階篇曲線擬合與迴歸分析張智星(Roger Jang)jangmirlab.orghttp:/mirlab.org/jang清大資工系 多媒體檢索實驗室資料擬和簡介n資料擬合（Data Fitting）n給定一組資料（含輸入及輸出），建立一個數學模型，來逼近此資料的輸入輸出特性n如果此資料包含一維輸入及輸出，則此數學模型可以表示成一條曲線，在此情況下又稱為曲線擬合（Curve Fitting）n迴歸分析（Regression Analysis）n使用統計的方法來進行資料擬和，並分析每一個變數的統計特性，此過程稱為迴歸分析曲線擬合簡介 n曲線擬合（Curve Fitting）n欲建立的數學模型是單輸入、單輸出（Single-input Single-output，簡稱SISO），其特性可用一條曲線來表示 n迴歸分析之分類n若模型是線性模型，則此類問題稱為線性迴歸（Linear Regression）n若模型是非線性模型，則稱為非線性迴歸（Nonlinear Regression）。n觀察資料是美國自 1790 至 1990 年（以 10 年為一單位）的總人口，此資料可由載入檔案 census.mat 得到：n範例10-1：censusPlot01.m 曲線擬合：美國人口範例 load census.mat%載入人口資料plot(cdate,pop,o);%cdate 代表年度，pop 代表人口總數xlabel(年度);ylabel(美國人口總數);n模型選取n由上圖資料分佈，可猜測這適合的曲線可能是二次拋物線n其中y為輸出，x為輸入，則為此模型的參數。由於參數相對於y呈線性關係，所以此模型為稱線性模型n目標n找出最好的參數值，使得模型輸出與實際資料越接近越好，此過程即稱為線性迴歸曲線擬合之模型選取 n曲線擬和的平方誤差n假設觀察資料可寫成 ,i=121。當輸入為 時，實際輸出為 。n模型的預測值為 n平方誤差：n總平方誤差 是參數 的函數，這也是我們要最小化的目標函數，可表示如下：曲線擬合之目標函數 n求得參數 、的最佳值n求出 對 、的導式，令其為零，即可解出 、的最佳值。n總平方誤差 為 的二次式n導式 、及 為 的一次式n令上述導式為零之後，我們可以得到一組三元一次線性聯立方程式，就可以解出參數 、的最佳值。目標函數之求解 n假設 21 個觀察點均通過此拋物線，將這 21 個點帶入拋物線方程式，得到下列21個等式：n亦可寫成n其中 、為已知，為未知向量。矩陣表示法 n觀察n上述21個方程式，只有 3 個未知數 ，所以通常不存在一組解來滿足這 21 個方程式。n在一般情況下，只能找到一組 ，使得等號兩邊的差異為最小，此差異可寫成 此即為前述的總平方誤差nMATLAB 提供一個簡單方便的左除（）指令，來解出最佳的 ，使得總平方誤差為最小。MATLAB的最小平方解n利用左除來算出最佳的參數值，並同時畫出具有最小平方誤差的二次曲線n範例10-2：census01.m曲線擬合運算範例 load census.mat%載入人口資料plot(cdate,pop,o);%cdate 代表年度，pop 代表人口總數A=ones(size(cdate),cdate,cdate.2;y=pop;theta=Ay;%利用左除，找出最佳的 theta 值plot(cdate,pop,o,cdate,A*theta,-);legend(實際人口數,預測人口數);xlabel(年度);ylabel(美國人口總數);曲線擬合結果 n由上述範例，我們可以找出最佳的 n因此具有最小平方誤差的拋物線可以寫成：提示：左除及右除n左除的概念，可記憶如下：原先的方程式是 A*theta=y，我們可將 A移項至等號右邊，而得到 theta=Ay。必須小心的是：原先 A 在乘式的第一項，所以移到等號右邊後，A 仍然必須是除式的第一項。n若我們要解的方程式是 theta*A=y，則同樣的概念可得到最小平方解 theta=y/A。n根據上拋物線數學模型，我們可以預測美國在 2000 年的人口總數為：n範例10-3：census02.m以模型預測人口總數load census.mat%載入人口資料A=ones(size(cdate),cdate,cdate.2;theta=Apop;%利用左除，找出最佳的 theta 值t=2000;pop2000=1,t,t2*theta;%在 2000 年美國人口線數預測值t=2010;pop2010=1,t,t2*theta;%在 2010 年美國人口線數預測值fprintf(美國人口在2000年的預測值=%g（百萬人）n,pop2000);fprintf(美國人口在2010年的預測值=%g（百萬人）n,pop2010);美國人口在2000年的預測值=274.622（百萬人）美國人口在2010年的預測值=301.824（百萬人）n上述例子推廣，得到一個 n 次多項式：n利用多項式的數學模型來進行曲線擬合，通稱為多項式擬合（Polynomial Fitting）n由於多項式擬和的應用面很廣，MATLAB 提供npolyfit 指令來找出多項式最佳參數nPolyval 指令來進行多項式的求值多項式擬和 n使用 polyfit&polyval 使程式碼更加簡潔n範例10-4：census03.mnpolyfit(cdate,pop,2)中的 2 代表用到的模型是 2 次多項式load census.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的範例 nMATLAB 下輸入census，可對 census 資料進行曲線擬合的結果，如下：n上述圖形可以看出，當多項式的次數越來越高時，外插常會出現不可信的結果。n這表示選用的模型參數太高，雖然誤差的平方和變小了，但是預測的可靠度也下降了。模型複雜度對預測的影響 n線性迴歸的成功與否，與所選取的模型息息相關n模型所含的參數越多，平方誤差會越小n若參數個數等於資料點個數，平方誤差會等於零，但這並不表示預測會最準，因為資料點含有雜訊n完全吻合資料的模型亦代表此模型受雜訊的影響最大，預測之準確度也會較差n如何選取模型，是線性迴歸的一個重大課題！nLeave-one-out method線性迴歸的模型選取n多輸入、單輸出的線性迴歸數學模型寫成n其中 為輸入，y 為輸出，、為此模型的參數，則是已知的函數，稱為基底函數（Basis Functions）n所給的資料點為 ，稱為取樣資料（Sample Data）或訓練資料（Training Data）。多輸入之線性迴歸n將上述資料點帶入模型後可得：n或可表示成矩陣格式：矩陣表示法 n由於 (即資料點個數遠大於可變參數個數)，欲使上式成立，須加上一誤差向量 e：n平方誤差則可寫成n求 的最佳值n直接取 對 的偏微分，並令其等於零，即可得到一組 n 元一次的線性聯立方程式n用矩陣運算來表示，的最佳值可表示成 n也可以使用 MATLAB 的左除來算出 的最佳值，即 。平方誤差和的最小化 n理論上，最佳的 值為 ，但是 容易造成電腦內部計算的誤差，MATLAB 實際在計算左除時，會依照矩陣 A 的特性而選用最佳的方法，因此可以得到較穩定且正確的數值解。n欲知如何推導最佳的 值，可參考筆者另一著作：NeuroFuzzy and Soft Computing，Prentice Hall，1997。提示n在 MATLAB 下輸入 peaks，可以畫出一個凹凸有致的曲面，如下：n此函數的方程式如下：曲面擬合範例（1/6）n在下列說明中，假設：n數學模型的基底函數已知n訓練資料包含正規分佈的雜訊n上述函數可寫成：n其中我們假設 、和 是未知參數，n 則是平均為零、變異為 1 的正規分佈雜訊。曲面擬合範例（2/6）n若要取得 100 筆訓練資料n範例10-5：peaks01.mnrandn 指令的使用即在加入正規分佈雜訊。上圖為我們收集到的訓練資料，由於雜訊很大，所以和原先未帶雜訊的圖形差異很大。曲面擬合範例（3/6）pointNum=10;xx,yy,zz=peaks(pointNum);zz=zz+randn(size(zz);%加入雜訊 surf(xx,yy,zz);axis tightn現在我們要用已知的基底函數，來找出最佳的 、和 n範例10-6：peaks02.mn由此找出的 值和最佳值 相當接近。曲面擬合範例（4/6）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=Az%最佳的 theta 值theta=3.0088 -10.0148 -0.2924 n根據上求得之參數，可以輸入較密的點，得到迴歸後的曲面n範例10-7：peaks03.m曲面擬合範例（5/6）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=Az;%最佳的 theta 值%畫出預測的曲面pointNum=31;xx,yy=meshgrid(linspace(-3,3,pointNum),linspace(-3,3,pointNum);n在上圖中，我們猜對了基底函數，因此得到非常好的曲面擬合。n只要基底函數正確，而且雜訊是正規分佈，那麼當資料點越來越多，上述的最小平方法就可以逼近參數的真正數值。曲面擬合範例（6/6）x=xx(:);y=yy(:);%轉為行向量 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 n非線性迴歸（Nonlinear Regression）是一個比較困難的問題，原因如下：n無法一次找到最佳解。n無法保證能夠找到最佳解。n須引用各種非線性最佳化的方法。n各種相關數學性質並不明顯。n以數學來描述，假設所用的數學模型是n其中 是輸入向量，是可變非線性函數，y 是輸出變數。n總平方誤差為非線性迴歸 n用一般最佳化（Optimization）的方法，來找出 的最小值，例如n梯度下降法（Gradient Descent）nSimplex下坡式搜尋（Simplex Downhill search）：此方法即為fminsearch指令所採用的方法n假設所用的數學模型為n其中、為線性參數，但1、2 為非線性參數n總平方誤差可表示：n欲找出使 為最小的 、1及2，需將 E 寫成一函式，並由其它最佳化的方法來求出此函式的最小值。非線性迴歸：誤差值的最小化n我們使用errorMeasure1.m來計算誤差值n範例10-8：errorMeasure01.mn其中 theta 是參數向量，包含了 、1 及 2，data 則是觀察到的資料點，傳回的值 則是總平方誤差。使用fminsearch的範例（1/3）function 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);n欲求出 的最小值，我們可使用 fminsearch 指令n範例10-9：nonlinearFit01.m使用fminsearch的範例（2/3）load 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)*exp(theta(4)*x);plot(x,y,ro,x,y2,b-);legend(Sample data,Regression curve);fprintf(誤差平方和=%dn,sum(y-y2).2);n上圖的曲線為 fminsearch 指令產生的迴歸曲線。nfminsearch 指令是一個使用 Simplex 下坡式搜尋法（Downhill Simplex Search）的最佳化方法，用來找出 errorMeasure1 的極小值，並傳回 theta 的最佳值。使用fminsearch的範例（3/3）計算時間=0.03誤差平方和=5.337871e-001n改良方向n上述方法把所有參數全部視為非線性參數。混成法將線性與非線性參數分開，各用不同的方法來處理。n上例數學模型為n 、線性參數：最小平方法，即左除或n1、2 非線性參數：Simplex 下坡式搜尋（即 fminsearch）n混成法的優點n最小平方法能夠在非線性參數固定的情況下，一次找到最好的線性參數的值，因為搜尋空間的維度由4降為2，最佳化會更有效率上述範例的改良 n使用上述混成（Hybrid）的方法，函式 errorMeasure1 須改寫成 errorMeasure2n範例10-10：errorMeasure2.mnlambda 是非線性參數向量，data 仍是觀察到的資料點，a 是利用最小平方法算出的最佳線性參數向量，傳回的 squareError 仍是總平方誤差混成法範例（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)*x);squaredError=sum(y-y2).2);n欲用此混成法求出誤差平方和的最小值n範例10-11：nonlinearFit02.m混成法範例（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-);legend(Sample data,Regression curve);fprintf(誤差平方和=%dn,sum(y-y2).2);n此種混成法可以產生較低的誤差平方和，同時所需的計算時間也比較短。混成法範例（3/3）計算時間=0.02誤差平方和=1.477226e-001n我們亦可利用變形法（Transformation），將一數學模型轉換成只包含線性參數的模型。n假設一模型為：n取自然對數，可得：n ln(a)及 b 成為線性參數，我們即可用最小平方法找出其最佳值n範例10-12：transformFit01.m變形法 load data2.txtx=data2(:,1);%已知資料點的 x 座標y=data2(:,2);%已知資料點的 y 座標A=ones(size(x)x;a=4.3282b=-1.8235誤差平方和=8.744185e-001變形法範例（1/2）theta=Alog(y);subplot(2,1,1)plot(x,log(y),o,x,A*theta);xlabel(x);ylabel(ln(y);title(ln(y)vs.x);legend(Actual value,Predicted value);a=exp(theta(1)%辨識得到之參數b=theta(2)%辨識得到之參數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經由變形法之後，此最小平方法所得到的最小總平方誤差是n而不是原模型的總平方誤差：n通常E為最小值時，E不一定是最小值，但亦離最小值不遠矣！變形法範例（2/2）n若要求取原始E的最小值，可再用 fminsearch，並以變形法得到的 a 及 b 為搜尋的起點n範例10-13：transformFit02.m變形法之改進範例（1/2）load data2.txtx=data2(:,1);%已知資料點的 x 座標y=data2(:,2);%已知資料點的 y 座標A=ones(size(x)x;theta=Alog(y);a=exp(theta(1)%辨識得到之參數b=theta(2)%辨識得到之參數theta0=a,b;%fminsearch 的啟始參數theta=fminsearch(errorMeasure3,theta0,data2);x=data2(:,1);y=data2(:,2);y2=theta(1)*exp(theta(2)*x);誤差平方和=1.680455e-001變形法之改進範例（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);n由上述範例可以看出，我們可以先使用變形法，先找出大略的參數值，再用 fminsearch 來對誤差平方和進行最小化，因此得到的誤差平方和，比只用變形法還要小。n以下是一些變形法可用的非線性模型，以及相關的轉換方法：可使用變形法的數學模型（1/3）編號非線性模型轉換後的線性模型參數轉換公式1234567可使用變形法的數學模型（2/3）8910可使用變形法的數學模型（3/3）n曲線擬合包含下列幾個步驟：n觀察資料，並剔除明顯不合理的資料（這些資料在統計學上稱為 Outliers）。n根據資料，選定數學模型及相關參數。n根據線性或非線性迴歸的各種方法，以及一組給定的訓練資料（Training Data），算出參數的最佳值。n觀察模型的誤差，以及模型對於其它測試資料（Test Data）的效能，以決定此模型的適用性。若適用，則停止此演算法。反之，若不適用，則根據模型的對於訓練及測試資料的誤差程度，重新修正模型，並回到步驟 3。曲線擬合工具箱的使用（1/4）n以上步驟，需要經驗，且必須反覆進行，可能耗費大量時間，有鑑於此，MathWorks 公司在 MATLAB 6.x 後，推出了曲線擬合工具箱（Curve Fitting Toolbox），讓使用者能以 GUI（圖形使用者介面）的方式，來進行曲線擬合，並能快速地檢視擬合的結果和成效。曲線擬合工具箱的使用（2/4）n說明曲線擬合工具箱的使用n首先，先載入 enso.mat，裡面包含兩個變數：nmonth：每個資料點發生的相對月份npressure：在復活島（Easter Island）和澳洲的達爾文（Darwin）兩地的大氣壓力差值，取其整個月的平均值。此差值會導引整個南半球的貿易風（Trade Winds）流向n根據這兩個變數，就可以呼叫曲線擬合工具箱來進行曲線的分析與擬合n範例10-14：cftool01.m曲線擬合工具箱的使用（3/4）load enso.mat%載入 month 和 pressure 變數cftool(month,pressure);%呼叫曲線擬合工具箱n此時此工具箱會將資料點畫出來，並將相關的操作介面顯示如上圖。曲線擬合工具箱的使用（4/4）