• 最小二乘法理論

最小二乘法（又稱最小平方法）是一種數學優化技術。它通過最小化誤差的平方和尋找數據的最佳函數匹配。利用最小二乘法可以簡便地求得未知的數據，并使得這些求得的數據與實際數據之間誤差的平方和為最小。是解決曲線拟合最常用的方法，其思路如下：

其中，是預選定的一組線性相關的函數，是待定系數，拟合準則是使與的距離的平方和最小，稱為最小二乘法準則。

最小二乘準則進行最小二乘平差計算的一個基本原則

• 代碼示例

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

結果如下：

最小二乘法很基本也很常用，其本質是插值，但是我發現在當數值非常大時它的效果卻不是很好，這時候就需要用到一些其他的方法，比如克裡金法等。

聲明：僅供參考