MATLAB的lsqcurvefit函数拟合非线性式子中的未知参数

往往在实际问题中都存在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

运行结果