最小二乘法線性和非線性擬合

關(guān)于最小二乘法線性和非線性擬合12實(shí)驗(yàn)?zāi)康膶?shí)驗(yàn)內(nèi)容2、掌握用數(shù)學(xué)軟件求解擬合問題。1、直觀了解擬合基本內(nèi)容。1、擬合問題引例及基本理論。4、實(shí)驗(yàn)作業(yè)。2、用數(shù)學(xué)軟件求解擬合問題。3、應(yīng)用實(shí)例第2頁,共54頁，2024年2月25日，星期天3擬合2.擬合的基本原理1.擬合問題引例第3頁,共54頁，2024年2月25日，星期天4擬合問題引例1溫度t(0C)20.532.751.073.095.7電阻R()7658268739421032已知熱敏電阻數(shù)據(jù)：求600C時的電阻R。

設(shè)

R=at+ba,b為待定系數(shù)第4頁,共54頁，2024年2月25日，星期天5擬合問題引例2

t(h)0.250.511.523468c(g/ml)19.2118.1515.3614.1012.899.327.455.243.01已知一室模型快速靜脈注射下的血藥濃度數(shù)據(jù)(t=0注射300mg)求血藥濃度隨時間的變化規(guī)律c(t).作半對數(shù)坐標(biāo)系(semilogy)下的圖形MATLAB(aa1)第5頁,共54頁，2024年2月25日，星期天6曲線擬合問題的提法已知一組（二維）數(shù)據(jù)，即平面上n個點(diǎn)（xi,yi)i=1,…n,

尋求一個函數(shù)（曲線）y=f(x),

使f(x)

在某種準(zhǔn)則下與所有數(shù)據(jù)點(diǎn)最為接近，即曲線擬合得最好。

+++++++++xyy=f(x)(xi,yi)

i

i為點(diǎn)（xi,yi)與曲線y=f(x)的距離第6頁,共54頁，2024年2月25日，星期天7擬合與插值的關(guān)系

函數(shù)插值與曲線擬合都是要根據(jù)一組數(shù)據(jù)構(gòu)造一個函數(shù)作為近似，由于近似的要求不同，二者的數(shù)學(xué)方法上是完全不同的。

實(shí)例：下面數(shù)據(jù)是某次實(shí)驗(yàn)所得，希望得到X和f之間的關(guān)系？MATLAB(cn)問題：給定一批數(shù)據(jù)點(diǎn)，需確定滿足特定要求的曲線或曲面解決方案：若不要求曲線（面）通過所有數(shù)據(jù)點(diǎn)，而是要求它反映對象整體的變化趨勢，這就是數(shù)據(jù)擬合，又稱曲線擬合或曲面擬合。若要求所求曲線（面）通過所給所有數(shù)據(jù)點(diǎn)，就是插值問題；第7頁,共54頁，2024年2月25日，星期天8最臨近插值、線性插值、樣條插值與曲線擬合結(jié)果：第8頁,共54頁，2024年2月25日，星期天9曲線擬合問題最常用的解法——線性最小二乘法的基本思路第一步:先選定一組函數(shù)

r1(x),r2(x),…rm(x),m<n,

令

f(x)=a1r1(x)+a2r2(x)+…+amrm(x)（1）其中

a1,a2,…am

為待定系數(shù)。

第二步:確定a1,a2,…am

的準(zhǔn)則（最小二乘準(zhǔn)則）：使n個點(diǎn)（xi,yi)與曲線y=f(x)的距離

i的平方和最小。記

問題歸結(jié)為，求

a1,a2,…am

使

J(a1,a2,…am)

最小。第9頁,共54頁，2024年2月25日，星期天10線性最小二乘法的求解：預(yù)備知識超定方程組：方程個數(shù)大于未知量個數(shù)的方程組即Ra=y其中超定方程一般是不存在解的矛盾方程組。

如果有向量a使得達(dá)到最小，則稱a為上述超定方程的最小二乘解。第10頁,共54頁，2024年2月25日，星期天11線性最小二乘法的求解

定理：當(dāng)RTR可逆時，超定方程組（3）存在最小二乘解，且即為方程組

RTRa=RTy的解：a=(RTR)-1RTy

所以，曲線擬合的最小二乘法要解決的問題，實(shí)際上就是求以下超定方程組的最小二乘解的問題。其中Ra=y（3）第11頁,共54頁，2024年2月25日，星期天12線性最小二乘擬合f(x)=a1r1(x)+…+amrm(x)中函數(shù){r1(x),…rm(x)}的選取

