曲线拟合实例

例如：已知数据队列buf=【5410。

】x取值1：n n是队列长度函数f（x）=a+b*sin(c*x+d) .avg 是队列平均值a b c d 为参数a范围（2/3,1）*avgb范围（0,1/3）*avgc的范围（0,24*pi）d （0,2*pi）1、首先定义目标函数function y=ga_curfit(x)global ydata nt=1:n;y=0;for i=1:ny=y+(ydata(i)-(x(:,1)+x(:,2).*sin(x(:,3).*t(i)+x(:,4)))).^2/n; endy=sqrt(y);end2、把数据b.txt放在工作空间目录中然后再命令窗口中输入clearglobal ydata nformat long gload b.txtydata=b';n=length(ydata);avg=sum(ydata)/n;LB=[2/3*avg 0 0 0];UB=[1*avg 1/3*avg 24*pi 2*pi];nvars=4;options=gaoptimset;options=gaoptimset(options,'PopulationSize',300); options=gaoptimset(options,'CrossoverFraction',0.8); options=gaoptimset(options,'MigrationFraction',0.1);options=gaoptimset(options,'Generations',500);options = gaoptimset(options,'TolFun', 1e-50);%options = gaoptimset(options,'InitialPopulation',final_pop);options = gaoptimset(options,'Display', 'final');options = gaoptimset(options,'PopInitRange', [LB;UB]);options = gaoptimset(options,'PlotFcns', @gaplotbestf);options=gaoptimset(options,'Vectorize','on');%目标函数向量化[x,fval,exitflag,output,final_pop,scores]=ga(@ga_curfit,nvars,[],[],[],[],LB,UB,[],options);t=1:n;plot(t,ydata,'r*');hold onplot(t,x(1)+x(2)*sin(x(3)*t+x(4)))legend('数据','拟合')引言之前曾发帖讨论过常微分方程参数拟合问题，就实际应用而言，还是以非线性代数方程参数拟合问题居多。

实现非线性代数方程参数拟合的软件有很多，比如MATLAB、Origin、SPSS和1stOpt等，特别是1stOpt，其代码简单易学，几乎不需调试，即可获得高质量的拟合结果，不过需要指出的是，该软件的拟合过程对操作者而言是一个黑箱。

拟合问题，实际上是一类最优化问题。

既然说到最优化问题，自然要涉及最优化算法、初值、程序代码。

初值的选取一直是困扰局部最优算法的问题，比如在调用MATLAB著名的lsqnonlin函数时，拟合结果对选用的初值具有很大的依赖性。

以具有全局搜索能力算法的计算结果，作为初值供局部最优化算法调用，是一种解决初值选取问题的方法。

本帖的主要目的正是讨论一种基于遗传算法拟合非线性代数方程参数的方法。

1. 本帖讨论的例子出处本帖讨论的问题与数据，取自2013年全国研究生数学建模竞赛E题，直观起见，现将问题表述如下：已知自变量p，因变量L，L和p满足以下函数式：L=(1-(1-p)^k3)^k1*(1-(1-p)^k4)^k2；且k1，k2≥0，k1+k2≥1 )；（需要指出的是，该公式中参数可进一步合并以防止过拟合，拟合结果不唯一，不过，在本帖中并未作合并处理）拟合出k1，k2，k3和k4。

2.结果与讨论运行MATLAB软件，调用lsqnonlin函数，参数初值按k1=1，k2=0，k3=1，k4=0作试探（先不管k1，k2的约束条件，实际上，k1，k2的约束常可自然满足），程序代码如下：clear all;clcformat longdata=;p=data(:,1); %pi Lexp=data(:,2); %Li k0=;lb=;ub=*1e6;%-------------------------------------------------------------------------% 使用函数lsqnonlin()进行参数估计OPTIONS=optimset('MaxFunEvals',1000,'TolFun',1e-12,'Algorithm','trust-region-reflective' );= ...lsqnonlin(@ObjFunc,k0,lb,ub,OPTIONS,p,Lexp); ci = nlparci(k,residual,jacobian); %residual;y=Objfit(p,k);fprintf('\n\n拟合结果:\n') fprintf(' \t残差平方和= %.6e\n\n',resnorm); fprintf('\n\t参数 a = %.9f',k(1)) fprintf('\n\t参数 b = %.9f',k(2)) fprintf('\n\t参数 c = %.9f',k(3)) fprintf('\n\t参数 d = %.9f',k(4)) n=length(p);R2=1-sum((Lexp-y).^2)./sum((Lexp-mean(Lexp)).^2);MSE=1/n*sum((Lexp-y).^2);MAE=1/n*sum(abs(Lexp-y));MAS=max(abs(Lexp-y));fprintf('\n\t决定系数R-Square = %.9f',R2); fprintf('\n\t均方误差MSE = %.9f',MSE); fprintf('\n\t平均绝对误差MAE = %.6f',MAE); fprintf('\n\t最大绝对误差MAS = %.6f',MAS); figure(2)plot(p,y,'b',p,Lexp,'or'),axis(),...text(0.05,0.95,),...text(0.05,0.85,),...text(0.05,0.75,),...text(0.05,0.65,),...xlabel('p'),ylabel('L'), legend('拟合结果','原始数据','Location','Best') %-------------------------------------------------------------------------function f = ObjFunc(k,p,Lexp) f=Objfit(p,k)-Lexp;%------------------------------------------------------------------------fun ction f = Objfit(p,k)f=(1-(1-p).^k(3)).^k(1).*(1-(1-p).^k(4)).^k(2);计算结果：残差平方和= 7.150814e-001参数 a = 1.000000003参数 b = 0.000000188参数 c = 1.000000000参数 d = 0.000000000决定系数R-Square = 0.622496369均方误差MSE = 0.042063615平均绝对误差MAE = 0.186062最大绝对误差MAS = 0.288096 图形结果见附图1。

