最小二乘法曲线拟合以及Matlab实现
在实际工程中，我们常会遇到这种问题：已知一组点的横纵坐标，需要绘制出一条尽可能逼近这些点的曲线（或直线），以进行进一步进行加工或者分析两个变量之间的相互关系。而获取这个曲线方程的过程就是曲线拟合。
目录
最小二乘法直线拟合原理
曲线拟合
Matlab实现代码
最小二乘法直线线拟合原理
首先，我们从曲线拟合的最简单情况——直线拟合来引入问题。如果待拟合点集近似排列在一条直线上时，我们可以设直线 $y=ax+b$ 为其拟合方程，系数 $A=[a,b]$ 为待求解项，已知：
   。
用矩阵形式表达为： Y=X0AY=X_{0}AY=X0​A，其中：


要求解A，可在方程两边同时左乘 X0X_{0}X0​ 的逆矩阵，如果它是一个方阵且非奇异的话。
但是，一般情况下 X0X_{0}X0​ 连方阵都不是，所以我们在此需要用 X0X_{0}X0​ 构造一个方阵，即方程两边同时左乘 X0X_{0}X0​ 的转置矩阵，得到方程:    X0TY=X0TX0AX_{0}^{T}Y=X_{0}^{T}X_{0}AX0T​Y=X0T​X0​A  。
此时，方程的系数矩阵 X0TX0X_{0}^{T}X_{0}X0T​X0​ 为方阵，所以两边同时左乘新系数矩阵 X0TX0X_{0}^{T}X_{0}X0T​X0​ 的逆矩阵，便可求得系数向量A ，即：(X0TX0)−1X0TY=A(X_{0}^{T}X_{0})^{-1}X_{0}^{T}Y=A(X0T​X0​)−1X0T​Y=A 。
方程A=(X0TX0)−1X0TYA=(X_{0}^{T}X_{0})^{-1}X_{0}^{T}YA=(X0T​X0​)−1X0T​Y 右边各部分均已知，所以可直接求解得到拟合直线的方程系数向量A。
曲线线拟合
当样本点的分布不为直线时，我们可用多项式曲线拟合，即拟合曲线方程为n阶多项式
y=∑i=0naixi=anxn+an−1xn−1+...+a1x+a0y=\sum_{i=0}^{n}a_ix^i=a_nx^n+a_{n-1}x^{n-1}+...+a_1x+a_0y=∑i=0n​ai​xi=an​xn+an−1​xn−1+...+a1​x+a0​ 。
用矩阵形式表示为： Y=X0AY=X_0AY=X0​A ，其中：

