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β+ε
β是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和 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涉及的文件 

（见附件）本贴录有视频，如需要私信。