地下水动力学中用到的最小二乘法
欢迎点评！如果公式乱码，请用浏览器访问，用鼠标右键 -> 加载映像 显示公式。
需要安装的库：Jupyter，Jupyterlab，numpy，scipy，matplotlib，ipympl，mpl_interactions
用 pip 安装：
pip install Jupyter Jupyterlab numpy scipy matplotlib ipympl mpl_interactions
以下内容源自地下水动力学课程教学内容，程序可在 Jupyter notebook 中运行。
最小二乘法是一种用途非常广泛的算法，原理简单易懂，高等数学中有介绍。实际上，诸如多元线性回归等方法都是依此为基础。
地下水动力学中多个知识模块都可以用最小二乘法解决，如 Jacob 井损计算， 曲线类型确定，Jacob 线性拟合求参，Theis 公式非线性拟合求参等。
最小二乘法原理
根据观测数据 ，求拟合直线
，使
取最小值。
令 ，
写成矩阵形式，
式中，
其中，,
，
，
，
解得
一般形式
考虑超定方程组（超定指未知数小于方程个数）
其中， 代表样本数，
代表参数维度，写成矩阵形式
为得到 的最佳估计
，将问题转化为如下的最值问题：
通过微分求最值，得
若 为非奇异矩阵，则
有唯一解。
可以看出，求解最小二乘问题的关键是构造方程组。numpy 库中的 numpy.linalg.solve 可用于求解形如 的方程组，
numpy.linalg.lstsq 是解超定方程组 的最小二乘法程序。
scipy 库的 scipy.linalg.lstsq 也是最小二乘法。
算例
一次模拟实验中，输入 （自变量）为
(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)，
观测到的输出 （因变量）为
(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)。
根据实验分析， 与
成线性关系，试确定关系表达式。
以下为 Python 程序，其中 numpy.array.T将矩阵转置，matplotlib 库的 matplotlib.pyplot 模块用于绘制图形。
%matplotlib widget
import matplotlib.pyplot as plt
import numpy as np
# 导入 numpy 与 matplotlib 库
# 控制小数的显示精度
np.set_printoptions(precision=4)
# 将原始数据转化为 numpy 的 array 数组
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）计算：
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]
显示计算结果并绘制图形
# 显示计算结果与误差
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() # 获取当前图片, 用 fig.savefig("out.png") 保存
plt.show()
可调用 np.linalg.solve 求解
# 求解方程组，beta 返回解
beta = np.linalg.solve(A, b)
# 显示计算结果与误差
print("beta =", beta, "residuals =",
"{:.4f}".format(sum((beta[0] + beta[1]*t - y)**2)))
也可调用 np.linalg.lstsq 求解，参数可参考在线帮助文档。
# 形成方程组
X = np.vstack([np.ones(len(t)), t]).T
# 调用 np.linalg.lstsq
beta = np.linalg.lstsq(X, y, rcond=None)
# 显示计算结果与误差
print("beta =", beta[0], "residuals =", "{:.4f}".format(beta[1][0]))
