最小二乘法及其python实现

所谓最小二乘法，即通过对数据进行拟合，使得拟合值与样本值的方差最小。

线性拟合

假设样本为 { x n } = x 1 , x 2 . . . x n , { y n } = y 1 , y 2 . . . y n \{x_n\}={x_1,x_2...x_n},\{y_n\}={y_1,y_2...y_n} {xn​}=x1​,x2​...xn​,{yn​}=y1​,y2​...yn​，其拟合之后的方程为$y=ax+b$。则拟合值与样本值之差即为误差，误差的平方和可以衡量总误差：
$$J(a,b)=\sum _{i=1}^n(y_i-a{x}_i-b{)}^2$$
对于误差函数，当其导数为0时有极值，故对误差函数求偏导数并使之为0：

$$\begin{array}{rlr}\frac{\partial J}{\partial a}& =\sum _{i=1}^n-2x_i(y_i-a{x}_i-b)& =0\\ \frac{\partial J}{\partial b}& =\sum _{i=1}^n-2(y_i-a{x}_i-b)& =0\end{array}$$可得

$$\begin{array}{rl}\sum _{i=1}^n{x}_i(y_i-a{x}_i-b)& =0\\ \sum _{i=1}^n(y_i-a{x}_i-b)& =0\end{array}$$

约定记号$S(xy)=\sum x_i{y}_i$, $S(x^2)=\sum x_i^2$, $S(x)=\sum x_i$, $E(x)=\frac{1}{n}\sum x_i$, $E(y)=\frac{1}{n}\sum y_i$,则上式化为

$$\begin{array}{rl}S(xy)-aS(x^2)-bS(x)& =0\\ E(y)-aE(x)-b& =0\end{array}$$
易得
$$\begin{array}{rl}a& =\frac{S(xy)-E(y)S(x)}{S(x^2)-E(x)S(x)}\\ b& =\frac{E(y)S(x^2)-S(xy)E(x)}{S(x^2)-E(x)S(x)}\end{array}$$
这个表达式还是非常简单的。

对于有些情况，我们往往选取自然序列作为自变量，这个时候在求自变量的取值时可以用到一些初等数学的推论，对于$x\in [m,n]$的自然序列来说，有
$$\begin{array}{rl}S(x)& =\frac{(m+n)(n-m+1)}{2}\\ S(x^2)& =\frac{n(n+1)(2n+1)-m(n-1)(2m-1)}{6}\end{array}$$代码为：

#文件名core.py
import numpy as np
def leastSquare(x,y):
    if len(x)==2:
    #此时x为自然序列
        sx = 0.5*(x[1]-x[0]+1)*(x[1]+x[0])
        ex = sx/(x[1]-x[0]+1)
        sx2 = ((x[1]*(x[1]+1)*(2*x[1]+1))
              -(x[0]*(x[0]-1)*(2*x[0]-1)))/6
        x = np.array(range(x[0],x[1]+1))
    else:
        sx = sum(x)
        ex = sx/len(x)
        sx2 = sum(x**2)
    
    sxy = sum(x*y)
    ey = np.mean(y)

    a = (sxy-ey*sx)/(sx2-ex*sx)
    b = (ey*sx2-sxy*ex)/(sx2-ex*sx)
    return a,b

测试一下

>>> x = np.arange(25)
>>> y = x*15+20+np.random.randn(len(x))*5	#randn生成正态分布噪声
>>> a,b = core.leastSquare(x,y)				
>>> plt.scatter(x,y)						#原始数据散点图
<matplotlib.collections.PathCollection object at 0x00000218DEBBEDC8>
>>> plt.plot(x,a*x+b)						#拟合直线
[<matplotlib.lines.Line2D object at 0x00000218E0314FC8>]
>>> plt.show()

得到

高阶多项式

对于高阶的多项式拟合，其思路与线性拟合是如出一辙的。对于样本 { x n } = x 1 , x 2 . . . x n , { y n } = y 1 , y 2 . . . y n \{x_n\}={x_1,x_2...x_n},\{y_n\}={y_1,y_2...y_n} {xn​}=x1​,x2​...xn​,{yn​}=y1​,y2​...yn​，假设其拟合之后的方程为$$y=\sum _{j=0}^m{a}_j{x}^j$$则相应地其误差方程组可表示为

$$J(a_j)=\sum _{i=0}^n(y_i-\sum _{j=0}^m{a}_j{x}_i^j{)}^2$$则其每个参数的偏导数可表示为

$$\frac{\partial J}{\partial a_k}=\sum _{i=0}^{n}2\cdot x_i^k(y_i-\sum _{j=0}^m{a}_j{x}_i^j)=0$$即

$$\sum _{i=0}^n{x}_i^k{y}_i-\sum _{i=0}^n\sum _{j=0}^m{a}_j{x}_i^{j+k}=0$$

和前面一样，约定

$$S_k=\sum _{i=0}^n{x}_i^k{y}_i,S_{kj}=\sum _{j=0}^m{x}_i^{j+k}$$，则对于任意$i$值，上式可变为
$$S_k-\cdot \sum _{j=0}^m{a}_{j}Skj=0$$

写成矩阵的形式即为
$$\left[\begin{array}{cccc}{S}_{00}& S_{01}& ...& S_{0m}\\ S_{10}& S_{11}& ...& S_{1m}\\ ...& ...& ...& ...\\ S_{m0}& S_{m1}& ...& S_{mm}\end{array}\right]\cdot \left[\begin{array}{c}{a}_0\\ a_1\\ ...\\ a_m\end{array}\right]=\left[\begin{array}{c}{S}_0\\ S_1\\ ...\\ S_m\end{array}\right]$$
代码如下

