这里写目录标题

• 原理
• 代码实现

原理

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

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

对于y=a0+a1x+ey=a_{0} + a_{1}x + ey=a0​+a1​x+e，在一组观测数据中拟合最好的模型，即使得残差eee的平方和最小。分别对a0、a1a_{0}、a_{1}a0​、a1​求偏导，可得：
∑i=1na0+∑i=1nxia1=∑i=1nyi∑i=1nxia0+∑i=1nxi2a1=∑i=1nxiyi\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} i=1∑n​a0​+i=1∑n​xi​a1​=i=1∑n​yi​i=1∑n​xi​a0​+i=1∑n​xi2​a1​=i=1∑n​xi​yi​​
则：
a1=n∑i=1nxiyi−∑i=1nxi∑i=1nyin∑i=1nxi2−(∑i=1nxi)2a_{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}\\ a1​=n∑i=1n​xi2​−(∑i=1n​xi​)2n∑i=1n​xi​yi​−∑i=1n​xi​∑i=1n​yi​​
a0=yˉ−a1xˉa_{0} = \text{\={y}} - a_{1}\text{\={x}} a0​=yˉ​−a1​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_sizeaverage_y = sum_y / points_sizek = (points_size*sum_xy - sum_x*sum_y)/(points_size*sum_sqare_x - sum_x*sum_x)b = average_y - average_x*kbest_model = np.zeros((2, 1), dtype=np.float32)best_model[0, 0] = bbest_model[1, 0] = kreturn 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_htheta = np.dot(np.dot(np.linalg.inv(np.dot(x_mat.T, x_mat)), x_mat.T), y.T)return theta

参考链接1
参考链接2