最小二乘参数估计(附完整MATLAB代码)

最小二乘法是一种广泛使用的参数估计方法。其基本思想是:给定一组数据点(x1,y1),(x2,y2),…,(xn,yn),通过找到一组参数,使实际观测到的y值和使用这组参数预测的y值之间的误差平方和最小化。

具体来说,我们假设数据符合以下模型:

yi = f(xi,β) + εi

这里f(xi,β)是一个含有未知参数β的函数,描述了y和x之间的关系;εi是随机误差。

我们要估计的参数就是β。最小二乘估计的目标函数是:

min Σ(yi - f(xi,β))2

通过求解这个最优化问题,我们可以得到参数β的估计值。

举一个简单的例子,如果我们假设关系为线性的:

yi = β1 + β2*xi + εi

那么问题就变成了线性回归,我们可以通过求解正常方程或使用梯度下降算法估计β1和β2。

总的来说,最小二乘法通过最小化预测值和实际值的误差来寻找最佳的参数估计,这使得它成为统计分析中一种非常有用和常被使用的方法。

MATLAB代码如下:

clc;clear all;close all;
%% 最小二乘非线性函数参数估计
% 非线性曲线拟合是已知输入向量xdata和输出向量ydata，
% 并且知道输入与输出的函数关系为ydata=F(x, xdata)，
% 但不知道系数向量x。今进行曲线拟合，求x使得输出的如下最小二乘表达式成立：
% min Σ(F(x,xdatai)-ydatai)^2
%% (1)输入数据
data=[2.85     -13.438017 
2.90     -13.988932 
2.95     -14.443324 
3.00     -14.810092 
3.05     -15.096124 
3.10     -15.308248 
3.15     -15.452979 
3.20     -15.536750 
3.25     -15.565603 
3.30     -15.544858 
3.35     -15.480221 
3.40     -15.376477 
3.45     -15.238199 
3.50     -15.069909 
3.55     -14.875819 
3.60     -14.659466 
3.70     -14.172289 
3.80     -13.631630];
%% (2)数据预处理
xdata =data(:,1).^3';%自变量
ydata =data(:,2)';%因变量

%% (3)开始拟合
x0 = [2
    -1
    1
    4];%设置初始估计值
LB=[-1000;
    -1000
    -1000
    -1000];%设置估计值的下限
UB=[1000;
    1000
    1000
    1000];%设置估计值的上限
[x,resnorm] = lsqcurvefit(@myfun,x0,xdata,ydata,LB,UB);%最小二乘法曲线拟合
V0=x(1);%获取V0的值
E0=x(2);%获取E0的值
B0=x(3)*160.2;%获取B0的值
B01=x(4);%获取B01的值
%% 输出结果
['V0=',num2str(V0)]
['E0=',num2str(E0)]
['B0=',num2str(B0)]
['B01=',num2str(B01)]

%% 测试拟合结果
Eoutput=myfun(x,xdata);%用Vi输入参数方程计算拟合的Ei,参数用 lsqcurvefit计算出来的值
x=x0;%参数用初值
Eoutput2=myfun(x,xdata);%用Vi输入参数方程计算拟合的Ei

%% 绘图
figure;
plot(xdata,ydata,'bo-');
xlabel('Vi');
ylabel('Ei');
legend('实际曲线','拟合曲线');
title('待拟合曲线');

figure;
plot(xdata,ydata,'bo-');
hold on;
plot(xdata,Eoutput2,'r*-');
xlabel('Vi');
ylabel('Ei');
legend('实际曲线','拟合初值的曲线');
title('拟合初值的拟合曲线与待拟合曲线比较');

figure;
plot(xdata,ydata,'bo-');
hold on;
plot(xdata,Eoutput,'r*-');
xlabel('V');
ylabel('E');
legend('实际曲线','拟合曲线');

%% 失败的例子
x0 = [-2
   1
    1
   2];%初始估计值
LB=[-1000;
    -1000
    -1000
    -1000];%下限
UB=[1000;
    1000
    1000
    1000];%上限
[x,resnorm] = lsqcurvefit(@myfun,x0,xdata,ydata,LB,UB);%最小二乘法曲线拟合
xdata=xdata;
Eoutput3=myfun(x,xdata);
x=x0;
Eoutput2=myfun(x,xdata);
figure;
plot(xdata,ydata,'bo-');
hold on;
plot(xdata,Eoutput2,'r*-');
xlabel('Vi');
ylabel('Ei');
legend('实际曲线','拟合初值的曲线');
title('拟合初值的拟合曲线与待拟合曲线比较');

figure;
plot(xdata,ydata,'bo-');
hold on;
plot(xdata,Eoutput3,'r*-');
xlabel('V');
ylabel('E');
legend('实际曲线','拟合曲线');
title('拟合曲线与待拟合曲线比较(失败)');

function E = myfun(x,xdata)
%% 函数表达式
%V0=体积
%E0=基态能量
%B0=体积模量
%B01=体积模量一撇
V0=x(1);
E0=x(2);
B0=x(3);
B01=x(4);
V=xdata;%自变量
E=E0+B0*V./B01.*((V0./V).^B01./(B01-1)+1)-B0*V0/(B01-1);%方程

程序结果如下:

V0=47.3058

E0=-15.5485

B0=72.0815

B01=3.4794