如何使用Matlab編程進行參數擬合

1 前言

之前幫瘋學網做過一個利用Matlab編程進行參數擬合 的教程，由於瘋學網好像倒閉了，希望之前做的工作不要白費，這裡拿出來分享下，希望能對蟲友的學習、科研工作有所幫助。其他的不多說，言歸正傳，下面從原理和實例對如何使用Matlab編程進行參數擬合進行講解。

2 基本概念和原理

所謂參數擬合，就是已知試驗或者真實數據，然後尋找一個模型對其規律進行模擬，求取模型中未知參數的一個過程。根據模型的複雜程度我分成以下幾類講解：

代數方程模型

線性模型

非線性模型

單變量，單目標問題

多變量，單目標問題

多變量，多目標問題，共享參數

複數模型擬合

微分方程模型

簡單微分方程參數擬合

複雜微分方程參數擬合

3 主要內容

3.1代數方程模型

y=f(x,β)

x, n維自變量x=[x1,x2,…,xn]'

β，p維參數向量β =[β1, β2,…, βn]'

y，m維因變量y=[y1,y2,…,yn]'

f，m維函數向量f=[f1,f2,…,fn]'

Matlab 求解函數

線性最小二乘法：ployfit， regress

非線性最小二乘法：fit， lsqnonlin，lsqcurvefit，nlinfitnlintool，fmincon

3.2微分方程模型

dx/dt=f(t,x,β,u)  x(t0)=x0

x, n維狀態變量x=[x1,x2,…,xn]'

x0, n維狀態變量初值x=[x10,x20,…,xn0]'

β，p維參數向量β =[β1, β2,…, βp]'

u，r維操作變量y=[u1,u2,…,um]'

f，ODE方程m維函數向量f=[f1,f2,…,fm]'

對於給定β，由龍格-庫塔積分可得x(t)

y(t)=g(x(t), β)

y為m維輸出量y=[y1,y2,…,yn]'

g為輸出量y與狀態變量x之間非線性關係的函數向量g=[g1,g2,…,gn]'

用導數法和積分法將ODE問題轉化為代數方程問題

3.3優化準則和參數初值確定方法

優化準則：最小二乘法，極大似然，概率密度函數

非線性模型必須採用非線性回歸的方法

參數初值確定方法

（1）根據模型結構 和本質

描述物理系統的模型的結構和本質可能指明未知參數的取值範圍（2）模型方程的漸進行為

例如指數衰減模型yi=k1+k2exp(-k3*xi)

xi-->∞，yi約為k1， xi-->0，yi=k1+k2(1-k3*xi)，利用接近0的x及對應y回歸，結合k1,就可得到k1 k2 k3初值（3）模型方程的變換

把非線性系統轉化為線性系統，線性最小二乘結果作為非線性估計的初值（4）直接搜索，全局算法

如果對參數不了解，可用直接搜索或者全局（多起點，遺傳等算法處理）

3.4決定性指標

回歸均值的總偏離

Ne：實驗點數目，

yei，yci，i次實驗值與計算值

由於公式比較難顯示，見附件ppt裡面內容

如何使用Matlab編程進行參數擬合

4 實例

4.1線性擬合函數：regress()

調用格式：b=regress(y,X)

[b,bint,r,rint,stats]= regress(y,X)

[b,bint,r,rint,stats]= regress(y,X,alpha)

該函數求解線性模型：y=Xβ+ε

β是p?1的參數向量；ε是服從標準正態分布的隨機幹擾的n?1的向量；y為n?1的因變量向量；X為n?p自變量矩陣。

bint返回β的95%的置信區間。r中為形狀殘差，rint中返回每一個殘差的95%置信區間。Stats向量包含R2統計量、回歸的F值和p值。

例 設y的值為給定的x的線性函數加服從標準正態分布的隨機幹擾值得到。即y=10+x+ε ；求線性擬合方程係數。

