往往在实际问题中都存在exp(x)、lnx、sinx等多种函数组合的非线性经验公式。对此我们就可以通过lsqcurvefit函数进行求解,该函数的方法被称为非线性最小二乘,损失函数一样,只不过类似于优化算法,给定参数初始值,然后优化参数,非线性最小二乘模型如下,即目标函数。
lsqcurvefit函数拟合格式
格式
x = lsqcurvefit(fun,a0,xdata,ydata)
x = lsqcurvefit(fun,a0,xdata,ydata,lb,ub)
x = lsqcurvefit(fun,a0,xdata,ydata,lb,ub,options)
[x,resnorm] = lsqcurvefit(…)
[x,resnorm,residual] = lsqcurvefit(…)
[x,resnorm,residual,exitflag] = lsqcurvefit(…)
[x,resnorm,residual,exitflag,output] = lsqcurvefit(…)
[x,resnorm,residual,exitflag,output,lambda]= lsqcurvefit(…)
[x,resnorm,residual,exitflag,output,lambda,jacobian] =lsqcurvefit(…)
参数说明:
在lsqcurvefit函数中,有trust-region-reflective和levenberg-marquardt两种算法可以求解参数,
其中默认采用信赖域(trust-region-reflective)算法。
a0为初始解向量, 因为求解是一个迭代的过程,需要先给定一个初始参数,再逐步修改参数的过程。
所以要对a0初始化,一般而言,可以随机,但是经验上取与解接近的值会提高计算速度。
xdata,ydata为满足关系ydata=F(a, xdata)的数据;
lb、ub为解向量的下界和上界lb≤a≤ub,若没有指定界,则lb=[ ],ub=[ ];
options为指定的优化参数;
fun为待拟合函数,计算x处拟合函数值,其定义为 function F = myfun(a,xdata)
resnorm=sum ((fun(a,xdata)-ydata).^2),即在a处残差的平方和;
residual=fun(a,xdata)-ydata,即在x处的残差;
exitflag为终止迭代的条件;
output为输出的优化信息;
lambda为解x处的Lagrange乘子;
jacobian为解x处拟合函数fun的jacobian矩阵。
f:符号函数句柄,如果是以m文件的形式调用的时候,别忘记加@.这里需要注意,f函数的返回值是和y
匹对的,即拟合参数的标准是(f-y)^2取最小值,具体看下面的例子
实例1
程序
clc;clear%清除变量
y=[100.3 101.1 102.1 101.1 101.6 104.4 102.5 102.1 103.9 103.9];
xdata=1:length(y);
a0=[1,3,7,5,7]; %初始估计值,随便写 这个是4次拟合 ,具体表达式可以随便改
options=optimset('Tolfun',1e-15); %方法设定
for i=1:1000
x=lsqcurvefit(@fun1,a0,xdata,y,[],[],options); %确定待定系数
a0=x;%以计算出的x为初值,循环迭代1000次
end
disp(x);%输出系数
yy = fun1(x,xdata);%利用已经拟合好的模型预测y值
len=[1:20];
len1 = fun1(x,len);%预测走势
result=[y;yy]%实际值与预测值
error=abs(y-yy);%误差
bfb=error./y%相对误差
errorsum=sum(error)/length(y)%平均误差
bfbsum=sum(bfb)/length(y)%平均相对误差
figure(1)
plot(xdata,y,'r-',xdata,yy,'b-')
legend('实际值','拟合值')
title('实际值与预测值的比较','fontsize',15)
ylabel('Y','fontsize',15)
xlabel('X','fontsize',15)
%axis([1 10 1 200]);%坐标范围
%set(gca,‘xtick’,[1 2 …10]);%X轴设定
%set(gca,‘ytick’,[1 20 40 …200]);%Y轴设定
figure(2)
plot(len,len1,'r-')
legend('拟合曲线')
title('拟合曲线图','fontsize',15)
ylabel('Y','fontsize',15)
xlabel('X','fontsize',15)
function y = fun1(a, x);
y = a(1) + a(2)*x + a(3)*x.^2 + a(4)*x.^3 + a(5)*x.^4;
end
运行结果
实例2
程序
clc;
clear all;
close all;
xdata = ...
[0.9 1.5 13.8 19.8 24.1 28.2 35.2 60.3 74.6 81.3];
ydata = ...
[455.2 428.6 124.1 67.3 43.2 28.1 13.1 -0.4 -1.3 -1.5];
fun = @(x,xdata)x(1)*exp(x(2)*xdata);
x0 = [100,-1];
x = lsqcurvefit(fun,x0,xdata,ydata)
times = linspace(xdata(1),xdata(end));
plot(xdata,ydata,'ko',times,fun(x,times),'b-')
legend('Data','Fitted exponential')
title('Data and Fitted Curve')
运行结果
实例3
用MATLAB程序拟合Logistic函数
程序
clear
clc
xdata=0:10:180;
ydata=[0 0 0.45 2.7 5.4 5.7 10.5 10.8 9.6 12.15 16.65 18.15 19.05 28.2 29.1 21.1 19.95 22.05 25.2];
%% 指定非线性函数拟合曲线
X0=[100 10 0.2];
[parameter,resnorm]=lsqcurvefit(@fun,X0,xdata,ydata); %指定拟合曲线
A=parameter(1);
B=parameter(2);
C=parameter(3);
fprintf('拟合Logistic曲线的参数A为:%.8f,B为:%.8f,C为:%.8f', A, B, C);
fit_y=fun(parameter,xdata);
figure(1)
plot(xdata, ydata, 'r*');
hold on
plot(xdata,fit_y,'b-');
xlabel('t');
ylabel('y');
legend('观测数据点','拟合曲线', 'Location', 'northwest');
saveas(gcf,sprintf('Logistic曲线.jpg'),'bmp');
%% Logistic函数
% y=A/(1+B*exp(-C*t))
function f=fun(X,t)
f=X(1)./(1+X(2).*exp(-X(3).*(t)));
end
运行结果