地下水動(dòng)力學(xué)中用到的最小二乘法
歡迎點(diǎn)評(píng)融撞！如果公式亂碼，請(qǐng)用瀏覽器訪問(wèn)粗蔚，用鼠標(biāo)右鍵 -> 加載映像 顯示公式尝偎。
需要安裝的庫(kù)：Jupyter，Jupyterlab鹏控，numpy致扯，scipy，matplotlib当辐，ipympl抖僵，mpl_interactions
用 pip 安裝：
pip install Jupyter Jupyterlab numpy scipy matplotlib ipympl mpl_interactions
以下內(nèi)容源自地下水動(dòng)力學(xué)課程教學(xué)內(nèi)容，程序可在 Jupyter notebook 中運(yùn)行缘揪。
最小二乘法是一種用途非常廣泛的算法耍群，原理簡(jiǎn)單易懂，高等數(shù)學(xué)中有介紹找筝。實(shí)際上蹈垢，諸如多元線性回歸等方法都是依此為基礎(chǔ)。
地下水動(dòng)力學(xué)中多個(gè)知識(shí)模塊都可以用最小二乘法解決袖裕，如 Jacob 井損計(jì)算曹抬， 曲線類(lèi)型確定，Jacob 線性擬合求參急鳄，Theis 公式非線性擬合求參等谤民。
最小二乘法原理
根據(jù)觀測(cè)數(shù)據(jù) ，求擬合直線
疾宏，使
取最小值张足。
令 ，
寫(xiě)成矩陣形式灾锯，
式中兢榨，
其中，,
顺饮，
，
凌那，
解得
一般形式
考慮超定方程組（超定指未知數(shù)小于方程個(gè)數(shù)）
其中兼雄， 代表樣本數(shù)，
代表參數(shù)維度帽蝶，寫(xiě)成矩陣形式
為得到 的最佳估計(jì)
赦肋，將問(wèn)題轉(zhuǎn)化為如下的最值問(wèn)題：
通過(guò)微分求最值，得
若 為非奇異矩陣，則
有唯一解佃乘。
可以看出囱井，求解最小二乘問(wèn)題的關(guān)鍵是構(gòu)造方程組。numpy 庫(kù)中的 numpy.linalg.solve 可用于求解形如 的方程組趣避，
numpy.linalg.lstsq 是解超定方程組 的最小二乘法程序庞呕。
scipy 庫(kù)的 scipy.linalg.lstsq 也是最小二乘法。
算例
一次模擬實(shí)驗(yàn)中程帕，輸入 （自變量）為
(0.00, 0.30, 0.60, 0.90, 1.08, 1.20, 1.30, 1.48, 1.60, 1.70, 1.78, 1.85, 1.90, 1.95, 2.00)住练，
觀測(cè)到的輸出 （因變量）為
(1.78, 1.91, 2.01, 2.12, 2.20, 2.22, 2.25, 2.32, 2.38, 2.41, 2.43, 2.47, 2.49, 2.48, 2.51)。
根據(jù)實(shí)驗(yàn)分析愁拭， 與
成線性關(guān)系讲逛，試確定關(guān)系表達(dá)式。
以下為 Python 程序岭埠，其中 numpy.array.T將矩陣轉(zhuǎn)置盏混，matplotlib 庫(kù)的 matplotlib.pyplot 模塊用于繪制圖形。
%matplotlib widget
import matplotlib.pyplot as plt
import numpy as np
# 導(dǎo)入 numpy 與 matplotlib 庫(kù)
# 控制小數(shù)的顯示精度
np.set_printoptions(precision=4)
# 將原始數(shù)據(jù)轉(zhuǎn)化為 numpy 的 array 數(shù)組
t = np.array([0.00, 0.30, 0.60, 0.90, 1.08, 1.20, 1.30,
1.48, 1.60, 1.70, 1.78, 1.85, 1.90, 1.95, 2.00])
y = np.array([1.78, 1.91, 2.01, 2.12, 2.20, 2.22, 2.25,
2.32, 2.38, 2.41, 2.43, 2.47, 2.49, 2.48, 2.51])
# 形成方程組
X = np.vstack([np.ones(len(t)), t]).T
A = np.dot(X.T, X)
b = np.dot(X.T, y)
按公式（1）計(jì)算：
beta = np.zeros(2)
# 求解方程組
beta[1] = (A[0,0]*b[1] - A[0,1]*b[0]) / (A[0,0]*A[1,1] - A[0,1]*A[0,1])
beta[0] = (b[0] - A[0,1]*beta[1]) / A[0,0]
顯示計(jì)算結(jié)果并繪制圖形
# 顯示計(jì)算結(jié)果與誤差
print("beta =", beta, "residuals =",
"{:.4f}".format(sum((beta[0] + beta[1]*t - y)**2)))
# 繪制圖形
plt.plot(t, y, "*")
plt.plot(t, beta[0]+beta[1]*t, "-", label="最小二乘法")
plt.xlabel(r"$t$")
plt.ylabel(r"$y$")
plt.legend(
prop={'family': 'Simsun', 'size': 10}, handlelength=2,
loc=4, title="圖例",
title_fontproperties={'family': 'SimHei', 'size': 12})
plt.grid(True)
plt.tight_layout()
fig = plt.gcf() # 獲取當(dāng)前圖片, 用 fig.savefig("out.png") 保存
plt.show()
可調(diào)用 np.linalg.solve 求解
# 求解方程組惜论，beta 返回解
beta = np.linalg.solve(A, b)
# 顯示計(jì)算結(jié)果與誤差
print("beta =", beta, "residuals =",
"{:.4f}".format(sum((beta[0] + beta[1]*t - y)**2)))
也可調(diào)用 np.linalg.lstsq 求解括饶，參數(shù)可參考在線幫助文檔。
# 形成方程組
X = np.vstack([np.ones(len(t)), t]).T
# 調(diào)用 np.linalg.lstsq
beta = np.linalg.lstsq(X, y, rcond=None)
# 顯示計(jì)算結(jié)果與誤差
print("beta =", beta[0], "residuals =", "{:.4f}".format(beta[1][0]))