x=[ones(10,1) (1:10)']

y=x*[10;1]+normrnd(0,0.1,10,1)

[b,bint, r,rint,stats]=regress(y,x,0.05)

rcoplot(r,rint)

4.2簡單線性模型-多項式擬合

多項式曲線擬合函數：polyfit( )

調用格式：p=polyfit(x,y,n)

[p,s]= polyfit(x,y,n)

x,y為數據點，n為多項式階數，返回p為冪次從高到低的多項式係數向量p。矩陣s用於生成預測值的誤差估計。

多項式曲線求值函數：polyval( )

調用格式：y=polyval(p,x)

[y,DELTA]=polyval(p,x,s)

y=polyval(p,x)為返回對應自變量x在給定係數P的多項式的值多項式曲線擬合的評價和置信區間函數：polyconf( )調用格式：[Y,DELTA]=polyconf(p,x,s)

[Y,DELTA]=polyconf(p,x,s,alpha)

多項式輸出 ps = poly2str(p,'x')

codefile Fit_polynomial

4.3穩健回歸函數：robustfit( )

穩健回歸是指此回歸方法相對於其他回歸方法而言，受異常值的影響較小。

調用格式：

b=robustfit(x,y)

[b,stats]=robustfit(x,y)

[b,stats]=robustfit(x,y,』wfun』,tune,』const』)例：演示一個異常數據點如何影響最小二乘擬合值與穩健擬合。首先利用函數y=10-2x加上一些隨機幹擾的項生成數據集，然後改變一個y的值形成異常值。調用不同的擬合函數，通過圖形觀查影響程度。

code File Robust_Fit

4.4 非線性問題使用matlab函數-lsqcurvefit，lsqnonlin

[x resnorm] = lsqcurvefit(fun,x0,xdata,ydata,...)fun   是我們需要擬合的函數，

x0    是我們對函數中各參數的猜想值，

xdata  則是自變量的值

ydata  是因變量的值，目標值

min  sum {(FUN(X,XDATA)-YDATA).^2}

x = lsqnonlin(fun,x0)

x0為初始解向量；fun為優化函數，fun返回向量值F，而不是平方和值，平方和隱含在算法中，fun的定義與前面相。

x = lsqnonlin(fun,x0,lb,ub)

lb、ub定義x的下界和上界

x = lsqnonlin(fun,x0,lb,ub,options)

options為指定優化參數，若x沒有界，則lb=[ ]，ub=[ ]

min  sum {FUN(X).^2}

4.5非線性問題使用matlab函數-fmincon

x= fmincon(fun,x0,A,b)

給定初值x0，求解fun函數的最小值x。fun函數的約束條件為A*x<= b，x0可以是標量或向量。

x= fmincon(fun,x0,A,b,Aeq,beq)

最小化fun函數，約束條件為Aeq*x= beq 和 A*x <= b。若沒有不等式線性約束存在，則設置A=[]、b=[]。

x= fmincon(fun,x0,A,b,Aeq,beq,lb,ub)

定義設計變量x的線性不等式約束下界lb和上界ub，使得總是有lb<= x <= ub。若無等式線性約束存在，則令Aeq=[]、beq=[]。

x= fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon)

在上面的基礎上，在nonlcon參數中提供非線性不等式c(x)或等式ceq(x)。

fmincon函數要求c(x) <= 0且ceq(x)= 0。

x= fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)用options參數指定的參數進行最小化。

x= fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options,P1,P2,...) 將問題參數P1, P2等直接傳遞給函數fun和nonlin。若不需要這些變量，則傳遞空矩陣到A, b, Aeq, beq, lb, ub, nonlcon和 optionsmin F(X) fun函數返回一向量或者標量

4.6非線性問題線性化

4.7單變量，單目標問題（cftool 的應用）

例子

滲流公式為y=A*(x-xc)^p

其中x為自變量，y為因變量，A、xc和p均為常數

為了測試模擬，設定A=18.5，xc=0.095，p=-2.3，得到以下數據x                y

---

-0.1001        3.5E+06

0.1002        3.3E+06

0.11           2.9E+05

0.12           9.0E+04

0.15           1.5E+04

0.2             3.3E+03

0.3             7.1E+02

0.4             2.8E+02

0.5             1.5E+02

0.6             8.9E+01

http://muchong.com/bbs/viewthread.php?tid=3866180

code file  percolation_fit

4.8 複數模型擬合

將模型分離成實部和虛部

求解如下優化模型

例子

http://muchong.com/bbs/viewthread.php?tid=4268659&target=self&page=1

模型

變量參數為：a,b,c,d1,m1,n1,d2,m2,n2,d3,m3,n3,擬合曲線為 e=e1+ie2=a-b/(x^2+icx)+m1/(n1-x^2-id1x)+m2/(n2-x^2-id2x)+m3/(n3-x^2-id3x)

Data

x                e1             e2

5.16957        1.854        4.5139

4.96279        1.9351        4.5777

4.77192        2.0221        4.4781

…

code file  complexfit

4.9簡單微分方程參數擬合

