原理

最小二乘法适用于对处理的一堆数据，不必精确的经过每一点，而是根据图像到每个数据点的距离和最小确定函数。需要注意的是，最小二乘是对全局进行拟合优化，对噪声比较敏感，所以如果有噪声比较大的观测值会影响拟合结果，此时建议用RANSAC算法拟合。

最小二乘法逼近的最简单的例子是根据一组观测值对(x1,y1),(x2,y2)…(xn,yn)来拟合一条直线。
在这里插入图片描述
对于 $y=a_{0} + a_{1}x + e$ ，在一组观测数据中拟合最好的模型，即使得残差 $e$ 的平方和最小。分别对 $a_{0}、a_{1}$ 求偏导，可得：
$\begin {matrix} \displaystyle\sum_{i=1}^na_{0} + \displaystyle\sum_{i=1}^nx_{i}a_{1} = \displaystyle\sum_{i=1}^ny_{i} \\ \displaystyle\sum_{i=1}^nx_{i}a_{0} + \displaystyle\sum_{i=1}^nx_{i}^2a_{1} = \displaystyle\sum_{i=1}^nx_{i}y_{i} \end {matrix}$
则：
$a_{1} = \cfrac {n\sum_{i=1}^nx_{i}y_{i} - \sum_{i=1}^n x_{i}\sum_{i=1}^n y_{i}}{n \sum_{i=1}^n x_{i}^2 - (\sum_{i=1}^n x_{i})^2}\\$
$a_{0} = \text{\={y}} - a_{1}\text{\={x}}$
对n次多项式同样适用。

代码实现

def leastsq_fit(points:np.ndarray):
    points_size = len(points)
    points_sum = np.sum(points, axis=0)
    sum_xy = np.sum(points[:, 0] * points[:, 1], axis=0)
    sum_x = points_sum[0]
    sum_y = points_sum[1]
    sum_sqare_x = np.sum(points ** 2, axis=0)[0]
    average_x = sum_x / points_size
    average_y = sum_y / points_size
    k = (points_size*sum_xy - sum_x*sum_y)/(points_size*sum_sqare_x - sum_x*sum_x)
    b = average_y - average_x*k
    best_model = np.zeros((2, 1), dtype=np.float32)
    best_model[0, 0] = b
    best_model[1, 0] = k

    return best_model

矩阵形式参考我另外的博客链接1，链接2

def leastsq_mutifunc(x, y, m):
    """
    多项式最小二乘法实现, 矩阵形式
    :param x:输入
    :param y:目标输出
    :param m:多项式阶数
    :return:多项式系数
    """
    x = np.array(x)
    y = np.array(y)

    assert m <= x.shape[0], f"the number of m({m}) need less than x's size({x.shape[0]})"
    assert x.shape[0] == y.shape[0], f"the size of x({x.shape[0]}) must equal to y's size({y.shape[0]}"
    x_mat = np.zeros((x.shape[0], m+1))
    for i in range(x.shape[0]):
        x_mat_h = np.zeros((1, m+1))
        for j in range(m+1):
            x_mat_h[0][j] = x[i] ** (m-j)
        x_mat[i] = x_mat_h
    theta = np.dot(np.dot(np.linalg.inv(np.dot(x_mat.T, x_mat)), x_mat.T), y.T)
    return theta