最小二乘法(ordinary least squares)趋势面拟合

• 最小二乘法理论

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

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

• 代码示例

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()

结果如下：

上例定义的多项式阶数为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()

结果如下：

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

声明：仅供参考