由給定的ODE參數初值結合ODE方程dC/dt=f（k，c，t），利用ode45積分可得對應時間點的濃度預測值Cp2.以濃度的預測值Cp與實驗值Ce的殘差平方和為優化目標，利用lsqnonlin或者fmincon進行搜索獲得最優的ODE參數，每搜一次，調用ode45計算預測濃度值，直到優化完成

http://muchong.com/bbs/viewthread.php?tid=4685934&target=self&page=1

問題

實驗數據為：（t,c）=（0,0.69）（2,0.645）（4,0.635）（8,0.62）（24,0.61）（48,0.61）.

其中t為時間，c為某離子的濃度。

動力學方程模型為：-dc/dt=k*(c0-c)^(1/3)*(c-c~).

其中c0為初始濃度可以取0.7，c~為平衡濃度取0.61.

code file ODE_parafit

4.10複雜微分方程參數擬合

在一個等溫間歇反應器中進行複雜反應，描述該反應系統的七個微分方程有5個參數k1, k2, k3, k4, k5, 初始狀態x(0)=[0.1883 0.2507 0.0467 0.0899 0.1804 0.1394 0.1046]',根據下表數據用最優化方法估計參數k

t            y1(x1)    y2(x4)    y3(x5)    y4(x6)

0        0.1883    0.0899    0.1804    0.1394

0.0100    0.2047    0.0866    0.1729    0.1297

0.0200    0.2181    0.0856    0.1680    0.1205…….

微分方程：

code file: KineticsEst5.m

4.11 多元非線性擬合，全局優化

例子

https://zhidao.baidu.com/question/168077392.html

matlab nlinfit多元非線性擬合 錯在哪裡？

需要擬合如下形式的模型：

y = beta(1)+beta(2)*V+exp(beta(3)+beta(4)*R)

我使用的代碼如下：

>> V=[47,69,54,65.6,46,61,66.5,50.8,55.9,44.5,53.3,54,56.8,51.9,60.1,67.3,62.4,61.1,61.7,76,51.7,50.7,57.1,65.6,53.1,60.6,59.1,46,49,47.3,75.3,52.5,58.3,54.2,46.8,55.1,47.8];

>> R=[47,1081,186,127,50,109,397,178,189,72,86,98,273,186,202,137,138,104,672,270,94,57,59,187,129,55,288,77,155,65,972,85,884,195,105,172,97];

>> x = [V;R]';

>> MSR=[18,11,11,16,20,23,26,2.9,11.1,12.2,22.1,19.8,8.6,4.3,17.6,35.9,27.2,30.7,11.9,42.3,16.7,31.9,41.6,28.2,12.2,50.9,12.2,12.6,1.9,20.7,34.6,21,5.1,7.7,5.4,10.9,9.2];

>> y = MSR';

