Matlab之数据回归分析 --- 最小二乘估计曲线拟合

1、最小二乘估计的概念

“二乘”意即平方的含义，所以这里也可以称为“最小平方估计”拟合，那么也就有了谁的平方的问题了，直观上理解就是测量值与统计估计值的偏差，工程上又叫残余误差。而“最小”意即所有的残余误差的平方和最小！

通过下图更直观理解，就是让所有红线分别平方，然后求和，所得的值最小！

2、Matlab直线拟合

2.1、直线拟合的数学推导

假设现有一组数据【x1，y1】【x2，y2】…【xn，yn】

设其拟合直线表达式为：$$y=a+bx\text{\hspace{1em}(1)}$$

对应残余误差表达式为：$$d_i=y_i-(a+b{x}_i)\text{\hspace{1em}(2)}$$

最小二乘登场：

$$D_i=\sum _{i=1}^n{d}_i^2=\sum _{i=1}^n(y_i-a-b{x}_i{)}^2\text{\hspace{1em}(3)}$$

因为使得平方和最小，意即求导为0，那么我们就分别对该式子的系数$a、b$求偏导可得到如下等式：

$$\{\begin{array}{l}\frac{\partial D_i}{\partial a}=\sum _{i=1}^{n}2(y_i-a-b{x}_i)(-1)=-2(\sum _{i=1}^n{y}_i-na-b\sum _{i=1}^n{x}_i)=0\\ \\ \frac{\partial D_i}{\partial b}=\sum _{i=1}^{n}2(y_i-a-b{x}_i)(-x_i)=-2(\sum _{i=1}^n{x}_i{y}_i-a\sum _{i=1}^n{x}_i-b\sum _{i=1}^n{x}_i^2)=0\end{array}\text{\hspace{1em}(4)}$$

整理可得，

$$\{\begin{array}{l}(\sum _{i=1}^n{y}_i-na-b\sum _{i=1}^n{x}_i)=0\\ \\ \sum _{i=1}^n{x}_i{y}_i-a\sum _{i=1}^n{x}_i-b\sum _{i=1}^n{x}_i^2)=0\end{array}⟹\{\begin{array}{l}a=\overline{y}-b\overline{x}\\ \\ b=\frac{\overline{xy}-\overline{x}\overline{y}}{\overline{x^2}-(\overline{x}{)}^2}\end{array}\text{\hspace{1em}(5)}$$

2.2、Matlab代码

clc;
clear;

%录入X轴数据
for a = 1:30 
    x(a) = a-1;
end
%录入Y轴数据
y=[1,2,3,8,6,9,5,4,8,5,9,19,16,12,15,24,22,36,40,40,32,32,36,39,52,52,56,57,62,69];

figure;
plot(x,y,'.');%画点
hold on
b = ( mean(x*y(:)) - mean(x(:)).*mean(y(:)) ) / (mean(x*x(:))-mean(x(:))^2 );
a = mean(y(:)) - b*mean(x(:));
y1 = a+b*x;
plot(x,y1);

2.3、绘制效果

3、多项式曲线拟合

3.1、曲线拟合的数学推导

同时，假设有一组数据【x1，y1】【x2，y2】…【xn，yn】

设拟合多项式表达式为：
$$y=a+a_{1}x+...+a_k{x}^k\text{\hspace{1em}(1)}$$

残余误差和表达式为：

$$D^2=\sum _{i=1}^n[y_i-(a_0+a_1{x}_i+...+a_k+x_i^k){]}^2\text{\hspace{1em}(2)}$$

为求平方和最小，我们只需对系数$a_0{a}_k$挨个求偏导，即可得到：

$$\{\begin{array}{l}-2\sum _{i=1}^n[y_i-(a_0+a_1{x}_i+...+a_k{x}_i^k)]=0\\ -2\sum _{i=1}^n[y_i-(a_0+a_1{x}_i+...+a_k{x}_i^k)]x_i=0\\ ......\\ -2\sum _{i=1}^n[y_i-(a_0+a_1{x}_i+...+a_k{x}_i^k)]x_i^k=0\end{array}\text{\hspace{1em}(3)}$$

上述等式继续化简可得：

$$\{\begin{array}{l}{a}_{0}n+a_1\sum _{i=1}^n{x}_i+...+a_k\sum _{i=1}^n{x}_i^k=\sum _{i=1}^n{x}_i^0{y}_i\\ a_0\sum _{i=1}^n{x}_i+a_1\sum _{i=1}^n{x}_i^2+...+a_k\sum _{i=1}^n{x}_i^{k+1}=\sum _{i=1}^n{x}_i^1{y}_i\\ ...\\ a_0\sum _{i=1}^n{x}_i^k+a_1\sum _{i=1}^n{x}_i^{k+1}+...+a_k\sum _{i=1}^n{x}_i^{2k}=\sum _{i=1}^n{x}_i^k{y}_i\end{array}\text{\hspace{1em}(4)}$$

将这些等式表达成矩阵的形式，可以得到如下矩阵:

$$\left[\begin{array}{cccc}n& \sum _{i=1}^n{x}_i& \cdots & \sum _{i=1}^n{x}_i^k\\ \sum _{i=1}^n{x}_i& \sum _{i=1}^n{x}_i^2& \cdots & \sum _{i=1}^n{x}_i^{k+1}\\ ⋮& ⋮& \ddots & ⋮\\ \sum _{i=1}^n{x}_k& \sum _{i=1}^n{x}_i^{k+1}& \cdots & \sum _{i=1}^n{x}_i^{2k}\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_k\end{array}\right]=\left[\begin{array}{c}\sum _{i=1}^n{y}_i\\ \sum _{i=1}^n{x}_i{y}_i\\ ⋮\\ \sum _{i=1}^n{x}_i^k{y}_i\end{array}\right]\text{\hspace{1em}(5)}$$

然后接着将这个范德蒙德行列式化简，可以得到：
$$\left[\begin{array}{cccc}1& x_1& \cdots & x_1^k\\ 1& x_2& \cdots & x_2^k\\ ⋮& ⋮& \ddots & ⋮\\ 1& x_n& \cdots & x_n^k\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_k\end{array}\right]=\left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_n\end{array}\right]\text{\hspace{1em}(6)}$$

接着我们的分析$XA=Y$，那么$A=X^{-1}Y$。由于这里需求$X$的逆矩阵，求逆矩阵的前提是满秩，但是这里可能没有办法满足这一点，所以我们利用广义矩阵可以得到，$A=(X^T\ast X{)}^{-1}\ast X^T\ast Y$,这样便得到了系数矩阵$A$，我们也就进一步得到了拟合曲线。

3.2、Matlab代码

clc;
clear;
%录入X轴数据
for a = 1:30 
    x(a) = a-1;
end
%录入Y轴数据
y=[1,2,3,6,6,7,8,9,8,10,9,19,16,14,15,24,28,36,40,40,42,41,37,39,52,52,56,57,62,69];

figure
plot(x,y,'.');%画点
hold on

k=5;%阶数  阶数可以在1-5之间更改看效果，记得每次更改了之后clear workspace然后在运行

%X数据录入
for a = 0:k
    for i = 1:30
        X(i,(a+1)) = x(i).^(a);
    end
end

Y = y';
A = (X'*X)^-1*X'*Y;%求矩阵系数A
A = A';%转置矩阵方便使用

z = 0:0.1:30;


if k==5
    y1 = A(1)+A(2).*z+A(3).*z.^2+A(4).*z.^3+A(5).*z.^4+A(6).*z.^5;%最后表达式用于绘图
elseif k==4
    y1 = A(1)+A(2).*z+A(3).*z.^2+A(4).*z.^3+A(5).*z.^4;%最后表达式用于绘图
elseif k==3
    y1 = A(1)+A(2).*z+A(3).*z.^2+A(4).*z.^3;%最后表达式用于绘图
elseif k==2
    y1 = A(1)+A(2).*z+A(3).*z.^2;%最后表达式用于绘图
elseif k==1
    y1 = A(1)+A(2).*z;%最后表达式用于绘图
end

plot(z,y1);
hold off

3.3、绘制效果

k=1

k=2

k=3

k=4

k=5