1.通過機(jī)理分析建立數(shù)學(xué)模型來確定f(x)；++++++++++++++++++++++++++++++f=a1+a2xf=a1+a2x+a3x2f=a1+a2x+a3x2f=a1+a2/xf=aebxf=ae-bx2.將數(shù)據(jù)(xi,yi)i=1,…n作圖，通過直觀判斷確定f(x)：第12頁,共54頁，2024年2月25日，星期天13用MATLAB解擬合問題1、線性最小二乘擬合2、非線性最小二乘擬合第13頁,共54頁，2024年2月25日，星期天14用MATLAB作線性最小二乘擬合1.作多項(xiàng)式f(x)=a1xm+…+amx+am+1擬合,可利用已有程序:a=polyfit(x,y,m)2.對超定方程組可得最小二乘意義下的解。，用3.多項(xiàng)式在x處的值y可用以下命令計算：

y=polyval（a，x）輸出擬合多項(xiàng)式系數(shù)a=[a1,…am,

am+1](數(shù)組)）輸入同長度的數(shù)組X，Y擬合多項(xiàng)式次數(shù)第14頁,共54頁，2024年2月25日，星期天15即要求出二次多項(xiàng)式:中的使得:例對下面一組數(shù)據(jù)作二次多項(xiàng)式擬合第15頁,共54頁，2024年2月25日，星期天161）輸入以下命令：x=0:0.1:1;y=[-0.4471.9783.286.167.087.347.669.569.489.3011.2];R=[(x.^2)'x'ones(11,1)]；

A=R\y'MATLAB(zxec1)解法1．用解超定方程的方法2）計算結(jié)果：Ａ=-9.810820.1293-0.0317第16頁,共54頁，2024年2月25日，星期天171）輸入以下命令：

x=0:0.1:1;y=[-0.4471.9783.286.167.087.347.669.569.489.3011.2];A=polyfit(x,y,2)z=polyval(A,x);plot(x,y,'k+',x,z,'r')%作出數(shù)據(jù)點(diǎn)和擬合曲線的圖形2）計算結(jié)果：Ａ=-9.810820.1293-0.0317解法2．用多項(xiàng)式擬合的命令MATLAB(zxec2)第17頁,共54頁，2024年2月25日，星期天181.lsqcurvefit已知數(shù)據(jù)點(diǎn)：xdata=（xdata1，xdata2，…，xdatan），

ydata=（ydata1，ydata2，…，ydatan）

用MATLAB作非線性最小二乘擬合Matlab的提供了兩個求非線性最小二乘擬合的函數(shù)：lsqcurvefit和lsqnonlin。兩個命令都要先建立M-文件fun.m，在其中定義函數(shù)f(x)，但兩者定義f(x)的方式是不同的,可參考例題.

lsqcurvefit用以求含參量x（向量）的向量值函數(shù)F(x,xdata)=（F（x，xdata1），…，F(xiàn)（x，xdatan））T中的參變量x(向量),使得

第18頁,共54頁，2024年2月25日，星期天19

輸入格式為:(1)x=lsqcurvefit(‘fun’,x0,xdata,ydata);(2)x=lsqcurvefit(‘fun’,x0,xdata,ydata,options);(3)x=lsqcurvefit(‘fun’,x0,xdata,ydata,options,’grad’);(4)[x,options]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);(5)[x,options,funval]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);(6)[x,options,funval,Jacob]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);fun是一個事先建立的定義函數(shù)F(x,xdata)的M-文件,自變量為x和xdata說明：x=lsqcurvefit(‘fun’,x0,xdata,ydata,options);迭代初值已知數(shù)據(jù)點(diǎn)選項(xiàng)見無約束優(yōu)化第19頁,共54頁，2024年2月25日，星期天20

lsqnonlin用以求含參量x（向量）的向量值函數(shù)

f(x)=(f1(x),f2(x),…,fn(x))T

中的參量x，使得

最小。其中fi（x）=f（x，xdatai，ydatai）

=F(x,xdatai)-ydatai

2.lsqnonlin已知數(shù)據(jù)點(diǎn)：xdata=（xdata1，xdata2，…，xdatan）

ydata=（ydata1，ydata2，…，ydatan）第20頁,共54頁，2024年2月25日，星期天21輸入格式為：

1）x=lsqnonlin（‘fun’，x0）；

2）x=lsqnonlin

（‘fun’，x0，options）；

3）x=lsqnonlin

（‘fun’，x0，options，‘grad’）；

