曲线拟合实例

例如：

已知数据队列buf=【5410。。。。。】

x取值1：n n是队列长度

函数f（x）=a+b*sin(c*x+d).

avg是队列平均值

a b c d为参数a范围（2/3,1）*avg

b范围（0,1/3）*avg

c的范围（0,24*pi）

d（0,2*pi）

1、首先定义目标函数

function y=ga_curfit(x)

global ydata n

t=1:n;

y=0;

for i=1:n

y=y+(ydata(i)-(x(:,1)+x(:,2).*sin(x(:,3).*t(i)+x(:,4)))).^2/n; end

y=sqrt(y);

end

2、把数据b.txt放在工作空间目录中

然后再命令窗口中输入

clear

global ydata n

format long g

load b.txt

ydata=b';

n=length(ydata);

avg=sum(ydata)/n;

LB=[2/3*avg000];

UB=[1*avg1/3*avg24*pi2*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 on

plot(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;clc

format long data=;

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)