拟合

一、曲线拟合的线性最小二乘法

（一）线性最小二乘法

1. 问题提出

（1）简介

已知一组（二维）数据，即平面上的n个点$(x_i,y_i)\ (i=1,2,\cdots,n)$，$x_i$互不相同，求一个函数（曲线）$y=f(x)$使得$f(x)$在某种准则下与所有数据点最为接近，即曲线拟合的最好

线性最小二乘法是解决曲线拟合最常见的方法，思路为：令

$f(x) = a_1r_1(x)+a_2r_2(x)+\cdots+a_mr_m(x)$

其中$r_k(x)$是事先选定的一组线性无关的函数，$a_k$是待定系数 $(k=1,2,\cdots,m\quad ,m<n)$

拟合的准则为使$y_i\ ,(i=1,2,\cdots,n)$与$f(x_i)$的距离$\delta_i$的平方和最小，称为最小二乘准则

（2）系数$a_k$的确定

记

$J(a_1,\cdots,a_m)=\sum^n_{i=1}\delta^2_i=\sum^n_{i=1}[f(x_i)-y_i]^2$

为求$a_1,\cdots,a_m$使得$J$达到最小，秩序利用极值的必要条件$\dfrac{\partial J}{\partial A_j}=0\ (j=1,\cdots,m)$，得到关于$a_1,\cdots,a_m$的线性方程组

$\sum^n_{i=1}r_j(x_i)[\sum^m_{k=1}a_kr_k(x_i)-y_i]=0\quad,(j=1,\cdots,m)$

即

$\sum^n_{i=1}a_k[\sum^m_{k=1}r_j(x_i)r_k(x_i)]=\sum^n_{i=1}r_j(x_i)y_i\quad,(j=1,\cdots,m)$

记

$R= \begin{bmatrix} r_1(x_1) & \cdots & r_m(x_1)\\ \vdots & \vdots & \vdots \\ r_1(x_n) & \cdots & r_m(x_n)\\ \end{bmatrix}_{n\times m}\\\\ A=[a_1,\cdots,a_m]^{\rm T}\ ,\ Y=[y_1,\cdots,y_n]^{\rm T}$

当${r_1(x),\cdots,r_m(x)}$线性无关时，$R$列满秩，$R^{\rm T}R$可逆，于是上方程组有唯一解

$A = (R^{\rm T}R)^{-1}R^{\rm T}Y$

（3）函数$r_k(x)$的选择

通常通过机理分析能够知道有什么样的函数关系，如果不清楚其关系，可以先对数据作图，直观的判断应该使用的拟合曲线

常见的拟合曲线选择

直线 $y=a_1x+a_2$
多项式（一般m=2，3，不宜太高）$y=a_1x^m+\cdots+a_mx+a{\small m+1}$
双曲线（一支）$y=\dfrac{a_1}{x}+a_2$
指数曲线 $y=a_1e^{a_2x}$

对于指数曲线，拟合前需做变量代换，化为对$a_1,a_2$的线性函数

可以选择多种曲线进行拟合，通过判断最后$J$值，选择最小的一条

也可以通过计算MSE、MAE等数值对拟合的结果进行评价

（二）MATLAB实现最小二乘法

1. 解方程法

由上述的可得 $J(a_1,\cdots,a_m) = ||RA-Y||^2$

MATLAB中的线性最小二乘标准型为 $\min\limits_A||RA-Y||^2_2$

CoreCode : A = R/Y

在此方法下，需要将用于训练的x0，y0数据初始化为列向量

其中A是返回的系数矩阵，R为拟合所用函数由x0列向量对应值构建的矩阵，Y为训练用的y0列向量

R的具体的构建方法如下：

将构建函数对应的每一项依次带入列向量x0，一个系数一列，所以需要将x0，y0数据初始化为列向量

特殊的，常数项可以用ones(size(x0))来替换，同时也可以用于倒数等位置Init用1的情况

最终得出的A数组每一位对应R矩阵的每一列所求的系数

Eg：

拟合函数：$a+bx^2+c\dfrac{1}{x}$ ，已知训练数据 $x0\ ,y0$ ，需要进行拟合的范围 $x$

R = [ones(size(x0)), x0.^2, ones(size(x0))./x]；
a = R(1)
b = R(2)
c = R(3)
y = a + b*x.^2 + c./x % 输出x所包含范围内的拟合y值

2. 多项式拟合Poly

如果取 ${r1(x),\cdots,r{m+1}(x)}={1,x,\cdots,x^m}$，即用m次多项式拟合给定的数据，可以用MATLAB中的poly系列函数执行

（1）拟合函数

Code : A = polyfit(x0, y0, <Count>)

其中<Count>为拟合多项式的次数，输出的参数A为拟合多项式的系数数组，从高次项向低次项的顺序

即在拟合函数$y=a(1)x^m+\cdots+a(m)x+a(m+1)$中，$a=[a(1),\cdots,a(m),a(m+1)]$

（2）数值计算

Code : y = polyval(A,x)

其中x为拟合需求的取值，当x换成x0的时候，可以用与求y_fit_test，并进行MSE检测

3. MSE检测拟合度

MSE(Mean Squared Error)是指最小平方误差，可以检验拟合结果和实际结果的MSE来判断拟合的程度

（1）MATLAB实现代码：

MSEValue = mse(<y_real> - <y_fit_test>)

如果使用解方程法拟合，应当首先对<y_real>进行转置，即 mse(<y_real>'-<y_fit_test>)

（2）y_fit_test的获取方法

拟合相应点
- 解方程法：y_fit_test = f(x0)
  - f为代入系数矩阵A的相应的函数求解表达式，其中自变量改为x0
- polyval法：y_fit_test = polyval(r,x0)
  - r为polyfit所得到的系数数组，进行ployval的自变量为x0
包含拟合后索引（x的步长需为等长）
- y_fit_test = y(1+(x0-min(x)./step))
  - 其中x为拟合要求的数据范围，step是x的步长
- y_fit_test = y(round(1+(x0-min(x))./(x(2)-x(1))))