如何使用Matlab编程进行参数拟合
1前言 2基本概念和原理 3主要内容 4实例 5涉及的文件 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，nlinfit nlintool，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.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决定性指标 回归均值的总偏离 Ne：实验点数目， yei，yci，i次实验值与计算值 由于公式比较难显示，见附件ppt里面内容  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β+ε β是p1的参数向量；ε是服从标准正态分布的随机干扰的n1的向量；y为n1的因变量向量；X为np自变量矩阵。 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和 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://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积分可得对应时间点的浓度预测值Cp 2.以浓度的预测值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_fit 结果如下     由图可全局算法估计参数结果优于一次非线性优化估计的参数 5涉及的文件 （见附件）本贴录有视频，如需要私信。 |