小木虫 --- 500万硕博科研人员喜爱的学术科研平台

首页 >> 计算模拟 >>查看话题

如何使用Matlab编程进行参数拟合


1前言
2基本概念和原理
3主要内容
4实例
5涉及的文件


1前言

之前帮疯学网做过一个利用Matlab编程进行参数拟合 的教程,由于疯学网好像倒闭了,希望之前做的工作不要白费,这里拿出来分享下,希望能对虫友的学习、科研工作有所帮助。其他的不多说,言归正传,下面从原理和实例对如何使用Matlab编程进行参数拟合进行讲解。

2基本概念和原理

所谓参数拟合,就是已知试验或者真实数据,然后寻找一个模型对其规律进行模拟,求取模型中未知参数的一个过程。根据模型的复杂程度我分成以下几类讲解:
代数方程模型
线性模型
非线性模型
单变量,单目标问题
多变量,单目标问题
多变量,多目标问题,共享参数
复数模型拟合
微分方程模型
简单微分方程参数拟合
复杂微分方程参数拟合

3主要内容


3.1代数方程模型

y=f(x,β)
x, n维自变量x='
β,p维参数向量β ='
y,m维因变量y='
f,m维函数向量f='
Matlab 求解函数
线性最小二乘法:ployfit, regress
非线性最小二乘法:fit, lsqnonlin,lsqcurvefit,nlinfit
nlintool,fmincon

3.2微分方程模型

dx/dt=f(t,x,β,u)  x(t0)=x0
x, n维状态变量x='
x0, n维状态变量初值x='
β,p维参数向量β ='
u,r维操作变量y='
f,ODE方程m维函数向量f='
对于给定β,由龙格-库塔积分可得x(t)
y(t)=g(x(t), β)
y为m维输出量y='
g为输出量y与状态变量x之间非线性关系的函数向量g='
用导数法和积分法将ODE问题转化为代数方程问题

3.2优化准则和参数初值确定方法

