13.數(shù)據(jù)插值與擬合實(shí)際中，通常需要處理實(shí)驗(yàn)或測(cè)量得到的離散數(shù)據(jù)（點(diǎn)）。插值與擬合方法就是要通過離散數(shù)據(jù)去確定一個(gè)近似函數(shù)（曲線或曲面），使其與已知數(shù)據(jù)有較高的擬合精度。.如果要求近似函數(shù)經(jīng)過所已知的所有數(shù)據(jù)點(diǎn)，此時(shí)稱為插值問題（不需要函數(shù)表達(dá)式）。.如果不要求近似函數(shù)經(jīng)過所有數(shù)據(jù)點(diǎn)，而是要求它能較好地反映數(shù)據(jù)變化規(guī)律，稱為數(shù)據(jù)擬合（必須有函數(shù)表達(dá)式）。插值與擬合都是根據(jù)實(shí)際中一組已知數(shù)據(jù)來構(gòu)造一個(gè)能夠反映數(shù)據(jù)變化規(guī)律的近似函數(shù)。區(qū)別是：【插值】不一定得到近似函數(shù)的表達(dá)形式，僅通過插值方法找到未知點(diǎn)對(duì)應(yīng)的值?！緮M合】要求得到一個(gè)具體的近似函數(shù)的表達(dá)式。因此，當(dāng)數(shù)據(jù)量不夠，但已知已有數(shù)據(jù)可信，需要補(bǔ)充數(shù)據(jù)，此時(shí)用【插值】。當(dāng)數(shù)據(jù)基本夠用，需要尋找因果變量之間的數(shù)量關(guān)系（推斷出表達(dá)式），進(jìn)而對(duì)未知的情形作預(yù)測(cè)，此時(shí)用【擬合工、數(shù)據(jù)插值根據(jù)選用不同類型的插值函數(shù)，逼近的效果就不同，一般有：（1）拉格朗日插值（lagrange插值）（2）分段線性插值（3）Hermite（4）三次樣條插值Matlab插值函數(shù)實(shí)現(xiàn)：（1）interp1（） 一維插值（2）intep2（）二維插值（3）interp3（） 三維插值（4）intern（）n維插值1.一維插值（自變量是1維數(shù)據(jù)）語法：yi=interp1（x0,y0,xi,‘method,）其中，x0,y0為原離散數(shù)據(jù)（x0為自變量，y0為因變量）；xi為需要插值的節(jié)點(diǎn)，method為插值方法。注：（1）要求x0是單調(diào)的，xi不超過x0的范圍；（2）插值方法有‘nearest’ 最鄰近插值；‘linear’ 線性插值；‘spline’ 三次樣條插值；‘cubic' 三次插值；默認(rèn)為分段線性插值。例1從1點(diǎn)12點(diǎn)的11小時(shí)內(nèi)，每隔1小時(shí)測(cè)量一次溫度，測(cè)得的溫度的數(shù)值依次為：5，8，9，15，25，29，31，30，22，25，27,24.試估計(jì)每隔1/10小時(shí)的溫度值。代碼：hours=1:12;temps=[589152529313022252724];h=1:0.1:12;t=interp1(hours,temps,h,lsplinel);plot(hours,temps,'+',h,t,hours,temps,'r:')%作圖，原散點(diǎn)數(shù)據(jù)用+標(biāo)記xlabel('Houri),ylabel(iDegreesCelsius')運(yùn)行結(jié)果：S^^PCS^^PCseemsu2.二維插值（自變量是2維數(shù)據(jù)）語法： zi=interp2（x0,y0,z0,xi,yi,‘method'）其中，x0,y0,z0為原離散數(shù)據(jù)（x0,y0為自變量，z0為因變量）；xi,yi為需要插值的節(jié)點(diǎn)，method為插值方法。注：（1）要求x0,y0是單調(diào)的，xi,yi不超過x0,y0的范圍，注意：要求xi取為行向量，yi取列向量；（2）插值方法有‘nearest’ 最鄰近插值；'linear’ 雙線性插值；‘cubic’一一雙三次插值；默認(rèn)為雙線性插值。例2山區(qū)地貌：在某山區(qū)測(cè)得一些地點(diǎn)的高度如下表：（平面區(qū)域?yàn)椋?200<=x<=4000,1200<=y<=3600）X12001600200024002800320036004000120011301250128012301040900500700160013201450142014001300700900850200013901500150014009001100106095024001500120011001350145012001150101028001500120011001550160015501380107032001500155016001550160016001600155036001480150015501510143013001200980試作出該山區(qū)的地貌圖和等高線圖。代碼：x0=1200:400:4000;y0=1200:400:3600;height=[11301250128012301040900500700;13201450142014001300700900850;139015001500140090011001060950;15001200110013501450120011501010;15001200110015501600155013801070;15001550160015501600160016001550;1480150015501510143013001200980];

mesh(x0,y0,height)title(，原數(shù)據(jù)散點(diǎn)圖')；pause %等待，按任意鍵切換到新圖形xi=1200:5:4000;yi=1200:5:3600;zi=interp2(x0,y0,height,xil,yi,lcubicl);mesh(xi,yi,zi)title('cubic插值圖');運(yùn)行結(jié)果：另外，二維插值函數(shù)還有g(shù)riddata函數(shù)，格式類似cz=griddata（x0,y0,z0,cx,cy,‘method'）插值方法有'nearest’ 最鄰近插值；‘linear' 雙線性插值；‘cubic’一一雙三次插值；‘v4’一一Matlab提供的插值方法；默認(rèn)為雙線性插值。與interp2區(qū)別是，interp2的插值數(shù)據(jù)必須是矩形域，即已知數(shù)據(jù)點(diǎn)（x,y）組成規(guī)則的矩陣（或稱之為柵格），可使用meshgid生成。而griddata函數(shù)的已知數(shù)據(jù)點(diǎn)（X，Y）不要求規(guī)則排列，特別是對(duì)試驗(yàn)中隨機(jī)沒有規(guī)律采取的數(shù)據(jù)進(jìn)行插值具有很好的效果。例3在某海域測(cè)得一些點(diǎn)（x,y）處的水深z由下表給出，船的吃水深度為5英尺，在矩形區(qū)域（75,200）X（-50,150）里的哪些地方船要避免進(jìn)入.xyz1297.54140141.58103.5236881478185.522.56195137.5810585.58157.5-6.59107.5-81977388156.58162-66.59162844117.5-33.59

作海底曲面圖，作出水深小于5的海域范圍，即z=5的等高線。代碼：x=[129.0140.0103.588.0185.5195.0105.5157.5107.577.081.0162.0162.0117.5];y=[7.5141.523.0147.022.5137.585.5-6.5-813.056.5-66.584.0-33.5];z=[48686889988949];x1=75:1:200;y1=-50:1:150;[x1,y1]=meshgrid(x1,y1);z1=griddata(x,y,z,x1,y1,lv4l);meshc(x1,y1,z1)pausez1(z1>=5)=nan;%將水深大于5的置為nan,這樣繪圖就不會(huì)顯示出來meshc(x1,y1,z1)運(yùn)行結(jié)果：1210-8-6-2：1501002005015020050100-5050

5一4.5-4-3.5：1505一4.5-4-3.5：150200-50 50二、數(shù)據(jù)擬合擬合最常用的方法為最小二乘法。原理是：因變量的實(shí)際值與擬合值之差，稱為殘差（擬合誤差），將所有殘差平方后相加，即得誤差平方和，“最好”擬合效果就是使誤差平方和最小，于是運(yùn)用極值原理，將“求最好擬合問題”轉(zhuǎn)換為“求誤差平方和最小”問題。Matlab實(shí)現(xiàn)：1.多項(xiàng)式擬合fX)=%Xm+???+amx+am+1語法:a=polyfit(x0,y0,m);y=polyval(a,x);注：（1）x,y的長度相同，m為擬合多項(xiàng)式次數(shù)，m=1時(shí)，一次多項(xiàng)式對(duì)應(yīng)一條直線，也稱為線性擬合（一元線性回歸）；（2）a返回的是擬合多項(xiàng)式的系數(shù)：[al,a2,…,am];（3）polyval函數(shù)用來計(jì)算擬合多項(xiàng)式在x點(diǎn)處的取值；（4）擬合多項(xiàng)式的次數(shù)一般為1-3次，不宜過高，否則有龍格現(xiàn)象（兩端出現(xiàn)很大誤差）。例41949年—1994年我國人口數(shù)據(jù)資料如下:年份x01949195419591964196919741979198419891994人口數(shù)y05.46.06.77.08.19.19.810.311.311.8分析我國人口增長的規(guī)律，預(yù)報(bào)2005年我國人口數(shù)。代碼：x0=[1949195419591964196919741979198419891994];y0=[5.46.06.77.08.19.19.810.311.311.8];a=polyfit(x0,y0,1);%線性擬合x1=[1949:0.5:1994];y1=polyval(a,x1);%或者直接用多項(xiàng)式計(jì)算：y1=a(1)*x1+a(2);y1_2005=polyval(a,2005)b=polyfit(x0,log(y0),1);%指數(shù)函數(shù)擬合y2=exp(b(2))*exp(b(1)*x1);y2_2005=exp(b(2))*exp(b(1)*2005)plot(x0,y0,,*,)holdonplot(x1,y1,'--r')holdonplot(x1,y2,'-k')legend('原數(shù)據(jù)散點(diǎn)圖'，'線性擬合曲線'，'指數(shù)擬合曲線')運(yùn)行結(jié)果：

y1_2005=13.5080y2_2005= 15.0596注：(1)實(shí)際2005年人口是13.3億，可見線性擬合效果更好；(2)人口指數(shù)增長模型為：y=a?如，兩邊取ln得lny=lna+bx令Y=lny，A=lna，B=b，則轉(zhuǎn)化為直線方程：Y=A+Bx故可以作線性擬合得到A，B值，從而得到y(tǒng)=eA.eBx.一般可以轉(zhuǎn)化為直線方程的曲線有：曲線類型變換直線方程幕函數(shù)：p=2Y-Iny,X=JnxV=Ina+bX指數(shù)函數(shù):y=aebxY=Iny,X=xY=lua+bX雙曲函婁bl. b--a-\——y xy xY=a+bXS型曲線+ 1y~a-beJY-LX-e-xY=a+hX（Logistic生長曲線，就是S型曲線）2.多元線性回歸y=P+Bx+L+Bx+eBj為回歸系數(shù)，￡為殘差。n=1時(shí)為線性回歸。Matlab實(shí)現(xiàn)：[b,bint,r,rint,stats]=regress（y,x,a）;注：（1）整體是一個(gè)假設(shè)檢驗(yàn)（是否符合多元線性關(guān)系）：原假設(shè)H0：B『0（即不存在線性關(guān)系）備擇假設(shè)H1：B.#0（回歸式成立）；（2）若有常數(shù)項(xiàng)B0,則x的第一列必須全為1；a為顯著性水平的概率，默認(rèn)是0.05;b返回B0邛值的點(diǎn)估計(jì)，bint返回B0,BbR值的區(qū)間估計(jì)，r返回殘差的點(diǎn)估計(jì)，rint返回殘差的區(qū)間估計(jì)；stats返回4個(gè)數(shù)：相關(guān)系數(shù)平方r2（越大越好，若小于0.5結(jié)果無價(jià)值）F統(tǒng)計(jì)量（用于計(jì)算概率）拒絕原假設(shè)的概率（＜a時(shí)拒絕）（越小越好）剩余方差=殘差平方和/[數(shù)據(jù)個(gè)數(shù)-（n+1）]（越小越好）其中，n+1為回歸系數(shù)B0,B1，…Bn的個(gè)數(shù)。例5已知某動(dòng)物的進(jìn)食量與體重?cái)?shù)據(jù)為進(jìn)食量1.5234.57.59.110.512體重5.66.67.27.810.110.813.516.5擬合直線y=B0+B1x.代碼：food=[1.5234.57.59.110.512];n=length(food);x=[ones(n,1),food'];y=[5.66.67.27.810.110.813.516.5]1;[b,bint,r,rint,stat]=regress(y,x,0.05)運(yùn)行結(jié)果：b= 4.15750.8950 ——>回歸方程：y=4.1575+0.8950xbint=2.4692 5.8458——>4.1575￡[2.4692,5.8458]0.6644 0.6644——>0.8950￡[0.6644,0.6644]r= 0.1000——>殘差1:5.6-(4.1575+0.8950*1.5戶0.10.65250.35750.38500.77011.50210.05511.6024rint= -2.12832.3283——>0.1000￡[-2.1283,2.3283]-1.52822.8332-2.00742.7223-2.84452.0744-3.14011.5999-3.29350.2893-2.35192.24170.4846 2.7201stat= 0.9376 90.1871 0.0001 1.0220>r2=0.9376>0.5(說明結(jié)果很好)0.0001<a=0.05(拒絕原假設(shè)，接受備擇假設(shè))【殘差r的平方和/[數(shù)據(jù)個(gè)數(shù)-(n+1)]]：(0.10002+0.65252+???+1.60242)/[8-(1+1)]=1.02203.一般曲線擬合(非線性回歸)曲線形式y(tǒng)=f(a,x)已知，但其中a=[a1,a2,…]為未知的待定系數(shù)。例如，已知曲線形式為y=f(a,x)=a1x2+a2sinx+a/3，以及原始數(shù)據(jù)x1,…,xn，y1,…,yn，利用最小二乘法原理，求出讓誤差平方和[F(a,x1)-y1]2H +[F(a,xn)-yn]2取最小值的a1，a2,a3.從而得到擬合曲線方程。Matlab實(shí)現(xiàn)：Isqcurvefit和Isqnonlin函數(shù)。Isqcurvefit函數(shù)語法：先定義形式函數(shù)f(a,x)的M-文件fun.m;再調(diào)用[x,resnorm,residual]=Isqcurvefit('fun’,x0,xdata,ydata,lb,ub,options);注：(1)返回值x為擬合出的待定系數(shù)，resnorm返回殘差的平方和(越小越好)，residual返回殘差；(2)‘fun’為預(yù)定義的形式函數(shù)，xdata,ydata為已知數(shù)據(jù)，x0為迭代初始值(初始值選取的不同將影響擬合結(jié)果)，lb，ub是可選參數(shù)，為x指定范圍“下界和上界”；(3)options是可選參數(shù)，指定采用的優(yōu)化參數(shù)，語法：OPTIONS=OPTIMSET('PARAM1',VALUE1,'PARAM2',VALUE2,...)示例：options=optimset('TolFun',1e-3);注意：形式函數(shù)是以待定系數(shù) a和已知數(shù)據(jù)x作為參數(shù)，lsqcurvefit函數(shù)的參數(shù)中包含已知數(shù)據(jù)x，y.例6用下面一組數(shù)據(jù)擬合C?)=a+be-0.02kt中的待定系數(shù)a，b，k.t.j1002003004005006007008009001000cj(*10-3)4.544.995.355.655.906.106.266.396.506.59該問題即解最優(yōu)化問題：滿足下式的的a，b，k:

minXa+be-0-02ktj一c- j-jT代碼:(Isqcurvefit函數(shù))%%%%%%%%%%函數(shù)單獨(dú)存為m文件％%%%%%%%%%%%%%%%%%%%%%%%functionf=curvefun1(x,tdata)%tdata是向量f=x(1)+x(2)*exp(-0.02*x(3)*tdata);%其中x(1)=a;x(2)=b;x(3)=k;end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%tdata=[100:100:1000]';cdata=1e-03*[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59]';x0=[0.2,0?05,0.05]';[x,resnorm,residual]=lsqcurvefit('curvefun1',x0,tdata,cdata)C_fit=curvefun1(x,tdata)plot(tdata,cdata,'o',tdata,f,'r')運(yùn)行結(jié)果：x=0.0069-0.00290.0809 >擬合曲線：C(t)=0.0069—0.0029e-0.02x0.08091resnorm=5.3938e-007（殘差平方和）residual=1.0e-003*（殘差）-0.0971-0.1739-0.2164resnorm=5.3938e-007（殘差平方和）residual=1.0e-003*（殘差）-0.0971-0.1739-0.2164-0.2463-0.2666-0.2713-0.2651-0.2538-0.2436-0.2313-0.2651-0.2538-0.2436-0.2313cfit=0.00440.00480.00510.00540.00560.00580.00600.0061 0.0063cfit=0.00440.00480.00510.00540.00560.00580.00600.0061 0.00630.006476.5676.565.554.54100 200 300 400 500 600 700 800 900 1000x10-3Isqnonlin函數(shù)與Isqcurvefit函數(shù)用法類似，不同在于形式函數(shù)只以待定系數(shù)a作為參數(shù)，但函數(shù)體里用到已知數(shù)據(jù)x,y.而在調(diào)用擬合函數(shù)時(shí)不再需要已知數(shù)據(jù)x,y作為參數(shù)：[x,resnorm,residual]=Isqnonlin(‘fun’,x0,lb,ub,options);對(duì)于例6,調(diào)用lsqnonlin函數(shù)實(shí)現(xiàn)的代碼如下：(運(yùn)行結(jié)果一樣)代碼：(lsqnonlin函數(shù))%%%%%%%%%%%%%%%%函數(shù)單獨(dú)存為m文件％%%%%%%%%%%%%%%%%%functionf=curvefun2(x)tdata=[100:100:1000]';cdata=1e-03*[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59]';f=x(1)+x(2)*exp(-0.02*x(3)*tdata)-cdata;end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%x0=[0.2,0?05,0.05]';[x,resnorm,residua]=lsqnonlin('curvefun2',x0)nlinfit函數(shù)非線性擬合的通用函數(shù)，適用面比Isqcurvefit函數(shù)廣，例如可以做加權(quán)最小二乘擬合。語法：回歸：[beta,r,J]=nlinfit(x,y,‘model’,beta0);預(yù)測(cè)：[Y,DELTA]=nlpredci(‘model’，x,beta,r,J);注：(1)beta返回?cái)M合的待定系數(shù)；r返回殘差；J返回?cái)M合函數(shù)的雅可比矩陣；x、y為已知數(shù)據(jù)，分別為nXm矩陣和n維列向量，對(duì)一元非線性回歸，x為n維列向量；‘model,為預(yù)先定義的形式函數(shù)；beta0為指定的待定系數(shù)初始值；nlpredic函數(shù)用來求nlinfit所得的擬合函數(shù)在x處的預(yù)測(cè)值Y及預(yù)測(cè)值的顯著性為1-alpha的置信區(qū)間Y±DELTA.例7混凝土的抗壓強(qiáng)度隨養(yǎng)護(hù)時(shí)間的延長而增加，現(xiàn)將一批混凝土作成12個(gè)試塊，下表記錄了養(yǎng)護(hù)日期x(日)及抗壓強(qiáng)度y(kg/cm2)的數(shù)據(jù)：養(yǎng)護(hù)時(shí)間x234579121417212856抗壓強(qiáng)度y(+r)354247535965687376828699建立非線性回歸模型:y=a+kemx+ke—mx對(duì)其中的待定系數(shù)a,k1,k2,m進(jìn)行擬合(回歸)，并對(duì)結(jié)果做檢驗(yàn)。注：表中的+r代表加上一個(gè)[-0.5,0.5]之間的隨機(jī)數(shù)。代碼：(使用了內(nèi)聯(lián)函數(shù)，也可以如前例將函數(shù)定義為m文件)x=[234579121417212856];r=rand(1,12)-0.5;y1=[354247535965687376828699];y=y1+r;myfunc=inline(lbeta(1)+beta(2)*exp(beta(4)*x)+beta(3)*exp(-beta(4)*x)','beta','x');[beta,r,J]=nlinfit(x,y,myfunc,[0.50.50.50.5])[y_fit,delta]=nlpredci(myfunc,xi,beta,r,J)plot(x,y,|Q|,x,yfit,iri)運(yùn)行結(jié)果：beta= 87.5790 0.0283-62.7214 0.1070r=-1.76570.13760.03782.2882 1.17550.9553 -2.6909-0.6823 -1.1195 0.6254 1.0755 -0.0368J=1.00001.23860.8074101.34891.00001.37840.7255136.62131.00001.53410.6518163.71261.00001.70730.5857183.92401.00002.11470.4729208.03871.00002.61920.3818216.18661.00003.61050.2770209.69251.00004.47190.2236198.13391.00006.16430.1622175.9447

1.00009.45660.1057144.91241.00009.45660.1057144.91241.0000 19.9978 0.0500103.69391.0000399.9116 0.0025643.7199y_fi