python:最小二乘法拟合原理及代码实现

原理

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

最小二乘法逼近的最简单的例子是根据一组观测值对(x1,y1),(x2,y2)…(xn,yn)来拟合一条直线。

对于$y=a_0+a_{1}x+e$，在一组观测数据中拟合最好的模型，即使得残差$e$的平方和最小。分别对$a_0、a_1$求偏导，可得：
$$\begin{array}{c}\sum _{i=1}^n{a}_0+\sum _{i=1}^n{x}_i{a}_1=\sum _{i=1}^n{y}_i\\ \sum _{i=1}^n{x}_i{a}_0+\sum _{i=1}^n{x}_i^2{a}_1=\sum _{i=1}^n{x}_i{y}_i\end{array}$$
则：
$$a_1=\frac{n{\sum }_{i=1}^n{x}_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
3
4
5
6
7
8
9
10
11
12
13
14
15
16

矩阵形式参考我另外的博客链接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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

参考链接1
参考链接2
参考链接3