第4講擬合與回歸

擬合與回歸四川師范大學(xué)趙凌曲線擬合已知一組（二維）數(shù)據(jù)，即平面上n個(gè)點(diǎn)（xi,yi)i=1,…n,尋求一個(gè)函數(shù)（曲線）y=f(x),使f(x)在某種準(zhǔn)則下與所有數(shù)據(jù)點(diǎn)最為接近，即曲線擬合得最好。

+++++++++xy一、擬合實(shí)例：下面數(shù)據(jù)是某次實(shí)驗(yàn)所得，希望得到X和f之間的關(guān)系？問題：給定一批數(shù)據(jù)點(diǎn)，需確定滿足特定要求的曲線或曲面解決方案：若不要求曲線（面）通過所有數(shù)據(jù)點(diǎn)，而是要求它反映對(duì)象整體的變化趨勢(shì)，這就是數(shù)據(jù)擬合，又稱曲線擬合或曲面擬合。若要求所求曲線（面）通過所給所有數(shù)據(jù)點(diǎn)，是插值問題；根據(jù)散點(diǎn)圖或者經(jīng)驗(yàn)公式，確定函數(shù)的形式。

函數(shù)的形式分成兩種：線性的（可化為線性的）和非線性的可化為線性的：確定一組函數(shù)r1(x),r2(x),…rm(x),設(shè)y=a1r1(x)+a2r2(x)+…+amrm(x)其中a1,a2,…am

