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

關于最小二乘法線性和非線性擬合1第1頁，課件共54頁，創(chuàng)作于2023年2月2實驗目的實驗內容2、掌握用數(shù)學軟件求解擬合問題。1、直觀了解擬合基本內容。1、擬合問題引例及基本理論。4、實驗作業(yè)。2、用數(shù)學軟件求解擬合問題。3、應用實例第2頁，課件共54頁，創(chuàng)作于2023年2月3擬合2.擬合的基本原理1.擬合問題引例第3頁，課件共54頁，創(chuàng)作于2023年2月4擬合問題引例1溫度t(0C)20.532.751.073.095.7電阻R()7658268739421032已知熱敏電阻數(shù)據(jù)：求600C時的電阻R。

設

R=at+ba,b為待定系數(shù)第4頁，課件共54頁，創(chuàng)作于2023年2月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ù)坐標系(semilogy)下的圖形MATLAB(aa1)第5頁，課件共54頁，創(chuàng)作于2023年2月6曲線擬合問題的提法已知一組（二維）數(shù)據(jù)，即平面上n個點（xi,yi)i=1,…n,

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

使f(x)

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

+++++++++xyy=f(x)(xi,yi)ii為點（xi,yi)與曲線y=f(x)的距離第6頁，課件共54頁，創(chuàng)作于2023年2月7擬合與插值的關系

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

實例：下面數(shù)據(jù)是某次實驗所得，希望得到X和f之間的關系？MATLAB(cn)問題：給定一批數(shù)據(jù)點，需確定滿足特定要求的曲線或曲面解決方案：若不要求曲線（面）通過所有數(shù)據(jù)點，而是要求它反映對象整體的變化趨勢，這就是數(shù)據(jù)擬合，又稱曲線擬合或曲面擬合。若要求所求曲線（面）通過所給所有數(shù)據(jù)點，就是插值問題；第7頁，課件共54頁，創(chuàng)作于2023年2月8最臨近插值、線性插值、樣條插值與曲線擬合結果：第8頁，課件共54頁，創(chuàng)作于2023年2月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

的準則（最小二乘準則）：使n個點（xi,yi)與曲線y=f(x)的距離i的平方和最小。記

問題歸結為，求

a1,a2,…am

使

J(a1,a2,…am)

最小。第9頁，課件共54頁，創(chuàng)作于2023年2月10線性最小二乘法的求解：預備知識超定方程組：方程個數(shù)大于未知量個數(shù)的方程組即Ra=y其中超定方程一般是不存在解的矛盾方程組。

如果有向量a使得達到最小，則稱a為上述超定方程的最小二乘解。第10頁，課件共54頁，創(chuàng)作于2023年2月11線性最小二乘法的求解

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

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

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

1.通過機理分析建立數(shù)學模型來確定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頁，創(chuàng)作于2023年2月13用MATLAB解擬合問題1、線性最小二乘擬合2、非線性最小二乘擬合第13頁，課件共54頁，創(chuàng)作于2023年2月14用MATLAB作線性最小二乘擬合1.作多項式f(x)=a1xm+…+amx+am+1擬合,可利用已有程序:a=polyfit(x,y,m)2.對超定方程組可得最小二乘意義下的解。，用3.多項式在x處的值y可用以下命令計算：

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

am+1](數(shù)組)）輸入同長度的數(shù)組X，Y擬合多項式次數(shù)第14頁，課件共54頁，創(chuàng)作于2023年2月15即要求出二次多項式:中的使得:例對下面一組數(shù)據(jù)作二次多項式擬合第15頁，課件共54頁，創(chuàng)作于2023年2月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）計算結果：Ａ=-9.810820.1293-0.0317第16頁，課件共54頁，創(chuàng)作于2023年2月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ù)點和擬合曲線的圖形2）計算結果：Ａ=-9.810820.1293-0.0317解法2．用多項式擬合的命令MATLAB(zxec2)第17頁，課件共54頁，創(chuàng)作于2023年2月181.lsqcurvefit已知數(shù)據(jù)點：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頁，創(chuàng)作于2023年2月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ù)點選項見無約束優(yōu)化第19頁，課件共54頁，創(chuàng)作于2023年2月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ù)點：xdata=（xdata1，xdata2，…，xdatan）

ydata=（ydata1，ydata2，…，ydatan）第20頁，課件共54頁，創(chuàng)作于2023年2月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迭代初值選項見無約束優(yōu)化第21頁，課件共54頁，創(chuàng)作于2023年2月22

例2用下面一組數(shù)據(jù)擬合中的參數(shù)a，b，k該問題即解最優(yōu)化問題：第22頁，課件共54頁，創(chuàng)作于2023年2月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頁，創(chuàng)作于2023年2月243）運算結果為：f=0.00430.00510.00560.00590.00610.00620.00620.00630.00630.0063x=0.0063-0.00340.25424）結論：a=0.0063,b=-0.0034,k=0.2542第24頁，課件共54頁，創(chuàng)作于2023年2月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ù)，故應將cdatatdata的值寫在curvefun2.m中第25頁，課件共54頁，創(chuàng)作于2023年2月263）運算結果為

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可以看出,兩個命令的計算結果是相同的.4）結論：即擬合得a=0.0063b=-0.0034k=0.2542第26頁，課件共54頁，創(chuàng)作于2023年2月27MATLAB解應用問題實例1、電阻問題2、給藥方案問題*3、水塔流量估計問題第27頁，課件共54頁，創(chuàng)作于2023年2月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.直接用結果相同。MATLAB(dianzu2)第28頁，課件共54頁，創(chuàng)作于2023年2月29一室模型：將整個機體看作一個房室，稱中心室，室內血藥濃度是均勻的。快速靜脈注射后，濃度立即上升；然后迅速下降。當濃度太低時，達不到預期的治療效果；當濃度太高，又可能導致藥物中毒或副作用太強。臨床上，每種藥物有一個最小有效濃度c1和一個最大有效濃度c2。設計給藥方案時，要使血藥濃度保持在c1~c2之間。本題設c1=10，c2=25(ug/ml).擬合問題實例2給藥方案——一種新藥用于臨床之前，必須設計給藥方案.