beta=nlinfit(x,y,myfun,[-100,1,-10,1000]

m-函數為：

function y  = myfun(beta,x)

y = beta(1)+beta(2)*x(:,1)+exp(beta(3)+beta(4)*x(:,2));錯誤的原因是：Input argument "beta" is undefined.beta沒有定義？

請問錯在哪裡？在線等待哦，謝謝！

在命令窗口輸入yy=myfun(beta,x),回車運行試試，也是不行的，

??? Error using ==> beta at 21

Not enough input arguments.

code file:Muti_var_fitCODE:

function Muti_var_fit

% matlab nlinfit多元非線性擬合錯在哪裡？

% 舉報|2010-07-19 11:59transtorseu | 分類：其他程式語言 | 瀏覽1500次% 需要擬合如下形式的模型：

%y = beta(1)+beta(2)*V+exp(beta(3)+beta(4)*R)%http://zhidao.baidu.com/link?url=7Jue1nv1dhSE2OpGMv6pKWOW7NX0zgmrpAZV6usGpavPA8PSUJSWY0Hp0AdgTkdbmdBTnYnLaIJXg4Z8wt2r8_myfun=@(x,xdata)x(1)+x(2)*xdata(:,1)+exp(x(3)+x(4)*xdata(:,2));V=[47,69,54,65.6,46,61,66.5,50.8,55.9,44.5,53.3,54,56.8,51.9,60.1,67.3,62.4,61.1,61.7,76,51.7,50.7,57.1,65.6,53.1,60.6,59.1,46,49,47.3,75.3,52.5,58.3,54.2,46.8,55.1,47.8]';R=[47,1081,186,127,50,109,397,178,189,72,86,98,273,186,202,137,138,104,672,270,94,57,59,187,129,55,288,77,155,65,972,85,884,195,105,172,97]';xdata = [V R];

MSR=[18,11,11,16,20,23,26,2.9,11.1,12.2,22.1,19.8,8.6,4.3,17.6,35.9,27.2,30.7,11.9,42.3,16.7,31.9,41.6,28.2,12.2,50.9,12.2,12.6,1.9,20.7,34.6,21,5.1,7.7,5.4,10.9,9.2];ydata = MSR';

%beta0=[-80 1 4 -0.01]

x0=[-20 1 1 0.01];

% [x,residual,jacobian,covb,resnorm]=nlinfit(xdata,ydata,myfun,x0);% ci = nlparci(x,residual,'covar',covb)

%% =======利用lsqcurvefit 非線性最小二乘法

[x,resnorm,residual,exitflag,output,lambda,jacobian]=lsqcurvefit(myfun,x0,xdata,ydata);ci = nlparci(x,residual,jacobian)

fprintf('\n\n使用函數lsqcurvefit()估計得到的參數值為:\n')fprintf('\t待求參數 k1 = %.6f\n',x(1))

fprintf('\t待求參數 k2 = %.6f\n',x(2))

fprintf('\t待求參數 k3 = %.6f\n',x(3))

fprintf('\t待求參數 k4 = %.6f\n',x(4))

fprintf('  The sum of the squares is: %.3e\n\n',resnorm)figure;plot(1:numel(ydata),ydata,'r-*',1:numel(ydata),myfun(x,xdata),'bo-')legend('real','model')

Text=['使用局部優化函數lsqcurvefit估計得到的參數'];title(Text)

%% ==========利用globalsearch和 fmincon=========tic

x0=[-10 1 1 0.1];

F =@(x)norm(x(1)+x(2)*xdata(:,1)+exp(x(3)+x(4)*xdata(:,2))-ydata);gs = GlobalSearch('Display','iter');

opts=optimset('Algorithm','trust-region-reflective');problem=createOptimProblem('fmincon','x0',x0,...

'objective',F,'lb',[-100,-100,-100,-1],'ub',[100,100,100,1],'options',opts);[xming,fming,flagg,outptg,manyminsg] = run(gs,problem)t1=toc

figure;plot(1:numel(ydata),ydata,'r-*',1:numel(ydata),myfun(xming,xdata),'bo-')legend('real','model')

Text1=['利用全局算法globalsearch和 fmincon估計得到的參數,耗時',num2str(t1),'s'];title(Text1)

%% =========全局算法 multistart和lsqcurvefit

tic

ms=MultiStart;

opts=optimset('Algorithm', 'trust-region-reflective');problem=createOptimProblem('lsqcurvefit','x0',x0,'xdata',...

xdata, 'ydata',ydata,'objective',myfun,'lb',[-100,-100,-100,-1],'ub',[100,100,100,1],'options',opts);[xminm,fminm,flagm,outptm,manyminsm]=run(ms,problem,20)t2=toc

figure;plot(1:numel(ydata),ydata,'r-*',1:numel(ydata),myfun(xminm,xdata),'bo-')legend('real','model')

Text2=['利用全局算法 multistart和lsqcurvefit估計得到的參數,耗時',num2str(t2),'s'];title(Text2)

%% =============遺傳算法的參數識別=======

tic

Fgen =@(x)norm(x(1)+x(2)*xdata(:,1)+exp(x(3)+x(4)*xdata(:,2))-ydata);options = gaoptimset('TolFun',1e-10);

[xgen,fval] = ga(Fgen,4,[],[],[],[],[-100;-100;-100;-1], [100;100;100;1])[xgen,resnorm,residual,exitflag,output,lambda,jacobian]=lsqcurvefit(myfun,xgen,xdata,ydata);ci = nlparci(xgen,residual,jacobian)

t3=toc

figure;plot(1:numel(ydata),ydata,'r-*',1:numel(ydata),myfun(xgen,xdata),'bo-')legend('real','model')

Text3=['利用全局算法遺傳算法的參數識別估計得到的參數,耗時',num2str(t3),'s'];title(Text3)

結果如下

由圖可全局算法估計參數結果優於一次非線性優化估計的參數

（文章來源：小木蟲  作者：dbb627）

#往期熱文推薦#

除了Sci-Hub，免費文獻下載還有誰？

如何在論文中畫出漂亮的插圖？

科研工作效率慢？這十款免費科研利器來幫你！

iSlide：你與 PPT 達人只有這一個插件的距離

常用電化學儀器及生產廠家