最小二乘法

最小二乘法原理与编程实现

背景

数据与数据(变量与变量)之间,很多时候是存在一些关系的,如线性关系和非线性关系。我们常常会希望找到数据之间的关系,用一个函数(或者一条曲线)去描述两个变量之间的关系。
然而因为各种原因,测量得到的数据是会存在误差的,于是我们要一种方法去减少误差带来的干扰,尽可能的描绘出数据之间的关系,这个方法就是最小二乘法,也称最小方差法。

原理

比如说,我们得到了四对数据,分别是(1,0.9),(2, 2.3),(3.5, 2.6),(4, 3.8)把它们画在平面直角坐标系上

从图像上我们很容易发现,这些点大致分布在一条直线上面,于是我们大胆的猜测y与x呈线性关系,于是我们很自然地想要用一条直线去拟合他们。
也就是

\[y = kx + b \]

然而,事实是它们只是看起来像一条直线,但实际上并不是一条直线。把方程组列出来,

\[\begin{cases} 1k + b = 0.9 \\ 2k + b = 2.3 \\ 3.5k + b = 2.6 \\ 4k + b = 3.8 \\ \end{cases} \]

求解一下就会发现,这个方程组无解!从图像上来看就是我们找不到一条直线通过所有的点。
那该怎么办?那我们只好退一步,不要求找一条完全经过所有点的直线,只要求找一条能够大致刻画它们关系的直线,并且使得误差最小。那误差怎么衡量呢?
我们随便画一条直线,比如说

\[y = kx + b\\(图例k=1,b=0) \]

因为有些点没有落到直线上,于是我们把误差定义为每一个观测点的y值和我们预测的真实值之间的距离的平方,也就是

\[e_i = (f(x_i)-y_i)^2=(kx_i + b - y_i)^2 \]

我们的目标是使得总体的误差最小,也就是

\[min(\sum_{i=1}^n{e_i}) \]

这是一条关于k和b的二元函数,我们求偏导数并找到导数为0的点就可以使其最小,也即是令

\[\begin{cases}\frac{\partial{e}}{\partial{k}}=2\sum{(kx_i + b-y_i)x_i} = 0 \\\frac{\partial{e}}{\partial{b}}=2\sum{(kx_i + b-y_i)} = 0 \\\end{cases} \]

我们把原来的方程组写成矩阵的形式(这里把k,b当成待求参数,写成\(x_1,x_2\))

\[\begin{bmatrix} 1 & 1 \\ 2 & 1 \\ 3.5 & 1 \\ 4 & 1 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} =\begin{bmatrix} 0.6\\3\\2.6\\3.8 \end{bmatrix}\]

然后我们把误差写成列向量,也就是

\[\begin{bmatrix} e_1\\e_2\\e_3\\e_4 \end{bmatrix} = \begin{bmatrix} 1 & 1 \\ 2 & 1 \\ 3.5 & 1 \\ 4 & 1 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} -\begin{bmatrix} 0.6\\3\\2.6\\3.8 \end{bmatrix} \]

\[E =\begin{bmatrix}e_1\\e_2\\e_3\\e_4\end{bmatrix}, A =\begin{bmatrix}1 & 1 \\2 & 1 \\3.5 & 1 \\4 & 1 \end{bmatrix}, X = \begin{bmatrix}x_1 \\ x_2\end{bmatrix}, b = \begin{bmatrix}0.6\\3\\2.6\\3.8\end{bmatrix} \]

则得到

\[E = AX - b \]

要让\(E^2\)最小,即

\[E^TE = (AX-b)^T(AX-b) \]

我们对X求导

\[\frac{\partial{E^2}}{\partial{X}}= \frac{\partial{(AX-b)^T(AX-b)}}{\partial{X}}\\= \frac{\partial{(AX-b)^T}}{\partial{X}} *(AX-b)+ \frac{\partial{(AX-b)^T}}{\partial{X}}*(AX-b)\\= 2\frac{\partial{(AX-b)^T}}{\partial{X}}*(AX-b)\\ =2\frac{\partial{(AX)}^T}{\partial{X}}*(AX-b)\\ =2A^T(AX-b) \]

