最小二乘法曲线拟合以及Matlab实现

最小二乘法曲线拟合以及Matlab实现最小二乘法曲线拟合以及Matlab实现在实际工程中,我们常会遇到这种问题:已知一组点的横纵坐标,需要绘制出一条尽可能逼近这些点的曲线(或直线),以进行进一步进行加工或者分析两个变量之间的相互关系。而获取这个曲线方程的过程就是曲线拟合。目录最小二乘法直线拟合原理曲线拟合Matlab实现代码最小二乘法直线线拟合原理首先,我们从曲线拟合的最简单情况——直线拟合来引……

大家好,又见面了,我是你们的朋友全栈君。

最小二乘法曲线拟合以及Matlab实现

在实际工程中,我们常会遇到这种问题:已知一组点的横纵坐标,需要绘制出一条尽可能逼近这些点的曲线(或直线),以进行进一步进行加工或者分析两个变量之间的相互关系。而获取这个曲线方程的过程就是曲线拟合。

目录

  • 最小二乘法直线拟合原理
  • 曲线拟合
  • Matlab实现代码

最小二乘法直线线拟合原理

首先,我们从曲线拟合的最简单情况——直线拟合来引入问题。如果待拟合点集近似排列在一条直线上时,我们可以设直线 y = a x + b y=ax+b y=ax+b 为其拟合方程,系数 A = [ a , b ] A=[a,b] A=[a,b] 为待求解项,已知:
这里写图片描述 。

用矩阵形式表达为: Y = X 0 A Y=X_{0}A Y=X0​A,其中:
这里写图片描述
要求解A,可在方程两边同时左乘 X 0 X_{0} X0​ 的逆矩阵,如果它是一个方阵且非奇异的话。

但是,一般情况下 X 0 X_{0} X0​ 连方阵都不是,所以我们在此需要用 X 0 X_{0} X0​ 构造一个方阵,即方程两边同时左乘 X 0 X_{0} X0​ 的转置矩阵,得到方程: X 0 T Y = X 0 T X 0 A X_{0}^{T}Y=X_{0}^{T}X_{0}A X0T​Y=X0T​X0​A 。

此时,方程的系数矩阵 X 0 T X 0 X_{0}^{T}X_{0} X0T​X0​ 为方阵,所以两边同时左乘新系数矩阵 X 0 T X 0 X_{0}^{T}X_{0} X0T​X0​ 的逆矩阵,便可求得系数向量A ,即: ( X 0 T X 0 ) − 1 X 0 T Y = A (X_{0}^{T}X_{0})^{-1}X_{0}^{T}Y=A (X0T​X0​)−1X0T​Y=A 。

方程 A = ( X 0 T X 0 ) − 1 X 0 T Y A=(X_{0}^{T}X_{0})^{-1}X_{0}^{T}Y A=(X0T​X0​)−1X0T​Y 右边各部分均已知,所以可直接求解得到拟合直线的方程系数向量A。

曲线线拟合

当样本点的分布不为直线时,我们可用多项式曲线拟合,即拟合曲线方程为n阶多项式

y = ∑ i = 0 n a i x i = a n x n + a n − 1 x n − 1 + . . . + a 1 x + a 0 y=\sum_{i=0}^{n}a_ix^i=a_nx^n+a_{n-1}x^{n-1}+…+a_1x+a_0 y=∑i=0n​ai​xi=an​xn+an−1​xn−1+...+a1​x+a0​ 。

用矩阵形式表示为: Y = X 0 A Y=X_0A Y=X0​A ,其中:

X0

待求解项为系数向量 A = [ a n , a n − 1 , . . . , a 2 , a 1 , a 0 ] T A=[a_n,a_{n-1},…,a_2,a_1,a_0]^T A=[an​,an−1​,...,a2​,a1​,a0​]T。

曲线拟合方程 Y = X 0 A Y=X_0A Y=X0​A 的求解方法与上面直线的求解方法一样,也是在方程 Y = X 0 A Y=X_0A Y=X0​A 两边同左乘 X 0 X_0 X0​的转置矩阵得到: X 0 T Y = X 0 T X 0 A X_{0}^{T}Y=X_{0}^{T}X_{0}A X0T​Y=X0T​X0​A,

再同时在新方程两边同时左乘 X 0 T X 0 X_{0}^{T}X_{0} X0T​X0​ 的逆矩阵,得到: ( X 0 T X 0 ) − 1 X 0 T Y = A (X_{0}^{T}X_{0})^{-1}X_{0}^{T}Y=A (X0T​X0​)−1X0T​Y=A

上式左边各部分均已知,所以可直接求解得拟合曲线方程的系数向量A。

Matlab实现代码