4）[x，options]=lsqnonlin

（‘fun’，x0，…）；

5）[x，options，funval]=lsqnonlin

（‘fun’，x0，…）；說明：x=lsqnonlin

（‘fun’，x0，options）；fun是一個事先建立的定義函數(shù)f(x)的M-文件，自變量為x迭代初值選項(xiàng)見無約束優(yōu)化第21頁,共54頁，2024年2月25日，星期天22

例2用下面一組數(shù)據(jù)擬合中的參數(shù)a，b，k該問題即解最優(yōu)化問題：第22頁,共54頁，2024年2月25日，星期天23MATLAB(fzxec1)

1）編寫M-文件curvefun1.m

functionf=curvefun1(x,tdata)f=x(1)+x(2)*exp(-0.02*x(3)*tdata)%其中x(1)=a;x(2)=b；x(3)=k;2）輸入命令tdata=100:100:1000cdata=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=lsqcurvefit('curvefun1',x0,tdata,cdata)f=curvefun1(x,tdata)

F(x，tdata)=，x=(a，b，k)解法1.用命令lsqcurvefit第23頁,共54頁，2024年2月25日，星期天243）運(yùn)算結(jié)果為：f=0.00430.00510.00560.00590.00610.00620.00620.00630.00630.0063x=0.0063-0.00340.25424）結(jié)論：a=0.0063,b=-0.0034,k=0.2542第24頁,共54頁，2024年2月25日，星期天25MATLAB(fzxec2)

解法2

用命令lsqnonlin

f(x)=F(x,tdata,ctada)=x=（a，b，k）1）編寫M-文件curvefun2.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)-cdata2）輸入命令:

x0=[0.2,0.05,0.05];x=lsqnonlin('curvefun2',x0)f=curvefun2(x)函數(shù)curvefun2的自變量是x，cdata和tdata是已知參數(shù)，故應(yīng)將cdatatdata的值寫在curvefun2.m中第25頁,共54頁，2024年2月25日，星期天263）運(yùn)算結(jié)果為