优化准则:最小二乘法,极大似然,概率密度函数
非线性模型必须采用非线性回归的方法
参数初值确定方法
(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决定性指标

{R^2}
回归均值的总偏离
Ne:实验点数目,
yei,yci,i次实验值与计算值
{R^2} = 1 - {{\sum\limits_{i = 1}^{{N_e}} {{{\left( {{y_{ei}} - {y_{ci}}} \right)}^2}} } \over {\sum\limits_{i = 1}^{{N_e}} {y_{ei}^2} }}
由于公式比较难显示,见附件ppt里面内容


4实例


4.1线性拟合函数:regress()

调用格式:b=regress(y,X)
= regress(y,X)
= 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=
y=x*+normrnd(0,0.1,10,1)
=regress(y,x,0.05)
rcoplot(r,rint)

4.2简单线性模型-多项式拟合

多项式曲线拟合函数:polyfit( )
调用格式:p=polyfit(x,y,n)
                = polyfit(x,y,n)
x,y为数据点,n为多项式阶数,返回p为幂次从高到低的多项式系数向量p。矩阵s用于生成预测值的误差估计。
多项式曲线求值函数:polyval( )
调用格式:y=polyval(p,x)
                =polyval(p,x,s)
y=polyval(p,x)为返回对应自变量x在给定系数P的多项式的值
多项式曲线拟合的评价和置信区间函数:polyconf( )
调用格式:=polyconf(p,x,s)
                =polyconf(p,x,s,alpha)
多项式输出 ps = poly2str(p,'x')   
codefile Fit_polynomial  
                     

4.3稳健回归函数:robustfit( )

稳健回归是指此回归方法相对于其他回归方法而言,受异常值的影响较小。
调用格式:       
b=robustfit(x,y)
=robustfit(x,y)
=robustfit(x,y,’wfun’,tune,’const’)
例:演示一个异常数据点如何影响最小二乘拟合值与稳健拟合。首先利用函数y=10-2x加上一些随机干扰的项生成数据集,然后改变一个y的值形成异常值。调用不同的拟合函数,通过图形观查影响程度。
code File Robust_Fit

4.4 非线性问题使用matlab函数-lsqcurvefit,lsqnonlin

= 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和 options
min 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://emuch.net/bbs/viewthread.php?tid=3866180
code file  percolation_fit

4.8 复数模型拟合

将模型分离成实部和虚部
求解如下优化模型
\
例子
http://emuch.net/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积分可得对应时间点的浓度预测值Cp
2.以浓度的预测值Cp与实验值Ce的残差平方和为优化目标,利用lsqnonlin或者fmincon进行搜索获得最优的ODE参数,每搜一次,调用ode45计算预测浓度值,直到优化完成
\min \sum\limits_{{\text{i}} = 1}^n {{{\left( {{C_{i,experiment}} - {{\text{C}}_{{\text{i,ODE calculated}}}}} \right)}^2}}
http://emuch.net/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)=',根据下表数据用最优化方法估计参数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
…….
微分方程:
{{{d{x_1}} \over {dt}} = {k_5} - q{x_1} - {k_1}{x_1}{x_2} - {k_4}{x_1}{x_6}\sqrt {0.9} }
{{{d{x_2}} \over {dt}} = 7.0 - qx_2 - {k_1}{x_1}{x_2} - 2{k_2}{x_2}{x_3}}
{{{d{x_3}} \over {dt}} = 1.75 - q{x_3} - {k_2}{x_2}{x_3}}
{{{d{x_4}} \over {dt}} = - q{x_4} + 2{k_1}{x_1}{x_2} - {k_3}{x_4}{x_5}}
{{{d{x_5}} \over {dt}} = - q{x_5} + 3{k_2}{x_2}{x_3} - {k_3}{x_4}{x_5}}
{{{d{x_6}} \over {dt}} = - q{x_6} + 2{k_3}{x_4}{x_5} - {k_4}{x_1}{x_6}\sqrt {0.9} }
{{{d{x_7}} \over {dt}} = - q{x_7} + 2{k_4}{x_1}{x_6}\sqrt {0.9} }
{q = 8.75 + {k_5}}
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=;
>> R=;
>> x = ';
>> MSR=;
>> y = MSR';
beta=nlinfit(x,y,myfun,
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_fit
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=';
R=';
xdata = ;
MSR=;
ydata = MSR';
%beta0=
x0=;
% =nlinfit(xdata,ydata,myfun,x0);
% ci = nlparci(x,residual,'covar',covb)
%% =======利用lsqcurvefit 非线性最小二乘法
=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=;
title(Text)
%% ==========利用globalsearch和 fmincon=========
tic
x0=;
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',,'ub',,'options',opts);
= run(gs,problem)
t1=toc
figure;plot(1:numel(ydata),ydata,'r-*',1:numel(ydata),myfun(xming,xdata),'bo-')
legend('real','model')
Text1=;
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',,'ub',,'options',opts);
=run(ms,problem,20)
t2=toc
figure;plot(1:numel(ydata),ydata,'r-*',1:numel(ydata),myfun(xminm,xdata),'bo-')
legend('real','model')
Text2=;
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);
= ga(Fgen,4,,,,,, )
=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=;
title(Text3)
结果如下




由图可全局算法估计参数结果优于一次非线性优化估计的参数

5涉及的文件

(见附件)
本贴录有视频,如需要私信。

用户评论

顶一下,感谢分享!

感谢分享,支持一下。

顶一下,感谢分享!

顶一下,感谢分享!

顶起,为版主加油!

顶一下,支持!!!

我也在用最小二乘法求二次多项式系数能加一下qq细聊吗?953058630:D

确实不错!支持一下!

顶一下,感谢分享!

楼主很专业!感谢分享!!!

顶一下,感谢分享!


研究生必备
与500万研究生在线互动!


扫描下载送金币