最小二乘多项式拟合曲线scipy.linalg.lstsq

python 的相关类库已经帮我们实现了具体算法，而我们只要学会使用就可以了。随着对技术的逐渐掌握及积累，当类库中的算法已经无法满足自身需求的时候，我们也可以尝试通过自己的方式实现各种算法

言归正传，什么是”最小二乘法”呢？

定义：最小二乘法（又称最小平方法）是一种数学优化技术，它通过最小化误差的平方和寻找数据的最佳函数匹配。

作用：利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。

原则：以 ”残差平方和最小” 确定直线位置(在数理统计中，残差是指实际观察值与估计值之间的差)。

基本思路：对于一元线性回归模型, 假设从总体中获取了 $n$ 组观察值（X1，Y1），（X2，Y2）， …，（Xn，Yn），对于平面中的这 $n$ 个点，可以使用无数条曲线来拟合。而线性回归就是要求样本回归函数尽可能好地拟合这组值，也就是说，这条直线应该尽可能的处于样本数据的中心位置。因此，选择最佳拟合曲线的标准可以确定为：使总的拟合误差（即总残差）达到最小。

案例

假定模型为 $y=2x+1$，我们一共有100组数据$(x_i,y_i)$，我们需要在这100组数据的基础上找出 $x_i$ 和 $y_i$ 的线性关系方程 $y=2x+1$。以下是代码实现：

import scipy.linalg as la
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
import warnings
warnings.filterwarnings('ignore')
matplotlib.use('TkAgg')
m = 100
x = np.linspace(-1, 1, m)
y_exact = 1 + 2 * x
xi = x + np.random.normal(0, 0.05, 100)
yi = 1 + 2 * xi + np.random.normal(0, 0.05, 100)
print(xi, "#xi")
print(yi, "#yi")
A = np.vstack([xi**0, xi**1])
sol, r, rank, s = la.lstsq(A.T, yi)
y_fit = sol[0] + sol[1] * x
print(sol, r, rank, s)
fig, ax = plt.subplots(figsize=(12, 4))
ax.plot(xi, yi, 'go', alpha=0.5, label='Simulated data')
ax.plot(x, y_exact, 'k', lw=2, label='True value $y = 1 + 2x$')
ax.plot(x, y_fit, 'b', lw=2, label='Least square fit')
ax.set_xlabel(r"$x$", fontsize=10)
ax.set_ylabel(r"$y$", fontsize=10)
ax.legend(loc=2)
plt.show()

输出

In[2]: runfile('E:/09-code/06-Turbulent_flow_spectrum/test/test1.py', wdir='E:/09-code/06-Turbulent_flow_spectrum/test')
[-0.93994976 -0.91969219 -0.83232459 -0.99153433 -0.93666776 -0.88241546
 -1.02403388 -0.89335338 -0.80696476 -0.79658244 -0.80776571 -0.77688263
								……
  0.67733455  0.73973969  0.83918006  0.83862599  0.82169548  0.8497527
  0.85869333  0.75000611  0.86614204  0.92745979  0.87372474  0.84313477
  0.9961282   0.93191607  0.9788632   0.90734897] #xi
  
[-0.92398102 -0.85616797 -0.62852475 -1.02604021 -0.95541186 -0.752893
 -1.07587482 -0.81073224 -0.61987295 -0.54452371 -0.58916095 -0.60001989
								……
  2.74031756  2.56113636  2.81197947  2.90151446  2.80454031  2.71162451
  3.07521601  2.94532927  3.00815703  2.87686634] #yi
  
[0.99643923 2.01948052] 0.2597864247412017 2 [10.00000063  5.8821223 ]

参考资料
[1] numpy np.polyfit()（最小二乘多项式拟合曲线）（有待进一步研究）；
[2] 详解Python中使用最小二乘法方法 2017.3；
[3] 机器学习（三）—python实现最小二乘法 2021.9
[4] 一文让你彻底搞懂最小二乘法（超详细推导） 2022.5；
[5] python最小二乘法一个例子（初级）2019.8；