%by hanlestudy@163.com
clear
clc
x=[2,4,5,6,6.8,7.5,9,12,13.3,15];
y=[-10,-6.9,-4.2,-2,0,2.1,3,5.2,6.4,4.5];
[~,k]=size(x);
for n=1:9
    X0=zeros(n+1,k);
    for k0=1:k           %构造矩阵X0
        for n0=1:n+1
            X0(n0,k0)=x(k0)^(n+1-n0);
        end
    end
    X=X0';
    ANSS=(X'*X)\X'*y';
    for i=1:n+1          %answer矩阵存储每次求得的方程系数,按列存储
       answer(i,n)=ANSS(i);
   end
    x0=0:0.01:17;
    y0=ANSS(1)*x0.^n    ;%根据求得的系数初始化并构造多项式方程
    for num=2:1:n+1     
        y0=y0+ANSS(num)*x0.^(n+1-num);
    end
    subplot(3,3,n)
    plot(x,y,'*')
    hold on
    plot(x0,y0)
end
suptitle('不同次数方程曲线拟合结果,从1到9阶')

运行结果

拟合曲线结果:

可以看出看来,当多项式的阶数过小是,曲线并不能很好地反映出样本点的分布情况;但阶数过高时,会出现过拟合的情况。
这里写图片描述

系数矩阵answer:

这里写图片描述

Matlab自带函数——polyfit

在matlab中,也有现成的曲线拟合函数polyfit,其也是基于最小二乘原理实现的,具体用法为:ans=polyfit(x,y,n). 其中x,y为待拟合点的坐标向量,n为多项式的阶数。

下面代码是用polyfit函数来做曲线拟合:

clear
clc
x=[2,4,5,6,6.8,7.5,9,12,13.3,15];
[~,k]=size(x);
y=[-10,-6.9,-4.2,-2,0,2.1,3,5.2,6.4,4.5];
for n=1:9
    ANSS=polyfit(x,y,n);  %用polyfit拟合曲线
    for i=1:n+1           %answer矩阵存储每次求得的方程系数,按列存储
       answer(i,n)=ANSS(i);
   end
    x0=0:0.01:17;
    y0=ANSS(1)*x0.^n    ; %根据求得的系数初始化并构造多项式方程
    for num=2:1:n+1     
        y0=y0+ANSS(num)*x0.^(n+1-num);
    end
    subplot(3,3,n)
    plot(x,y,'*')
    hold on
    plot(x0,y0)
end
suptitle('不同次数方程曲线拟合结果,从1到9阶')

运行结果:

用polyfit拟合的结果与第一份代码运行的结果基本一样
这里写图片描述

申明

本文为本人原创,转载请注明出处!

版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 举报,一经查实,本站将立刻删除。

发布者:全栈程序员-站长,转载请注明出处:https://javaforall.net/152976.html原文链接:https://javaforall.net

赞 (0)
全栈程序员-站长的头像全栈程序员-站长


相关推荐

  • offset size_c语言sizeof求结构体长度

    offset size_c语言sizeof求结构体长度sizeof与offsetof在程序中经常遇到,但在面试中其应用使得许多小伙伴吃闭门羹,被面试官问得哑口无言。接下来对两者的应用做详细介绍。关于sizeof定义:sizeof乃C/C++中的一个操作符(operator),简单的说其作用就是返回一个对象或者类型所占的内存字节数。(MSDN上的解释为:Thesizeofkeywordgivestheamountofstorage,in

    2022年8月22日
  • DELPHI程序员招聘_招聘java程序员

    DELPHI程序员招聘_招聘java程序员北京地区招聘Delphi程序员,要求工作经验2年以上,熟悉Delphi7+SQL有PB开发经验优先可全职或外派工作地点:西城区六铺炕联系QQ:408390946

    2022年8月1日
  • springaop实现原理面试_springmvc模式的工作原理

    springaop实现原理面试_springmvc模式的工作原理SpringAop实现原理什么是AOPAOP(Aspect-OrientedProgramming,面向方面编程),可以说是OOP(Object-OrientedPrograming,面向对象编程)的补充和完善。OOP引入封装、继承和多态性等概念来建立一种对象层次结构,用以模拟公共行为的一个集合。当我们需要为分散的对象引入公共行为的时候,OOP则显得无能为力。也就是说,OOP允许你定义从上到下的

    2022年9月18日
  • platform device

    platform deviceplatformdevice================= 头文件:linux/platform_device.h  为什么使用platformdevice?————————–    从硬件的角度来说,集成在嵌入式芯片内部的外设离CPU最近,它们不依附于GPIO或者PCI,I2C此类的 总线,它们挂接在soc内存空

    2022年7月24日
  • 原地算法矩阵置0_矩阵归一化处理步骤

    原地算法矩阵置0_矩阵归一化处理步骤给定一个 m x n 的矩阵,如果一个元素为 0 ,则将其所在行和列的所有元素都设为 0 。请使用 原地 算法。进阶:一个直观的解决方案是使用 O(mn) 的额外空间,但这并不是一个好的解决方案。一个简单的改进方案是使用 O(m + n) 的额外空间,但这仍然不是最好的解决方案。你能想出一个仅使用常量空间的解决方案吗?示例 1:输入:matrix = [[1,1,1],[1,0,1],[1,1,1]]输出:[[1,0,1],[0,0,0],[1,0,1]]示例 2:输入:matrix

    2022年8月8日
  • CLion和vscode快捷键设置[通俗易懂]

    CLion和vscode快捷键设置[通俗易懂]ctrl+shift+F7:把某个单词高亮(),ESC:退出高亮。

    2022年6月11日

发表回复

您的邮箱地址不会被公开。 必填项已用 * 标注

关注全栈程序员社区公众号