#传入参数格式为np.array,n为阶数
def leastSquareMulti(x,y,n):
    X = [np.sum(x**i) for i in range(2*n+1)]
    Y = np.array([[np.sum(y*x**i)] for i in range(n+1)])
    S = np.array([X[i:i+n+1] for i in range(n+1)])
    return np.linalg.solve(S,Y)		#

经测试结果如下：

>>> x = np.arange(25)
>>> y = x**3+3*x**2+2*x+12
>>> import core
>>> core.leastSquareMulti(x,y,3)
array([[12.],		#此为常数项
       [ 2.],
       [ 3.],
       [ 1.]])

多自变量

对于样本 { x 1 n } = x 11 , x 12 . . . x 1 n , { x 2 n } = x 21 , x 22 . . . x 2 n , . . . . . . . { x m n } = x m 1 , x m 2 . . . x m n , { y n } = y 1 , y 2 . . . y n \begin{aligned} \{x_{1n}\}&=x_{11},x_{12}&...&x_{1n},\\ \{x_{2n}\}&=x_{21},x_{22}&...&x_{2n},\\ &....&...&\\ \{x_{mn}\}&=x_{m1},x_{m2}&...&x_{mn},\\ \{y_n\}&=y_1,y_2&...&y_n \end{aligned} {x1n​}{x2n​}{xmn​}{yn​}​=x11​,x12​=x21​,x22​....=xm1​,xm2​=y1​,y2​​...............​x1n​,x2n​,xmn​,yn​​
假设其拟合之后的方程为$$y=\sum _{j=1}^m{a}_j{x}_j$$则相应地其误差方程组可表示为

$$J(a_j)=\sum _{i=1}^n(y_i-\sum _{j=1}^m{a}_j{x}_{ji}{)}^2$$则其每个参数的偏导数可表示为

$$\frac{\partial J}{\partial a_k}=\sum _{i=1}^{n}2\cdot x_{ki}(y_i-\sum _{j=1}^m{a}_j{x}_{ji})=0$$即

$$\sum _{i=1}^n{x}_{ki}{y}_i-\sum _{i=1}^n\sum _{j=1}^m{a}_j{x}_{ji}{x}_{ki}=0$$

约定$Y_k={\sum }_{i=1}^n{x}_{ki}{y}_i$, $X_{jk}={\sum }_{i=1}^n{x}_{ji}{x}_{ki}$,其矩阵形式为
$$\left[\begin{array}{cccc}{X}_{11}& X_{12}& ...& X_{1m}\\ X_{21}& X_{22}& ...& X_{2m}\\ ...& ...& ...& ...\\ X_{m1}& X_{m2}& ...& X_{mm}\end{array}\right]\cdot \left[\begin{array}{c}{a}_1\\ a_2\\ ...\\ a_m\end{array}\right]=\left[\begin{array}{c}{Y}_1\\ Y_2\\ ...\\ Y_m\end{array}\right]$$

如果最终的拟合方程需要常数项，则只需对$x$增添一组值为1的样本即可，其对应的$a_{m+1}$即为常数项。

在具体的编程中，假设其输入的自变量为一个矩阵$X$，每行代表某一自变量的不同取值，列表示每一组取值的不同自变量。那么上式左侧的系数矩阵可以表示为$X\cdot X^T$。

指数函数

一般来说，对于形如$y=a{e}^{bx}$这样的函数来说，只需左右取对数，便可得到形如$lny=bx+lna$这样的线性形式，通过简单的坐标变换，即可得到$b$和$lna$的值。

然而，对于形如$y=a_1{e}^{b_{1}x}+a_2{e}^{b_{2}x}$的函数，便无能为力了。

这时，如果$x$是一个自然序列，或者间距恒定，那么我们可以通过上述表达式构建一个线性关系。设$x$的间距为$\delta$，约定$$\begin{array}{rl}{y}_0& =a_1{e}^{b_{1}x}+a_2{e}^{b_{2}x}=y\\ y_1& =a_1{e}^{b_1(x+\delta )}+a_2{e}^{b_2(x+\delta )}\\ y_2& =a_1{e}^{b_1(x+2\delta )}+a_2{e}^{b_2(x+2\delta )}\end{array}$$

对于上式，可以得到关系$y_2=y_1\cdot (e^{b_1\delta }+e^{b_2\delta })-y_0\cdot e^{b_1\delta }\cdot e^{b_2\delta }$，即可通过最小二乘法求出$(e^{b_1\delta }+e^{b_2\delta })$和$e^{b_1\delta }\cdot e^{b_2\delta }$。

则由$y_2/y_1,y_2/y_0$得到的$a,b$可组成一元二次方程$x^2-a\cdot x-b=0$的解即为$e^{b_1\delta }$和$e^{b_2\delta }$。由于$\delta$是定值，故可得到$b_1、b_2$的值。

这时，我们可以选取两种不同的技术方案，其一是将其转化为多元最小二乘拟合，令$x_1=e^{b_{1}x},x_2=e^{b_{2}x}$。另一种则是令$X=e^{(b_1-b_2)x},Y=\frac{y}{e^{b_{2}x}}$，则拟合方程变为$Y=a_{1}X+a_2$。

令$\delta$为1，则其代码为

def expFit(x,y):
    y0 = y[0:-3]
    y1 = y[1:-2]
    y2 = y[2:-1]

    B,C = leastSquare(y2/y0,y1/y0)
    b1 = np.log((B-np.sqrt(B**2+4*C))/2)
    b2 = np.log((B+np.sqrt(B**2+4*C))/2)

    X = np.exp(b1-b2)*x
    Y = y/np.exp(b2*x)

    a1,a2 = leastSquare(X,Y)
    return a1,a2,b1,b2