我们令其等于0,就可以解出X

\[X = (A^TA)^{-1}A^Tb \]

这是经过计算分析得到的,那有没有一些更加直觉化的解析呢?

有!我们可以通过向量的角度去理解它的原理!

刚刚那个矩阵方程可以写成向量方程的形式,

\[\begin{bmatrix}1 & 1 \\2 & 1 \\3.5 & 1 \\4 & 1 \end{bmatrix}\begin{bmatrix}x_1 \\ x_2\end{bmatrix}=x_1\begin{bmatrix}1 \\2 \\3.5 \\4 \\\end{bmatrix}+x_2\begin{bmatrix}1 \\1 \\1 \\1 \end{bmatrix}=\begin{bmatrix}0.6\\3\\2.6\\3.8\end{bmatrix} \]

我们记

\[v_1 = \begin{bmatrix} 1 \\ 2 \\ 3.5 \\ 4 \\ \end{bmatrix}, v_2 = \begin{bmatrix} 1 \\ 1 \\ 1 \\ 1 \end{bmatrix} b = \begin{bmatrix} 0.6\\3\\2.6\\3.8 \end{bmatrix} \]

容易看到,\(v_1\)\(v_2\)是线性无关的,也就是它们的线性组合可以张成四维空间里的一个平面。

我们想要找的是它的某种组合使得其等于b,但这种组合找不到,于是我们就去找那个平面里的一条线,使得它到b的距离最短。容易想象那条线就是b在span(\(v_1,v_2\))下的投影,我们记该投影向量为c ,则b和c之间的距离就是\(||b-c||\),我们记为e,这时e最小,且满足e垂直与span(\(v_1,v_2\))。所以e肯定也垂直于\(v_1,v_2\),即是

\[ev_1=0,ev_2=0\\ev_1+ev_2=0\\\begin{bmatrix}1\\2\\3.5\\4\end{bmatrix}*e+\begin{bmatrix}1\\1\\1\\1\end{bmatrix}*e=0\\ \]

这里是点乘,写成矩阵相乘就是

\[\begin{bmatrix}1&2&3.5&4\\1&1&1&1\end{bmatrix}e=0\\ \]

注意到这里是A的转置\(A^T\),再把e=(AX-b)代入得

\[A^T(AX-b)=0 \]

结果是和上面推导的一样的

以上的方法对于二次,三次拟合都是成立的,具体可以参考其他资料。

拟合的结果

一次拟合

###二次拟合
###三次拟合
##代码(Python)
# -*- coding: utf-8 -*-
"""
Created on Mon Feb 17 13:51:43 2020
最小二乘法
@author: urahyou
"""
import numpy as np
import matplotlib.pyplot as plt

x = np.array([3.0, 5, 6, 8, 10])
y = np.array([5.0, 2, 1, 2, 4])

p1 = plt.scatter(x,y,c='red')

def LSD(x, y, n):
    N = x.size  #获取方程组个数
    A = np.ones(N)
    
    for i in range(1, n+1):
        A = np.vstack((A,x**i))  #垂直拼接
    A = A.T  #转置回来
    #求解
    B = np.linalg.inv(A.T@A)@A.T
    #求出解系数
    sol = np.dot(B, y)
    return sol

def poly(x,sol):
    y =  np.zeros_like(x)  #每一个x,对应一个y
    n = sol.size
    for i in range(n):
        y += sol[i]*x**i
    return y


sol = LSD(x,y,2)
X = np.arange(0, 14, 0.1)
Y = poly(X,sol)
p2 = plt.plot(X,Y)

posted @ 2020-02-23 14:40  裏表異体  阅读(68)  评论(0编辑  收藏