如何使用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简单微分方程参数拟合

1. 由给定的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）

#往期热文推荐#