由结果可见，决定系数R-Square = 0.622，拟合结果相当不理想，于是以其他初值作试探，结果如下：（k1 k2 k3 k4）决定系数R-Square （0 0 0 0） -3.722589381 （ 1 0 0 0） 0.999998738 （0 1 0 0） 0.999836583 （0 0 1 0） -3.723048781 （0 0 0 1） -3.723048781 （ 1 1 0 0） -1.136614242 （ 1 0 1 0） 0.622496369 …… ……可见拟合结果对初值的依赖性明显。

对于调用MATLAB自带的遗传算法GA函数，以其计算结果作为初值，代码如下：clear all;close all;clc data=;p=data(:,1); %pi Lexp=data(:,2); %Li %-------------遗传算法----------------------------------------------------- options = gaoptimset('Generations',1000,'StallGenLimit',300,...'StallTimeLimit',50,'TolFun',1e-12,'TolCon',1e-12); =ga(@objfun,4,options);options = gaoptimset('InitialPopulation',final_pop,'Generations',1000,'StallGenLimit',300,...'StallTimeLimit',50,'TolFun',1e-12,'TolCon',1e-12); =ga(@obj fun,4,options);fprintf('\n\n遗传算法的初始估计数值:\n'); fprintf('\n\t参数a0 = %.9f',k2(1)); fprintf('\n\t参数b0 = %.9f',k2(2)); fprintf('\n\t参数c0 = %.9f',k2(3)); fprintf('\n\t参数d0 = %.9f',k2(4)); %------------------------------最小二乘法---------------------------------- k0=k2;lb=;ub=*1e6;OPTIONS=optimset('MaxFunEvals',1000,'TolFun',1e-12,'Algorithm','trust-region-reflective', 'Display','Off');= ...lsqnonlin(@ObjFunc,k0,lb,ub,OPTIONS,p,Lexp); %-----------------------------------------------------------------%------------------结果输出与图形化---------------------------------------------- k=k3;y=Lorentz(p,k);fprintf('\n\n最终拟合结果:\n') fprintf(' \t残差平方和= %.6e\n\n',resnorm); fprintf('\n\t参数 a = %.9f',k(1)) fprintf('\n\t参数 b = %.9f',k(2)) fprintf('\n\t参数 c = %.9f',k(3)) fprintf('\n\t参数 d = %.9f',k(4)) n=length(p);R2=1-sum((Lexp-y).^2)./sum((Lexp-mean(Lexp)).^2);MSE=1/n*sum((Lexp-y).^2);MAE=1/n*sum(abs(Lexp-y));MAS=max(abs(Lexp-y));fprintf('\n\t决定系数R-Square = %.9f',R2); fprintf('\n\t均方误差MSE = %.6f',MSE); fprintf('\n\t平均绝对误差MAE = %.6f',MAE); fprintf('\n\t最大绝对误差MAS = %.6f',MAS); figureplot(p,y,'b',p,Lexp,'or'),axis(),...text(0.05,0.95,),...text(0.05,0.85,),...text(0.05,0.75,),...text(0.05,0.65,),...xlabel('p'),ylabel('L'), legend('拟合结果','原始数据','Location','Best') %-------------------------------------------------------------------------function f = ObjFunc(k,p,Lexp) f=Lorentz(p,k)-Lexp;%----------------构造拟合目标函数------------------------------------------- function fun=objfun(k) data=;p=data(:,1); %pi Lexp=data(:,2); %Li n=length(p);g1=k(1)+k(2);if (k(1)<0||k(2)<0||g1<1)fun=inf;elsefor i=1:nFF(i)=(Lexp(i)-Lorentz(p(i),k))^2; endfun=sum(FF);end%------------------------------构造待拟合函数------------------------------ fun ction f = Lorentz(p,k) f=(1-(1-p).^k(3)).^k(1).*(1-(1-p).^k(4)).^k(2);计算结果：遗传算法的初始估计数值:参数a0 = 0.942351724参数b0 = 0.889975854参数c0 = 0.655520744参数d0 = 1.049009027 最终拟合结果:残差平方和= 2.390252e-006参数 a = 1.494321551参数 b = 0.522975583参数 c = 0.714822496参数 d = 2.216764819决定系数R-Square = 0.999998738均方误差MSE = 0.000000 平均绝对误差MAE = 0.000315最大绝对误差MAS = 0.000798>> 图形结果见附图2。