為待定系數(shù)。線性的：y=a1x1+a2x2+…+amxm,要求m<n(樣本容量）

其中a1,a2,…am

為待定系數(shù)。非線性的：y=f(x)通常選擇的六類曲線如下：(7)多項(xiàng)式線性化方法兩端取對(duì)數(shù)得：lny

=ln

+

x令：y'=lny，則有y'

=ln

+

x基本形式：圖像

<指數(shù)函數(shù)線性化方法兩端取對(duì)數(shù)得：lny

=ln

+/x令：y'

=lny，x'

=1/x，則有y'

=ln

+

x'基本形式：圖像負(fù)指數(shù)函數(shù)

<線性化方法兩端取對(duì)數(shù)得：lg

y=lg

+

lgx令：y'=lgy，x'=lg

x，則y'

=lg

+x'基本形式：圖像0<<1

1

=1-1<

<0

<-1

=-1冪函數(shù)線性化方法令：y'=1/y，x'=1/x,則有y'

=

+

x'基本形式：圖像

<0

>0雙曲線函數(shù)線性化方法x'=lgx

,則有y'

=

+

x'基本形式：圖像

0

<0對(duì)數(shù)函數(shù)線性化方法令：y'=1/y，x'=e-x,則有y'

=

+

x'基本形式：圖像

S

型曲線線性化方法令：y'=y，x1=x,x2=x2,…,xn=xn則有基本形式：圖像多項(xiàng)式曲線polyfit

polyval

對(duì)于不能化為線性模型的非線性模型，應(yīng)直接用非線性最小二乘法處理曲線擬合的Matlab實(shí)現(xiàn)

b=regress(Y,X)1、確定回歸系數(shù)的點(diǎn)估計(jì)值：線性回歸：regress非線性回歸函數(shù)nlinfit1.

利用nlinfit函數(shù)作非線性擬合調(diào)用格式未知參數(shù)事先用m-文件定義的非線性函數(shù)[beta,r,J,COVB,mse]=nlinfit(X,y,fun,b0,options)回歸系數(shù)初值優(yōu)化屬性設(shè)置雅可比矩陣殘差lsqnonlin，lsqcurvefit此處有鏈接案例

根據(jù)經(jīng)驗(yàn)，人口增長的預(yù)測(cè)模型通常采用Logistic函數(shù)其中y(t)為t時(shí)刻人口數(shù)，A，B，C為常數(shù)。試根據(jù)1975-2005年的中國人口數(shù)據(jù)（見下頁表），得出中國人口增長預(yù)測(cè)模型。年份時(shí)間人口（萬人）人口（億人）年份時(shí)間人口（萬人）人口（億人）19750924209.24219911611582311.582319761937179.371719921711717111.717119772949749.497419931811851711.851719783962599.625919941911985011.98519794975429.754219952012112112.112119805987059.870519962112238912.23891981610007210.007219972212362612.36261982710165410.165419982312476112.47611983810300810.300819992412578612.57861984910435710.435720002512674312.674319851010585110.585120012612762712.762719861110750710.750720022712845312.845319871210930010.9320032812922712.922719881311102611.102620042912998812.998819891411270411.270420053013075613.075619901511433311.4333原始數(shù)據(jù)散點(diǎn)和折線圖year=renkou_data(:,1);t=renkou_data(:,2);y=renkou_data(:,4);fun=@(beta,t)[beta(1)./(1+beta(2)*exp(beta(3)*t))];[beta,resid,J]=nlinfit(t,y,fun,[15,1,1]);beta=16.16340.7712-0.0408部分結(jié)果Logistic函數(shù)表達(dá)式為：建模案例酶促反應(yīng)

問題研究酶促反應(yīng)（酶催化反應(yīng)）中嘌呤霉素的反應(yīng)速度與底物（反應(yīng)物）濃度之間關(guān)系的影響

建立數(shù)學(xué)模型，反映該酶促反應(yīng)的速度與底物濃度以及經(jīng)嘌呤霉素處理與否之間的關(guān)系

設(shè)計(jì)了兩個(gè)實(shí)驗(yàn)：酶經(jīng)過嘌呤霉素處理；酶未經(jīng)嘌呤霉素處理。實(shí)驗(yàn)數(shù)據(jù)見下表:

方案底物濃度(ppm)0.020.060.110.220.561.10反應(yīng)速度處理764797107123139159152191201207200未處理6751848698115131124144158160/可選的模型有：Michaelis-Menten模型指數(shù)增長模型線性化模型

經(jīng)嘌呤霉素處理后實(shí)驗(yàn)數(shù)據(jù)的估計(jì)結(jié)果

參數(shù)參數(shù)估計(jì)值（×10-3）置信區(qū)間（×10-3）

15.107[3.5396.676]

20.247[0.1760.319]R2=0.8557F=59.2975p=0.0000對(duì)

1

,2非線性

對(duì)

1,

2線性

線性化模型結(jié)果分析

x較大時(shí)，y有較大偏差1/x較小時(shí)有很好的線性趨勢(shì)，1/x較大時(shí)出現(xiàn)很大的起落

參數(shù)估計(jì)時(shí)，x較?。?/x很大）的數(shù)據(jù)控制了回歸參數(shù)的確定

1/y1/xxyx=[0.020.020.060.060.110.110.220.220.560.561.101.10];y=[764797107123139159152191201207200];%建立自定義函數(shù)zl=@(beta,t)[beta(1)*t./(beta(2)+t)];%在命令窗口輸入命令bo=[195.80270.04841];[beta,R,J]=nlinfit(x,y,zl,bo')直接用非線性擬合公式nlinfit求參數(shù)估計(jì)混合反應(yīng)模型

在同一模型中考慮嘌呤霉素處理的影響x=[0.02001.00000.02001.00000.06001.00000.06001.00000.11001.00000.11001.00000.22001.00000.22001.00000.56001.00000.56001.00001.10001.00001.100000.020000.020000.060000.060000.110000.110000.220000.220000.560000.560001.100001.10000]y=[7647971071231391591521912012072006751848698115131124144158160162]functionY=zl2(beta,t)Y=(beta(1)*ones(length(t(:,1)),1)...+beta(3)*t(:,2)).*t(:,1)./...((beta(2)*ones(length(t(:,1)),1)+beta(4)*t(:,2))+t(:,1));bo=[1700.05600.01];[beta,R,j]=nlinfit(x,y,@zl2,bo');Nlintool(x,y,@zl2,bo)二、回歸曲線擬合問題的特點(diǎn)是，根據(jù)得到的若干有關(guān)變量的一組數(shù)據(jù)，尋找因變量與自變量之間的一個(gè)函數(shù)，使這個(gè)函數(shù)對(duì)那組數(shù)據(jù)擬合得最好。通常，函數(shù)的形式可以由經(jīng)驗(yàn)、先驗(yàn)知識(shí)或?qū)?shù)據(jù)的直觀觀察決定，要作的工作是由數(shù)據(jù)用最小二乘法計(jì)算函數(shù)中的待定系數(shù)。從計(jì)算的角度看，問題似乎已經(jīng)完全解決了，還有進(jìn)一步研究的必要嗎?回歸的步驟：一、確定函數(shù)形式（散點(diǎn)圖或先驗(yàn)知識(shí)）二、確定待定參數(shù)（曲線擬合）三、擬合優(yōu)度檢驗(yàn)四、顯著性檢驗(yàn)回歸分析就是對(duì)擬合問題作的統(tǒng)計(jì)分析。二、回歸32

第一節(jié)多元線性回歸模型

一、多元線性回歸模型

一般形式：對(duì)于有K-1個(gè)解釋變量的線性回歸模型

33

二、多元線性回歸模型的矩陣表示

多個(gè)解釋變量的多元線性回歸模型的n組樣本觀測(cè)值，可表示為

用矩陣表示

33

第二節(jié)多元線性回歸模型的估計(jì)

一、普通最小二乘法（OLS）原則：尋求剩余平方和最小的參數(shù)估計(jì)式

即求偏導(dǎo)，并令其為0其中即3435用矩陣表示的正規(guī)方程偏導(dǎo)數(shù)因?yàn)闃颖净貧w函數(shù)為兩邊左乘根據(jù)最小二乘原則則正規(guī)方程為回歸系數(shù)的區(qū)間估計(jì)不含原點(diǎn)，則此變量線性關(guān)性顯著回歸系數(shù)的點(diǎn)估計(jì)和區(qū)間估計(jì)、并檢驗(yàn)回歸模型：

[b,bint,r,rint,stats]=regress(Y,X,alpha)回歸系數(shù)的點(diǎn)估計(jì)殘差用于檢驗(yàn)回歸模型的統(tǒng)計(jì)量，有三個(gè)數(shù)值：相關(guān)系數(shù)r2、F值、與F對(duì)應(yīng)的概率p置信區(qū)間顯著性水平（缺省時(shí)為0.05）37第三節(jié)

多元線性回歸模型的檢驗(yàn)

一、多元回歸的擬合優(yōu)度檢驗(yàn)

多重可決系數(shù)：在多元回歸模型中，由各個(gè)解釋變量聯(lián)合起來解釋了的Y的變差，在Y的總變差中占的比重，用表示

多元回歸中多重可決系數(shù)可表示為

37回歸系數(shù)的區(qū)間估計(jì)不含原點(diǎn)，則此變量線性關(guān)性顯著回歸系數(shù)的點(diǎn)估計(jì)和區(qū)間估計(jì)、并檢驗(yàn)回歸模型：

[b,bint,r,rint,stats]=regress(Y,X,alpha)回歸系數(shù)的區(qū)間估計(jì)殘差用于檢驗(yàn)回歸模型的統(tǒng)計(jì)量，有三個(gè)數(shù)值：相關(guān)系數(shù)r2、F值、與F對(duì)應(yīng)的概率p置信區(qū)間顯著性水平（缺省時(shí)為0.05）線性關(guān)系檢驗(yàn)－F檢驗(yàn)檢驗(yàn)因變量與所有自變量之間的是否顯著也被稱為總體的顯著性檢驗(yàn)檢驗(yàn)方法是將回歸離差平方和(SSR)同剩余離差平方和(SSE)加以比較，應(yīng)用F檢驗(yàn)來分析二者之間的差別是否顯著如果是顯著的，因變量與自變量之間存在線性關(guān)系如果不顯著，因變量與自變量之間不存在線性關(guān)系線性關(guān)系檢驗(yàn)提出假設(shè)H0：

1

2

p=0線性關(guān)系不顯著H1：

1，

2，，

p至少有一個(gè)不等于02.計(jì)算檢驗(yàn)統(tǒng)計(jì)量F3.確定顯著性水平和分子自由度p、分母自由度n-p-1找出臨界值F

4.作出決策：若F>F

，拒絕H0回歸系數(shù)的區(qū)間估計(jì)不含原點(diǎn)，則此變量線性關(guān)性顯著回歸系數(shù)的點(diǎn)估計(jì)和區(qū)間估計(jì)、并檢驗(yàn)回歸模型：

[b,bint,r,rint,stats]=regress(Y,X,alpha)回歸系數(shù)的區(qū)間估計(jì)殘差用于檢驗(yàn)回歸模型的統(tǒng)計(jì)量，有三個(gè)數(shù)值：相關(guān)系數(shù)r2、F值、與F對(duì)應(yīng)的概率p置信區(qū)間顯著性水平（缺省時(shí)為0.05）回歸系數(shù)的檢驗(yàn)－t檢驗(yàn)線性關(guān)系檢驗(yàn)通過后，對(duì)各個(gè)回歸系數(shù)有選擇地進(jìn)行一次或多次檢驗(yàn)應(yīng)用t檢驗(yàn)統(tǒng)計(jì)量回歸系數(shù)的檢驗(yàn)

(步驟)提出假設(shè)H0：bi=0(自變量xi與

因變量y沒有線性關(guān)系)H1：bi

0(自變量xi與

因變量y有線性關(guān)系)計(jì)算檢驗(yàn)的統(tǒng)計(jì)量t確定顯著性水平，并進(jìn)行決策

t>t

，拒絕H0；t<t

，不拒絕H0回歸系數(shù)的區(qū)間估計(jì)不含原點(diǎn)，則此變量線性關(guān)性顯著2、求回歸系數(shù)的點(diǎn)估計(jì)和區(qū)間估計(jì)、并檢驗(yàn)回歸模型：

[b,bint,r,rint,stats]=regress(Y,X,alpha)回歸系數(shù)的區(qū)間估計(jì)殘差用于檢驗(yàn)回歸模型的統(tǒng)計(jì)量，有三個(gè)數(shù)值：相關(guān)系數(shù)r2、F值、與F對(duì)應(yīng)的概率p置信區(qū)間顯著性水平（缺省時(shí)為0.05）

3、案例

考察15名不同程度的煙民的每日抽煙量、飲酒量（啤酒）與其心電圖指標(biāo)(zb)的對(duì)應(yīng)數(shù)據(jù)，試建立心電圖指標(biāo)關(guān)于日抽煙量和日飲酒量的適合的回歸模型。組別(g)日抽煙量(x)/支日飲酒量(y)/升心電圖指標(biāo)(zb)130102801251126013513330140144001451441022012170218112102251228022513300223132903401441034515420348164253501845035519470xyz=[。。。。。。];x=xyz(:,1);y=xyz(:,2);z=xyz(:,3);n=size(x,1);xy=[ones(n,1),x,y];[b,bint,r,rint,stats]=regress(z,xy)b= r= rint= 66.0944 -17.7298 -63.672 28.21246.9774 -5.0743 -62.5486 52.40012.2314 -9.3109 -68.5999 49.9781 23.5708 -33.7319 80.8734bint= -1.3161 -56.5055 53.8732-38.5544 170.7431 -62.4187 -98.9014 -25.9364.3205 9.6342 -6.2326 -61.3311 48.866-10.4242 14.8869 12.6943 -45.363 70.7517 30.463 -22.6035 83.5294 34.4177 -15.1129 83.9484 33.5708 -21.5381 88.6797 6.4525 -51.6618 64.5667 -11.7111 -68.6822 45.26 -5.1286 -57.4971 47.2399 -22.2469 -68.6258 24.132結(jié)果stats= 0.924673.57410.0000751.6477回歸方程r2F值p值誤差方差估計(jì)rcoplot(r,rint)殘差分析xy=[ones(n,1),x,y,x.^2,x.*y,y.^2];[b,bint,r,rint,stats]=regress(z,xy)調(diào)用regress函數(shù)作回歸分析的Matlab程序2部分結(jié)果b= bint= -283.1555 -744.3545178.043617.7399 -11.037 46.516822.681 -93.3294 138.6914-0.1906 -0.7285 0.34720.0601 -4.1065 4.2267-0.4046 -9.6039 8.7946

stats= 0.954537.75250604.8832一室模型：將整個(gè)機(jī)體看作一個(gè)房室，稱中心室，室內(nèi)血藥濃度是均勻的?？焖凫o脈注射后，濃度立即上升；然后迅速下降。當(dāng)濃度太低時(shí)，達(dá)不到預(yù)期的治療效果；當(dāng)濃度太高，又可能導(dǎo)致藥物中毒或副作用太強(qiáng)。臨床上，每種藥物有一個(gè)最小有效濃度c1和一個(gè)最大有效濃度c2。設(shè)計(jì)給藥方案時(shí)，要使血藥濃度保持在c1~c2之間。本題設(shè)c1=10，c2=25(ug/ml).擬合問題實(shí)例2給藥方案——一種新藥用于臨床之前，必須設(shè)計(jì)給藥方案.藥物進(jìn)入機(jī)體后血液輸送到全身，在這個(gè)過程中不斷地被吸收、分布、代謝，最終排出體外，藥物在血液中的濃度，即單位體積血液中的藥物含量，稱為血藥濃度。案例

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

t(h)0.250.511.523468c(g/ml)19.2118.1515.3614.1012.899.327.455.243.01要設(shè)計(jì)給藥方案,必須知道給藥后血藥濃度隨時(shí)間變化的規(guī)律。從實(shí)驗(yàn)和理論兩方面著手：給藥方案1.在快速靜脈注射的給藥方式下，研究血藥濃度（單位體積血液中的藥物含量）的變化規(guī)律。tc2cc10

問題2.給定藥物的最小有效濃度和最大治療濃度，設(shè)計(jì)給藥方案：每次注射劑量多大；間隔時(shí)間多長分析理論：用一室模型研究血藥濃度變化規(guī)律

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

在此，d=300mg，t及c（t）在某些點(diǎn)處的值見前表，需經(jīng)擬合求出參數(shù)k、v用線性最小二乘擬合c(t)計(jì)算結(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))程序：給藥方案設(shè)計(jì)cc2c10

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

血藥濃度c(t)

應(yīng)c1

c(t)

c2初次劑量D0應(yīng)加大給藥方案記為：2、1、計(jì)算結(jié)果：給藥方案：c1=10,c2=25k=0.2347v=15.02故可制定給藥方案：即:

首次注射375mg，其余每次注射225mg，注射的間隔時(shí)間為4小時(shí)。估計(jì)水塔的流量2、解題思路3、算法設(shè)計(jì)與編程1、問題某居民區(qū)有一供居民用水的園柱形水塔，一般可以通過測(cè)量其水位來估計(jì)水的流量，但面臨的困難是，當(dāng)水塔水位下降到設(shè)定的最低水位時(shí)，水泵自動(dòng)啟動(dòng)向水塔供水，到設(shè)定的最高水位時(shí)停止供水，這段時(shí)間無法測(cè)量水塔的水位和水泵的供水量．通常水泵每天供水一兩次，每次約兩小時(shí).水塔是一個(gè)高12.2米，直徑17.4米的正園柱．按照設(shè)計(jì)，水塔水位降至約8.2米時(shí)，水泵自動(dòng)啟動(dòng)，水位升到約10.8米時(shí)水泵停止工作．表1是某一天的水位測(cè)量記錄，試估計(jì)任何時(shí)刻（包括水泵正供水時(shí)）從水塔流出的水流量，及一天的總用水量．流量估計(jì)的解題思路擬合水位~時(shí)間函數(shù)確定流量~時(shí)間函數(shù)估計(jì)一天總用水量擬合水位~時(shí)間函數(shù)

測(cè)量記錄看，一天有兩個(gè)供水時(shí)段（以下稱第1供水時(shí)段和第2供水時(shí)段），和3個(gè)水泵不工作時(shí)段（以下稱第1時(shí)段t=0到t=8.97，第2次時(shí)段t=10.95到t=20.84和第3時(shí)段t=23以后）．對(duì)第1、2時(shí)段的測(cè)量數(shù)據(jù)直接分別作多項(xiàng)式擬合，得到水位函數(shù)．為使擬合曲線比較光滑，多項(xiàng)式次數(shù)不要太高，一般在3~6．由于第3時(shí)段只有3個(gè)測(cè)量記錄，無法對(duì)這一時(shí)段的水位作出較好的擬合．

2、確定流量~時(shí)間函數(shù)

對(duì)于第1、2時(shí)段只需將水位函數(shù)求導(dǎo)數(shù)即可，對(duì)于兩個(gè)供水時(shí)段的流量，則用供水時(shí)段前后（水泵不工作時(shí)段）的流量擬合得到，并且將擬合得到的第2供水時(shí)段流量外推，將第3時(shí)段流量包含在第2供水時(shí)段內(nèi)．3、一天總用水量的估計(jì)

總用水量等于兩個(gè)水泵不工作時(shí)段和兩個(gè)供水時(shí)段用水量之和，它們都可以由流量對(duì)時(shí)間的積分得到。算法設(shè)計(jì)與編程1、擬合第1、2時(shí)段的水位，并導(dǎo)出流量2、擬合供水時(shí)段的流量3、估計(jì)一天總用水量4、流量及總用水量的檢驗(yàn)

1、擬合第1時(shí)段的水位，并導(dǎo)出流量設(shè)t，h為已輸入的時(shí)刻和水位測(cè)量記錄（水泵啟動(dòng)的4個(gè)時(shí)刻不輸入），第1時(shí)段各時(shí)刻的流量可如下得：1）c1=polyfit（t（1：10），h（1：10），3）；

%用3次多項(xiàng)式擬合第1時(shí)段水位，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時(shí)刻的流量

4）流量函數(shù)為：

2、擬合第2時(shí)段的水位，并導(dǎo)出流量設(shè)t，h為已輸入的時(shí)刻和水位測(cè)量記錄（水泵啟動(dòng)的4個(gè)時(shí)刻不輸入），第2時(shí)段各時(shí)刻的流量可如下得：1）c2=polyfit(t(10.9:21),h(10.9:21),3)；

%用3次多項(xiàng)式擬合第2時(shí)段水位，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時(shí)刻的流量4）流量函數(shù)為：

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

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

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

xx12=[xx1xx2]；

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

tp12=9：0.1：11；

x12=polyval（c12，tp12）；%x12輸出第1供水時(shí)段各時(shí)刻的流量擬合的流量函數(shù)為：在第2供水時(shí)段之前取t=20，20.8兩點(diǎn)的流水量，在該時(shí)刻之后（第3時(shí)段）僅有3個(gè)水位記錄，我們用差分得到流量，然后用這4個(gè)數(shù)值擬合第2供水時(shí)段的流量如下：

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

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

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

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

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

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

t3=20.8：0.1：24；

x3=polyval（c3，tp3）；%x3輸出第2供水時(shí)段（外推至t=24）各時(shí)刻的流量擬合的流量函數(shù)為：

3、一天總用水量的估計(jì)第1、2時(shí)段和第1、2供水時(shí)段流量的積分之和，就是一天總用水量．雖然諸時(shí)段的流量已表為多項(xiàng)式函數(shù)，積分可以解析地算出，這里仍用數(shù)值積分計(jì)算如下：

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

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

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

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

y=（y1+y2+y12+y3）*237.8*0.01；%一天總用水量（）計(jì)算結(jié)果：y1=146.2,y2=266.8,y12=47.4,y3=77.3，y=1250.4

4、流量及總用水量的檢驗(yàn)計(jì)算出的各時(shí)刻的流量可用水位記錄的數(shù)值微分來檢驗(yàn)．用水量y1可用第1時(shí)段水位測(cè)量記錄中下降高度968-822=146來檢驗(yàn)，類似地，y2用1082-822=260檢驗(yàn)．供水時(shí)段流量的一種檢驗(yàn)方法如下：供水時(shí)段的用水量加上水位上升值260是該時(shí)段泵入的水量，除以時(shí)段長度得到水泵的功率（單位時(shí)間泵入的水量），而兩個(gè)供水時(shí)段水泵的功率應(yīng)大致相等．第1、2時(shí)段水泵的功率可計(jì)算如下：

p1=(y12+260)/2；%第1供水時(shí)段水泵的功率（水量仍以高度計(jì)）

tp4=20.8：0.1：23；

xp2=polyval（c3，tp4）；%xp2輸出第2供水時(shí)段各時(shí)刻的流量

p2=(0.1*trapz(xp2)+260)/2.2；%第2供水時(shí)段水泵的功率（水量仍以高度計(jì)）計(jì)算結(jié)果：p1=154.5，p2=140.1計(jì)算結(jié)果流量函數(shù)為：流量曲線見圖n=(3,4)n=(5,6)練習(xí)1用給定的多項(xiàng)式，如y=x3-6x2+5x-3，產(chǎn)生一組數(shù)據(jù)(xi,yi，i=1,2,…,n),再在yi上添加隨機(jī)干擾(可用rand產(chǎn)生(0,1)均勻分布隨機(jī)數(shù),或用rands產(chǎn)生N(0,1)分布隨機(jī)數(shù))，然后用xi和添加了隨機(jī)干擾的yi作的3次多項(xiàng)式擬合，與原系數(shù)比較。如果作2或4次多項(xiàng)式擬合，結(jié)果如何？

練習(xí)2、用電壓V=10伏的電池給電容器充電，電容器上t時(shí)刻的電壓為，其中V0是電容器的初始電壓，是充電常數(shù)。試由下面一組t，V數(shù)據(jù)確定V0，。作業(yè)：農(nóng)作物施肥效果分析電力市場(chǎng)的輸電阻塞管理（題目）電力市場(chǎng)的輸電阻塞管理（論文一）電力市場(chǎng)的輸電阻塞管理（論文二）用非線性最小二乘擬合c(t)-用lsqcurvefit2、主程序lihe2.m如下cleartdata=[0.250.511.523468];cdata=[19.2118.1515.3614.1012.899.327.455.243.01];x0=[10,0.5];x=lsqcurvefit('curvefun3',x0,tdata,cdata);f=curvefun3(x,tdata)x1、用M-文件curvefun3.m定義函數(shù)functionf=curvefun3(x,tdata)d=300f=(x(1)\d)*exp(-x(2)*tdata)%x(1)=v;x(2)=k

1.lsqcurvefit已知數(shù)據(jù)點(diǎn)：

xdata=（xdata1，xdata2，…，xdatan），

ydata=（ydata1，ydata2，…，ydatan）

用MATLAB作非線性最小二乘擬合

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

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

輸入格式為:(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是一個(gè)事先建立的定義函數(shù)F(x,xdata)

的M-文件,自變量為x和xdata說明：x=lsqcurvefit(‘fun’,x0,xdata,ydata,options);迭代初值已知數(shù)據(jù)點(diǎn)選項(xiàng)見無約束優(yōu)化

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）輸入格式為：

1）x=lsq