• 最小二乘法理论

最小二乘法（又称最小平方法）是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。是解决曲线拟合最常用的方法，其思路如下：

其中，是预选定的一组线性相关的函数，是待定系数，拟合准则是使与的距离的平方和最小，称为最小二乘法准则。

最小二乘准则进行最小二乘平差计算的一个基本原则

• 代码示例

import pandas as pd import numpy as np from sklearn.linear_model import LinearRegression import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D #---------------------------------------------------------------------------------------------------------------------- # Step1:创建需要被拟合的目标。三维空间中，定义4x4的网格，首先定义z值，而网格点上的z值不一样，我们所要做的就是根据这个z值去拟 # 合这个面上所有点的值。 #---------------------------------------------------------------------------------------------------------------------- np.random.seed(0) dim = 4 Z = (np.ones((dim, dim)) * np.arange(1, dim 1, 1))**3 np.random.rand(dim, dim) * 200 x = np.arange(1, dim 1).reshape(-1, 1) y = np.arange(1, dim 1).reshape(1, -1) X, Y = np.meshgrid(x, y) #---------------------------------------------------------------------------------------------------------------------- # Step2:自定义一组线性相关的函数, 3阶。 #---------------------------------------------------------------------------------------------------------------------- features = {} features['x^0*y^0'] = np.matmul(x**0, y**0).flatten() features['x*y'] = np.matmul(x, y).flatten() features['x*y^2'] = np.matmul(x, y**2).flatten() features['x^2*y^0'] = np.matmul(x**2, y**0).flatten() features['x^2*y'] = np.matmul(x**2, y).flatten() features['x^3*y^2'] = np.matmul(x**3, y**2).flatten() features['x^3*y'] = np.matmul(x**3, y).flatten() features['x^0*y^3'] = np.matmul(x**0, y**3).flatten() dataset = pd.DataFrame(features) #---------------------------------------------------------------------------------------------------------------------- # Step3:将选定函数与目标值带入SkLearn包中的线性回归拟合模块，它可以使平方和最小，结果返回截距和斜率。 #---------------------------------------------------------------------------------------------------------------------- reg = LinearRegression().fit(dataset.values, Z.flatten()) # reg.intercept_为截距, reg.coef_为斜率 z_pred = reg.intercept_ np.matmul(dataset.values, reg.coef_.reshape(-1, 1)).reshape(dim, dim) #---------------------------------------------------------------------------------------------------------------------- # Step4:可视化。 #---------------------------------------------------------------------------------------------------------------------- fig = plt.figure(figsize=(5, 5)) ax = Axes3D(fig) ax.plot_surface(X, Y, z_pred, label='prediction', cmap=plt.get_cmap('rainbow')) ax.scatter(X, Y, Z, c='r', label='datapoints') plt.show()

结果如下：

图1

上例定义的多项式阶数为3，对于大多数问题已经足够了，如果想定义更高阶数，则可参考如下代码：

import itertools import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D #----------------------------------------------------------------------------------------------------------------------- # Step1: 创建线性相关的函数，阶数自己定义，结果为拟合系数； #----------------------------------------------------------------------------------------------------------------------- def polyfit2d(x, y, z, order): ncols = (order 1)**2 G = np.zeros((x.size, ncols)) ij = itertools.product(range(order 1), range(order 1)) for k, (i, j) in enumerate(ij): G[:, k] = x**i * y**j m, _, _, _ = np.linalg.lstsq(G, z, rcond=-1) # lstsq的输出包括四部分：回归系数、残差平方和、自变量X的秩、X的奇异值 return m #----------------------------------------------------------------------------------------------------------------------- # Step2: 创建拟合函数，将欲拟合值和拟合系数待入，返回预测值； #----------------------------------------------------------------------------------------------------------------------- def polyval2d(x, y, m): order = int(np.sqrt(len(m))) - 1 # 根据多项式的列数反算阶数 ij = itertools.product(range(order 1), range(order 1)) z = np.zeros_like(x) for a, (i, j) in zip(m, ij): z = a * x**i * y**j return z #----------------------------------------------------------------------------------------------------------------------- # Step3: 示例； #----------------------------------------------------------------------------------------------------------------------- x = np.array([4, 5, 5, 4]) y = np.array([2, 3, 4, 5]) z = np.array([2, 3, 4, 7]) N_ORDER = 4 m = polyfit2d(x, y, z, N_ORDER) N_MESH = 10 xx, yy = np.meshgrid( np.linspace(x.min(), x.max(), N_MESH), np.linspace(y.min(), y.max(), N_MESH)) zz = polyval2d(xx, yy, m) #----------------------------------------------------------------------------------------------------------------------- # Step4: 可视化； #----------------------------------------------------------------------------------------------------------------------- fig = plt.figure(figsize=(5, 5)) ax = Axes3D(fig) ax.plot_surface(xx, yy, zz, label='prediction', cmap=plt.get_cmap('rainbow')) ax.scatter(x, y, z, c='r', label='datapoints') plt.show()

结果如下：

最小二乘法很基本也很常用，其本质是插值，但是我发现在当数值非常大时它的效果却不是很好，这时候就需要用到一些其他的方法，比如克里金法等。

声明：仅供参考