藥物進入機體后血液輸送到全身，在這個過程中不斷地被吸收、分布、代謝，最終排出體外，藥物在血液中的濃度，即單位體積血液中的藥物含量，稱為血藥濃度。第29頁，課件共54頁，創(chuàng)作于2023年2月30

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

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

要設計給藥方案,必須知道給藥后血藥濃度隨時間變化的規(guī)律。從實驗和理論兩方面著手：第30頁，課件共54頁，創(chuàng)作于2023年2月給藥方案1.在快速靜脈注射的給藥方式下，研究血藥濃度（單位體積血液中的藥物含量）的變化規(guī)律。tc2cc10問題2.給定藥物的最小有效濃度和最大治療濃度，設計給藥方案：每次注射劑量多大；間隔時間多長。分析

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

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

在此，d=300mg，t及c（t）在某些點處的值見前表，需經(jīng)擬合求出參數(shù)k、v第32頁，課件共54頁，創(chuàng)作于2023年2月用線性最小二乘擬合c(t)MATLAB(lihe1)計算結果：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頁，創(chuàng)作于2023年2月給藥方案設計cc2c10t

設每次注射劑量D,間隔時間

血藥濃度c(t)

應c1c(t)c2

初次劑量D0應加大給藥方案記為：2、1、計算結果：給藥方案：c1=10,c2=25k=0.2347v=15.02第34頁，課件共54頁，創(chuàng)作于2023年2月35故可制定給藥方案：即:

首次注射375mg，其余每次注射225mg，注射的間隔時間為4小時。第35頁，課件共54頁，創(chuàng)作于2023年2月36估計水塔的流量2、解題思路3、算法設計與編程1、問題第36頁，課件共54頁，創(chuàng)作于2023年2月37

某居民區(qū)有一供居民用水的園柱形水塔，一般可以通過測量其水位來估計水的流量，但面臨的困難是，當水塔水位下降到設定的最低水位時，水泵自動啟動向水塔供水，到設定的最高水位時停止供水，這段時間無法測量水塔的水位和水泵的供水量．通常水泵每天供水一兩次，每次約兩小時.水塔是一個高12.2米，直徑17.4米的正園柱．按照設計，水塔水位降至約8.2米時，水泵自動啟動，水位升到約10.8米時水泵停止工作．表1是某一天的水位測量記錄，試估計任何時刻（包括水泵正供水時）從水塔流出的水流量，及一天的總用水量．第37頁，課件共54頁，創(chuàng)作于2023年2月38第38頁，課件共54頁，創(chuàng)作于2023年2月39流量估計的解題思路擬合水位~時間函數(shù)確定流量~時間函數(shù)估計一天總用水量第39頁，課件共54頁，創(chuàng)作于2023年2月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ù)直接分別作多項式擬合，得到水位函數(shù)．為使擬合曲線比較光滑，多項式次數(shù)不要太高，一般在3~6．由于第3時段只有3個測量記錄，無法對這一時段的水位作出較好的擬合．第40頁，課件共54頁，創(chuàng)作于2023年2月412、確定流量~時間函數(shù)

對于第1、2時段只需將水位函數(shù)求導數(shù)即可，對于兩個供水時段的流量，則用供水時段前后（水泵不工作時段）的流量擬合得到，并且將擬合得到的第2供水時段流量外推，將第3時段流量包含在第2供水時段內．第41頁，課件共54頁，創(chuàng)作于2023年2月423、一天總用水量的估計

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

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

%a1輸出多項式（系數(shù)為c1）導數(shù)的系數(shù)

3）tp1=0：0.1：9；

x1=-polyval（a1，tp1）；

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

MATLAB(llgj1)4）流量函數(shù)為：第44頁，課件共54頁，創(chuàng)作于2023年2月452、擬合第2時段的水位，并導出流量設t，h為已輸入的時刻和水位測量記錄（水泵啟動的4個時刻不輸入），第2時段各時刻的流量可如下得：1）c2=polyfit(t(10.9:21),h(10.9:21),3)；

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

%a2輸出多項式（系數(shù)為c2）導數(shù)的系數(shù)

3）tp2=10.9:0.1:21;x2=-polyval(a2,tp2);%x2輸出多項式（系數(shù)為a2）在tp2點的函數(shù)值（取負后邊為正值），即tp2時刻的流量MATLAB(llgj2)4）流量函數(shù)為：第45頁，課件共54頁，創(chuàng)作于2023年2月46

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

tp12=9：0.1：11；

x12=polyval（c12，tp12）；%x12輸出第1供水時段各時刻的流量MATLAB(llgj3)擬合的流量函數(shù)為：第46頁，課件共54頁，創(chuàng)作于2023年2月47

在第2供水時段之前取t=20，20.8兩點的流水量，在該時刻之后（第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次多項式

t3=20.8：0.1：24；

x3=polyval（c3，tp3）；%x3輸出第2供水時段（外推至t=24）各時刻的流量MATLAB(llgj4)擬合的流量函數(shù)為：第47頁，課件共54頁，創(chuàng)作于2023年2月483、一天總用水量的估計第1、2時段和第1、2供水時段流量的積分之和，就是一天總用水量．雖然諸時段的流量已表為多項式函數(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