f=1.0e-003*(0.2322-0.1243-0.2495-0.2413-0.1668-0.07240.02410.11590.20300.2792x=0.0063-0.00340.2542可以看出,兩個命令的計算結(jié)果是相同的.4）結(jié)論：即擬合得a=0.0063b=-0.0034k=0.2542第26頁,共54頁，2024年2月25日，星期天27MATLAB解應(yīng)用問題實(shí)例1、電阻問題2、給藥方案問題*3、水塔流量估計問題第27頁,共54頁，2024年2月25日，星期天28MATLAB(dianzu1)電阻問題溫度t(0C)20.532.751.073.095.7電阻R()7658268739421032例.由數(shù)據(jù)擬合R=a1t+a2方法1.用命令polyfit(x,y,m)得到a1=3.3940,a2=702.4918方法2.直接用結(jié)果相同。MATLAB(dianzu2)第28頁,共54頁，2024年2月25日，星期天29一室模型：將整個機(jī)體看作一個房室，稱中心室，室內(nèi)血藥濃度是均勻的。快速靜脈注射后，濃度立即上升；然后迅速下降。當(dāng)濃度太低時，達(dá)不到預(yù)期的治療效果；當(dāng)濃度太高，又可能導(dǎo)致藥物中毒或副作用太強(qiáng)。臨床上，每種藥物有一個最小有效濃度c1和一個最大有效濃度c2。設(shè)計給藥方案時，要使血藥濃度保持在c1~c2之間。本題設(shè)c1=10，c2=25(ug/ml).擬合問題實(shí)例2給藥方案——一種新藥用于臨床之前，必須設(shè)計給藥方案.

藥物進(jìn)入機(jī)體后血液輸送到全身，在這個過程中不斷地被吸收、分布、代謝，最終排出體外，藥物在血液中的濃度，即單位體積血液中的藥物含量，稱為血藥濃度。第29頁,共54頁，2024年2月25日，星期天30

在實(shí)驗(yàn)方面,對某人用快速靜脈注射方式一次注入該藥物300mg后,在一定時刻t(小時)采集血藥,測得血藥濃度c(ug/ml)如下表:

t(h)0.250.511.523468c(g/ml)19.2118.1515.3614.1012.899.327.455.243.01

要設(shè)計給藥方案,必須知道給藥后血藥濃度隨時間變化的規(guī)律。從實(shí)驗(yàn)和理論兩方面著手：第30頁,共54頁，2024年2月25日，星期天給藥方案1.在快速靜脈注射的給藥方式下，研究血藥濃度（單位體積血液中的藥物含量）的變化規(guī)律。tc2cc10

問題2.給定藥物的最小有效濃度和最大治療濃度，設(shè)計給藥方案：每次注射劑量多大；間隔時間多長。分析

理論：用一室模型研究血藥濃度變化規(guī)律

實(shí)驗(yàn)：對血藥濃度數(shù)據(jù)作擬合，符合負(fù)指數(shù)變化規(guī)律第31頁,共54頁，2024年2月25日，星期天3.血液容積v,t=0注射劑量d,血藥濃度立即為d/v.2.藥物排除速率與血藥濃度成正比，比例系數(shù)k(>0)模型假設(shè)1.機(jī)體看作一個房室，室內(nèi)血藥濃度均勻——一室模型模型建立

在此，d=300mg，t及c（t）在某些點(diǎn)處的值見前表，需經(jīng)擬合求出參數(shù)k、v第32頁,共54頁，2024年2月25日，星期天用線性最小二乘擬合c(t)MATLAB(lihe1)計算結(jié)果：d=300;t=[0.250.511.523468];c=[19.2118.1515.3614.1012.899.327.455.243.01];y=log(c);a=polyfit(t,y,1)k=-a(1)v=d/exp(a(2))程序：用非線性最小二乘擬合c(t)第33頁,共54頁，2024年2月25日，星期天給藥方案設(shè)計cc2c10

t

設(shè)每次注射劑量D,間隔時間

血藥濃度c(t)

應(yīng)c1

c(t)

c2

初次劑量D0應(yīng)加大給藥方案記為：2、1、計算結(jié)果：給藥方案：c1=10,c2=25k=0.2347v=15.02第34頁,共54頁，2024年2月25日，星期天35故可制定給藥方案：即:

首次注射375mg，其余每次注射225mg，注射的間隔時間為4小時。第35頁,共54頁，2024年2月25日，星期天36估計水塔的流量2、解題思路3、算法設(shè)計與編程1、問題第36頁,共54頁，2024年2月25日，星期天37

某居民區(qū)有一供居民用水的園柱形水塔，一般可以通過測量其水位來估計水的流量，但面臨的困難是，當(dāng)水塔水位下降到設(shè)定的最低水位時，水泵自動啟動向水塔供水，到設(shè)定的最高水位時停止供水，這段時間無法測量水塔的水位和水泵的供水量．通常水泵每天供水一兩次，每次約兩小時.水塔是一個高12.2米，直徑17.4米的正園柱．按照設(shè)計，水塔水位降至約8.2米時，水泵自動啟動，水位升到約10.8米時水泵停止工作．表1是某一天的水位測量記錄，試估計任何時刻（包括水泵正供水時）從水塔流出的水流量，及一天的總用水量．第37頁,共54頁，2024年2月25日，星期天38第38頁,共54頁，2024年2月25日，星期天39流量估計的解題思路擬合水位~時間函數(shù)確定流量~時間函數(shù)估計一天總用水量第39頁,共54頁，2024年2月25日，星期天40

擬合水位~時間函數(shù)測量記錄看，一天有兩個供水時段（以下稱第1供水時段和第2供水時段），和3個水泵不工作時段（以下稱第1時段t=0到t=8.97，第2次時段t=10.95到t=20.84和第3時段t=23以后）．對第1、2時段的測量數(shù)據(jù)直接分別作多項(xiàng)式擬合，得到水位函數(shù)．為使擬合曲線比較光滑，多項(xiàng)式次數(shù)不要太高，一般在3~6．由于第3時段只有3個測量記錄，無法對這一時段的水位作出較好的擬合．第40頁,共54頁，2024年2月25日，星期天412、確定流量~時間函數(shù)

對于第1、2時段只需將水位函數(shù)求導(dǎo)數(shù)即可，對于兩個供水時段的流量，則用供水時段前后（水泵不工作時段）的流量擬合得到，并且將擬合得到的第2供水時段流量外推，將第3時段流量包含在第2供水時段內(nèi)．第41頁,共54頁，2024年2月25日，星期天423、一天總用水量的估計

總用水量等于兩個水泵不工作時段和兩個供水時段用水量之和，它們都可以由流量對時間的積分得到。第42頁,共54頁，2024年2月25日，星期天43算法設(shè)計與編程1、擬合第1、2時段的水位，并導(dǎo)出流量2、擬合供水時段的流量3、估計一天總用水量4、流量及總用水量的檢驗(yàn)第43頁,共54頁，2024年2月25日，星期天441、擬合第1時段的水位，并導(dǎo)出流量設(shè)t，h為已輸入的時刻和水位測量記錄（水泵啟動的4個時刻不輸入），第1時段各時刻的流量可如下得：1）c1=polyfit（t（1：10），h（1：10），3）；

%用3次多項(xiàng)式擬合第1時段水位，c1輸出3次多項(xiàng)式的系數(shù)2）a1=polyder（c1）；

%a1輸出多項(xiàng)式（系數(shù)為c1）導(dǎo)數(shù)的系數(shù)

3）tp1=0：0.1：9；

x1=-polyval（a1，tp1）；

%x1輸出多項(xiàng)式（系數(shù)為a1）在tp1點(diǎn)的函數(shù)值（取負(fù)后邊為正值），即tp1時刻的流量

MATLAB(llgj1)4）流量函數(shù)為：第44頁,共54頁，2024年2月25日，星期天452、擬合第2時段的水位，并導(dǎo)出流量設(shè)t，h為已輸入的時刻和水位測量記錄（水泵啟動的4個時刻不輸入），第2時段各時刻的流量可如下得：1）c2=polyfit(t(10.9:21),h(10.9:21),3)；

%用3次多項(xiàng)式擬合第2時段水位，c2輸出3次多項(xiàng)式的系數(shù)2）a2=polyder(c2)；

%a2輸出多項(xiàng)式（系數(shù)為c2）導(dǎo)數(shù)的系數(shù)

3）tp2=10.9:0.1:21;x2=-polyval(a2,tp2);%x2輸出多項(xiàng)式（系數(shù)為a2）在tp2點(diǎn)的函數(shù)值（取負(fù)后邊為正值），即tp2時刻的流量MATLAB(llgj2)4）流量函數(shù)為：第45頁,共54頁，2024年2月25日，星期天46

3、擬合供水時段的流量在第1供水時段（t=9~11）之前（即第1時段）和之后（即第2時段）各取幾點(diǎn)，其流量已經(jīng)得到，用它們擬合第1供水時段的流量．為使流量函數(shù)在t=9和t=11連續(xù)，我們簡單地只取4個點(diǎn)，擬合3次多項(xiàng)式（即曲線必過這4個點(diǎn)），實(shí)現(xiàn)如下：

xx1=-polyval（a1，[89]）；%取第1時段在t=8,9的流量

xx2=-polyval（a2，[1112]）；%取第2時段在t=11,12的流量

xx12=[xx1xx2]；

c12=polyfit（[891112]，xx12，3）；%擬合3次多項(xiàng)式

tp12=9：0.1：11；

x12=polyval（c12，tp12）；%x12輸出第1供水時段各時刻的流量MATLAB(llgj3)擬合的流量函數(shù)為：第46頁,共54頁，2024年2月25日，星期天47

在第2供水時段之前取t=20，20.8兩點(diǎn)的流水量，在該時刻之后（第3時段）僅有3個水位記錄，我們用差分得到流量，然后用這4個數(shù)值擬合第2供水時段的流量如下：

dt3=diff（t(22：24)）；%最后3個時刻的兩兩之差

dh3=diff（h(22：24)）；%最后3個水位的兩兩之差

dht3=-dh3./dt3；%t(22)和t(23)的流量

t3=[2020.8t(22)t(23)]；

xx3=[-polyval(a2，t3(1：2)，dht3)]；%取t3各時刻的流量

c3=polyfit（t3，xx3，3）；%擬合3次多項(xiàng)式

t3=20.8：0.1：24；

x3=polyval（c3，tp3）；%x3輸出第2供水時段（外推至t=24）各時刻的流量MATLAB(llgj4)擬合的流量函數(shù)為：第47頁,共54頁，2024年2月25日，星期天483、一天總用水量的估計第1、2時段和第1、2供水時段流量的積分之和，就是一天總用水量．雖然諸時段的流量已表為多項(xiàng)式函數(shù)，積分可以解析地算出，這里仍用數(shù)值積分計算如下：

y1=0.1*trapz(x1)；%第1時段用水量（仍按高度計），0.1為積分步長

y2=0.1*trapz(x2)；%第2時段用水量

y12=0.1*trapz(x12)；%第1供水時段用水量

y3=0.1*trapz(x3)；%第2供水時段用水量

y=（y1+y2+y12+y3）*237.8*0.01；%一天總用