MATLAB之最小二乘法

一、算法原理

给定一些列点x1,x2,.....xn,对应的函数值为y1,y2,......yn。若拟合曲线为y=ax+b,根据条件可写出线形方程组为:

[x1 1;x2 1;...;x3;1]*[a;b]=[y1;...;yn]即A*[a;b]=[y2;...;yn]

因为A不是方阵,无法求逆。故做如下变形:

方法一:

\Rightarrow A'*A*[a;b]=A'*[y1;...;yn]

\Rightarrow [a;b]=inv(A'*A)*A'*[y1,...,yn]

方法二:

对矩阵A做QR分解A=Q*R

\Rightarrow Q*R*W=[y1,...,yn]

\Rightarrow R*w=Q'*[y1,...,yn]

\Rightarrow w=inv(R)*Q'*[y1,...,yn]

二、matlab程序

%% 最小二乘法插值
%设拟合直线为y=ax+b   ax1+b=y1...axn+b=yn
%写成[x1 1;x2 1;...;x3;1]*[a;b]=[y2;...;yn]  即A*[a;b]=[y2;...;yn]
clear
clc
close
x=1:6; 
y=[1 4 5 8 10 11];%一系列插值点
plot(x,y,'o');%绘制散点图
hold on
axis([-5 10 -2 15]);%坐标轴范围

A=[1 1;2 1;3 1;4 1;5 1;6 1];%A*[a;b]=[y2;...;yn]
y=y';
% y=[1;4;5;8;10;11];  %法方方程组 A*W=Y  即 A'*A*W=A'*Y  (
w=inv(A'*A)*A'*y;      %W=inv(A'*A)*A'*Y   求得系数W=[a,b]
y1=w(1).*x+w(2);
plot(x,y1,':');
% QR分解法做最小二乘拟合
Q=orth(A);  %QR分解法  A*W=Y  A=Q*R   Q*R*W=Y  W=inv(R)*Q'*Y
R=Q'*A;
a1=inv(R)*Q'*y;
x=0:8;
y2=a1(1).*x+a1(2);
plot(x,y2);

三、matlab自带最小二乘拟合函数  lsqcurvefit

1、函数简介

 x = lsqcurvefit(fun,a0,xdata,ydata,lb,ub,options)

非线性曲线拟合是已知输入向量xdata和输出向量ydata,
并且知道输入与输出的函数关系为ydata=F(x, xdata),即已知拟合曲线的形式(一次函数,二次函数,...)
但不知道系数向量a。进行曲线拟合,求x使得输出的最小二乘表达式成立:
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(…)
参数说明:
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),即在x处残差的平方和;
residual=fun(a,xdata)-ydata,即在x处的残差;
exitflag为终止迭代的条件;
output为输出的优化信息;
lambda为解a处的Lagrange乘子;
jacobian为解a处拟合函数fun的jacobian矩阵。

2、matlab程序

close
clc
clear
xdata=linspace(0,2*pi,15);
y=5*sin(xdata)+2*xdata+xdata.^2;
y=y+2*rand(1,15);
plot(xdata,y,'o') %绘制出原始数据的散点图
hold on

fun=@(x,xdata) x(1)*sin(xdata)+x(2)*xdata+x(3)*xdata.^2;%待求拟合函数的形式
x=lsqcurvefit(fun,[0 0 0],xdata,y);% [0 0 0]为插值多项式初始系数a0,
%xdata为输入原始数据x坐标,y为原始数据纵坐标,返回值x为待求系数矩阵
xx=linspace(0,2*pi,150);%拟合函数曲线的x坐标
yy=fun(x,xx);%拟合函数曲线
plot(xx,yy)%绘制拟合函数曲线

lb=[-1 -1 -1];%系数下限
ub=[6 3 2]; %系数上限
x=lsqcurvefit(fun,[0 0 0],xdata,y,lb,ub);
xx=linspace(0,2*pi,150);
yy=fun(x,xx);
plot(xx,yy)
%选项设置
%options = optimset(Name,Value,...)不可以将求解器名称作为第一个参数
%options = optimoptions(SolverName,Name,Value,...)可以将求解器名称作为第一个参数
% options=optimoptions('lsqcurvefit','Algorithm','levenberg-marquardt');%算法设置
options=optimoptions('lsqcurvefit','Display','final');%显示设置
lb=[-1 -1 -1];
ub=[6 3 2];
[x,~,~,exitflag,~,~,jacobian]=lsqcurvefit(fun,[0 0 0],xdata,y,lb,ub,options)
xx=linspace(0,2*pi,150);
yy=fun(x,xx);
plot(xx,yy)

  • 13
    点赞
  • 89
    收藏
  • 打赏
    打赏
  • 1
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
©️2022 CSDN 皮肤主题:1024 设计师:我叫白小胖 返回首页
评论 1

打赏作者

天涯铭

你的鼓励将是我创作的最大动力

¥2 ¥4 ¥6 ¥10 ¥20
输入1-500的整数
余额支付 (余额:-- )
扫码支付
扫码支付:¥2
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值