待求解项为系数向量A=[an,an−1,...,a2,a1,a0]TA=[a_n,a_{n-1},...,a_2,a_1,a_0]^TA=[an​,an−1​,...,a2​,a1​,a0​]T。
曲线拟合方程Y=X0AY=X_0AY=X0​A 的求解方法与上面直线的求解方法一样，也是在方程Y=X0AY=X_0AY=X0​A 两边同左乘X0X_0X0​的转置矩阵得到:    X0TY=X0TX0AX_{0}^{T}Y=X_{0}^{T}X_{0}AX0T​Y=X0T​X0​A，
再同时在新方程两边同时左乘X0TX0X_{0}^{T}X_{0}X0T​X0​ 的逆矩阵，得到：(X0TX0)−1X0TY=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://blog.csdn.net/zhuoyueljl/article/details/53675460
原
最小二乘法拟合球面源代码（matlab）
2016年12月15日 20:13:56
阅读数：3200
1、源代码
[html] view plain copy
<span style="font-size:14px;">[x,y,z]=sphere(5);  
x = [0,0,1,0,0,0.7].';   %自己的数据  
y = [0,1,0,0,-1,0.7].';  %自己的数据  
z = [1,0,0,-1,0,0].';    %自己的数据  
% 拟合  
x = x(:); y = y(:); z = z(:);  
data=unique([x(:)-0.1,y(:)+0.2,z(:)],'rows');  
f=@(p,data)(data(:,1)-p(1)).^2+(data(:,2)-p(2)).^2+(data(:,3)-p(3)).^2-p(4)^2;  
p=nlinfit(data,zeros(size(data,1),1),f,[0 0 0 1]')%拟合的参数  
hold on  
plot3(data(:,1),data(:,2),data(:,3),'o')  
[X,Y,Z]=meshgrid(linspace(-14,14));  
V=(X-p(1)).^2+(Y-p(2)).^2+(Z-p(3)).^2-p(4)^2;  
isosurface(X,Y,Z,V,0);  
alpha .5;camlight;axis equal;grid on;view(3);  
title(sprintf('(x-%f)^2+(y-%f)^2+(z-%f)^2=%f',p(1),p(2),p(3),p(4)^2))</span>  

在matlab的输入窗口输入源代码后按Enter键如下图所示


2、实验结果

MATLAB最小二乘法拟合高次曲面
前言
1. **函数文件源码** ：
2. **解决上述问题**
3. **生成源代码**

前言
引用来引用去实在没意思（http://blog.sina.com.cn/s/blog_8702e2b60102x4qg.html），看到的很多最小二乘法拟合曲面方程基本都是基于这样的一个方程，代码也没有什么大的改动。只是应用的时候确实存在很多问题，不太适合实际的问题。
简单分析一下上述参考的源码存在的一些问题。
1. 函数文件源码 ：
function [a0, a1, a2, a3, a4, a5] = least_square_surface(x,y,z)
% 初始化矩阵
A = zeros(6,6);
B = zeros(6,1);
% 矩阵赋值（根据最小二乘法对最小二乘矩阵赋值）
for i=1:length(x)
    for j = 1:length(y)
        A(1,1) = 1+A(1,1);
        A(1,2) = x(i,j)+A(1,2);
        A(1,3) = y(i,j)+A(1,3);
        A(1,4) = x(i,j)^2+A(1,4);
        A(1,5) = x(i,j)*y(i,j)+A(1,5);
        A(1,6) = y(i,j)^2+A(1,6);
        
        A(2,1) = x(i,j)+A(2,1);
        A(2,2) = x(i,j)^2+A(2,2);
        A(2,3) = x(i,j)*y(i,j)+A(2,3);
        A(2,4) = x(i,j)^3+A(2,4);
        A(2,5) = x(i,j)^2*y(i,j)+A(2,5);
        A(2,6) = x(i,j)*y(i,j)^2+A(2,6);
        
        A(3,1) = y(i,j)+A(3,1);
        A(3,2) = x(i,j)*y(i,j)+A(3,2);
        A(3,3) = y(i,j)^2+A(3,3);
        A(3,4) = x(i,j)^2+A(3,4);
        A(3,5) = x(i,j)*y(i,j)^2+A(3,5);
        A(3,6) = y(i,j)^3+A(3,6);
        
        A(4,1) = x(i,j)^2+A(4,1);
        A(4,2) = x(i,j)^3+A(4,2);
        A(4,3) = x(i,j)^2*y(i,j)+A(4,3);
        A(4,4) = x(i,j)^4+A(4,4);
        A(4,5) = x(i,j)^3*y(i,j)+A(4,5);
        A(4,6) = x(i,j)^2*y(i,j)^2+A(4,6);
        
        
        A(5,1) = x(i,j)*y(i,j)+A(5,1);
        A(5,2) = x(i,j)^2*y(i,j)+A(5,2);
        A(5,3) = x(i,j)^3+A(5,3);
        A(5,4) = x(i,j)^3*y(i,j)+A(5,4);
        A(5,5) = x(i,j)^2*y(i,j)^2+A(5,5);
        A(5,6) = x(i,j)*y(i,j)^3+A(5,6);
        
        A(6,1) = y(i,j)^2+A(6,1);
        A(6,2) = x(i,j)*y(i,j)^2+A(6,2);
        A(6,3) = y(i,j)^3+A(6,3);
        A(6,4) = x(i,j)^2*y(i,j)^2+A(6,4);
        A(6,5) = x(i,j)*y(i,j)^3+A(6,5);
        A(6,6) = y(i,j)^4+A(6,6);
        
        B(1,1) = z(i,j)+B(1,1);
        B(2,1) = z(i,j)*x(i,j);
        B(3,1) = z(i,j)*y(i,j);
        B(4,1) = z(i,j)*x(i,j)^2;
        B(5,1) = z(i,j)*x(i,j)*y(i,j);
        B(6,1) = z(i,j)*y(i,j)^2;
    end
end
C = inv(A)*B;
a0 = C(1);
a1 = C(2);
a2 = C(3);
a3 = C(4);
a4 = C(5);
a5 = C(6);


源码的基本思想是基于数学方程计算 a0~a5 这6个参数的值，并将这些参数以行向量形式返回。我们观察函数文件，看到两层 for 循环的循环体中有 x(i,j) 这样的参数，这要求传入参数 x 需要是一个 length(x)*length(y) 的矩阵，而我们在实际应用的时候，可能传入参数并不满足这样的矩阵形式，更多的可能是个向量，或者[x,y,z] 这样的矩阵，直接以输入参数代入可能会报“位置 1 的索引超出数组范围(不能超过 1)。”的错误。
2. 解决上述问题
找了个比较简单的办法来解决上述问题。我用的 MATLAB 版本是 2018a ，我们知道有一个拟合工具：
>> cftool



Fit name: 自定义；

x data， y data,  z data：都需要从工作空间加载，所以你的这些数据需要先导入到 workspace（工作区）；

上部中间的部分是自己设置的参数，第一行 Polynomial 为多项式拟合，拟合算法在类 fittype里面，自己查参数意义；

当然多项式拟合可以选择 x 和 y 的最高次项的幂。
3. 生成源代码
有些同学可能在想：“我要源代码，这样才能往里面传参数啊！！！”

步骤：

文件（File）->Generate Code

它就会自动生成一个函数文件，你直接调用就行了，不过里面还有一些数据类型的问题待解决，也就是函数的返回参数：

[fitresult, gof] = createFit(x, y, Z)


可能你直接用这个做还会遇到有些问题，不过终归不是什么难题了。

还是补充上怎么获取方程系数吧，参考http://www.ilovematlab.com/thread-574266-1-1.html

MATLAB的 coeffvalues（）函数。help 文档：“coeffvalues（FUN）返回CFIT对象FUN的系数值作为行向量。”

为了很好的进行椭圆方程拟合，本文先对椭圆基本知识进行复习，后进行非标准椭圆方程拟合公式推导，最后有matlab代码的实现。
1.
用最小二乘法做椭圆拟合
1.1.
椭圆标准方程
对椭圆印象最深的就是高中时教过的，一条绳子，两个钉子，一支笔，就可以绘制出一个椭圆。固定两个钉子，让钉子之间的距离小于绳子的长度，然后用绳子的两端分别固定在两个钉子上，放一支笔在绳子的任意位置，拉紧线进行划线，画出来的完成图形就行椭圆。
然后用数学知识来解释，这两个钉子做为两个焦点，绳子的长度为2a(这里用2a，而不是a，仅是为了后续写标准的时候好看而已)，钉子之间的距离为2c。为了更好的写方程，建立一个标准方程的坐标系。以两个钉子的连线中点为坐标原点，钉子连线所在的方向为x方向，绘制xoy坐标系。
其中绳子的长度为
两个钉子(焦点)的距离为
两个焦点的坐标分别为(-c,0)和(c,0)。
然后利用绳子长度的定义进行标准方程证明，具体如下：
设椭圆上任意一点M坐标为(x，y)，则根据椭圆定义，写方程如下：
即
继续证明如下：
1.2.
椭圆非标准方程
对于椭圆的非标准方程，可以写成以下形式：
以上，是推导过程。由于是在word中编辑的文档，为了省事，这里直接截图放置在博客中，分享给大家。
以下是matlab封装的最小二乘法拟合椭圆方程的函数
function ellipse_paras =
ellipsefit(X,Y,n)
% 椭圆方程：x^2+Axy+By^2+Cx+Dy+E=0
% 采用最小二乘法进行拟合椭圆
% Input: X --- a vector of x measurements
% Y --- a vector of y measurements
% n --- the number of measurements
% author:Sigrid,2018-03-19
% M*[A B C D E]' = N
%初始化椭圆方程结果
ellipse_paras.A=0;
ellipse_paras.B=0;
ellipse_paras.C=0;
ellipse_paras.D=0;
ellipse_paras.E=0;
x=X;
y=Y;
xy=x.*y;
x2=x.*x;
y2=y.*y;
x3=x2.*x;
y3=y2.*y;
x2y=x2.*y;
xy2=x.*y2;
x4=x2.*x2;
y4=y2.*y2;
x3y=x3.*y;
xy3=x.*y3;
x2y2=x2.*y2;
% Construct M
M=[ sum(x2y2) sum(xy3) sum(x2y) sum(xy2) sum(xy);
sum(xy3) sum(y4) sum(xy2) sum(y3) sum(y2);
sum(x2y) sum(xy2) sum(x3) sum(xy) sum(x) ;
sum(xy2) sum(y3) sum(xy) sum(y2) sum(y) ;
sum(xy) sum(y2) sum(x) sum(y) n ];
% Construct N
N = -[sum(x3y) sum(x2y2) sum(x3) sum(x2y) sum(x2)]';
P=M\N;
A=P(1);
B=P(2);
C=P(3);
D=P(4);
E=P(5);