在科学实验和统计分析中，常常需要根据一组离散数据点 $(x_i, y_i)$（$i = 0, 1, \ldots, n$）来确定变量 $x$ 和 $y$ 之间的函数关系 $y = f(x)$。由于数据通常存在测量误差，我们不要求曲线经过所有数据点，而是寻求一条能反映数据总体趋势的曲线，使得各数据点到曲线的偏差（在某种意义下）最小。这就是曲线拟合问题。
最常用的方法是最小二乘法（Least Squares Method），它使各点偏差的平方和最小。
给定数据 $(x_i, y_i)$（$i = 0, 1, \ldots, n$），求一次函数： $$y = a + bx$$
使得残差平方和： $$Q(a, b) = \sum_{i=0}^{n} (y_i - a - bx_i)^2$$
最小。
由极值必要条件： $$\frac{\partial Q}{\partial a} = -2\sum_{i=0}^{n} (y_i - a - bx_i) = 0$$ $$\frac{\partial Q}{\partial b} = -2\sum_{i=0}^{n} x_i(y_i - a - bx_i) = 0$$
整理得法方程组（正规方程组）：
$$\begin{cases}
(n+1)a + b\sum x_i = \sum y_i
a\sum x_i + b\sum x_i^2 = \sum x_i y_i
\end{cases}$$
写成矩阵形式：
$$\begin{bmatrix} n+1 & \sum x_i
\sum x_i & \sum x_i^2 \end{bmatrix} \begin{bmatrix} a
b \end{bmatrix} = \begin{bmatrix} \sum y_i
\sum x_i y_i \end{bmatrix}$$
引入记号：
则： $$b = \frac{l_{xy}}{l_{xx}}, \quad a = \bar{y} - b\bar{x}$$
残差平方和：$Q = \sum (y_i - a - bx_i)^2$
相关系数： $$r = \frac{l_{xy}}{\sqrt{l_{xx} l_{yy}}} \in [-1, 1]$$
决定系数：$R^2 = r^2$，表示拟合优度。
求 $m$ 次多项式（$m < n$）： $$P_m(x) = a_0 + a_1 x + a_2 x^2 + \cdots + a_m x^m$$
使得： $$Q(a_0, a_1, \ldots, a_m) = \sum_{i=0}^{n} \left(y_i - \sum_{j=0}^{m} a_j x_i^j\right)^2$$
最小。
$$\sum_{j=0}^{m} \left(\sum_{i=0}^{n} x_i^{k+j}\right) a_j = \sum_{i=0}^{n} x_i^k y_i, \quad k = 0, 1, \ldots, m$$
矩阵形式：
$$\begin{bmatrix}
s_0 & s_1 & \cdots & s_m
s_1 & s_2 & \cdots & s_{m+1}
\vdots & \vdots & \ddots & \vdots
s_m & s_{m+1} & \cdots & s_{2m}
\end{bmatrix}
\begin{bmatrix} a_0
a_1
\vdots
a_m \end{bmatrix} =
\begin{bmatrix} t_0
t_1
\vdots
t_m \end{bmatrix}$$
其中 $s_k = \sum_{i=0}^{n} x_i^k$，$t_k = \sum_{i=0}^{n} x_i^k y_i$。
系数矩阵是Hankel矩阵，当 $m$ 较大时病态严重。
为避免病态，可用正交多项式 $\{\varphi_j(x)\}$： $$P_m(x) = c_0 \varphi_0(x) + c_1 \varphi_1(x) + \cdots + c_m \varphi_m(x)$$
其中 $\varphi_j(x)$ 满足离散正交性： $$\sum_{i=0}^{n} \varphi_k(x_i) \varphi_j(x_i) = 0 \quad (k \neq j)$$
则： $$c_j = \frac{\sum_{i=0}^{n} y_i \varphi_j(x_i)}{\sum_{i=0}^{n} \varphi_j^2(x_i)}$$
许多非线性模型可通过变量替换化为线性形式：
| 原模型 | 变换 | 线性形式 | |
| ——– | —— | ———- | |
| $y = ae | {bx}$ | $Y = \ln y$ | $Y = \ln a + bx$ |
|---|---|---|---|
| $y = ax | b$ | $Y = \ln y$, $X = \ln x$ | $Y = \ln a + bX$ |
| $y = \frac{x}{a + bx}$ | $Y = \frac{1}{y}$, $X = \frac{1}{x}$ | $Y = aX + b$ | |
| $y = \frac{a}{b + x}$ | $Y = \frac{1}{y}$ | $Y = \frac{b}{a} + \frac{1}{a}x$ | |
| $y = ae | {b/x}$ | $Y = \ln y$, $X = \frac{1}{x}$ | $Y = \ln a + bX$ |
对于一般形式 $y = f(x; \alpha_1, \ldots, \alpha_m)$，最小化： $$Q(\alpha_1, \ldots, \alpha_m) = \sum_{i=0}^{n} [y_i - f(x_i; \alpha_1, \ldots, \alpha_m)]^2$$
这需要求解非线性方程组（由 $\frac{\partial Q}{\partial \alpha_k} = 0$ 得到），通常用迭代法如Gauss-Newton法或Levenberg-Marquardt法。
求 $z = a + bx + cy$，最小化： $$Q(a, b, c) = \sum_{i=0}^{n} (z_i - a - bx_i - cy_i)^2$$
法方程：
$$\begin{cases}
(n+1)a + b\sum x_i + c\sum y_i = \sum z_i
a\sum x_i + b\sum x_i^2 + c\sum x_i y_i = \sum x_i z_i
a\sum y_i + b\sum x_i y_i + c\sum y_i^2 = \sum y_i z_i
\end{cases}$$
求 $y = \beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p$。
矩阵形式： $$\mathbf{y} = X\boldsymbol{\beta} + \boldsymbol{\varepsilon}$$
其中 $X$ 是 $n \times (p+1)$ 的设计矩阵。最小二乘解： $$\hat{\boldsymbol{\beta}} = (X^T X)^{-1} X^T \mathbf{y}$$
当各数据点的可靠性不同时，引入权因子 $\omega_i > 0$： $$Q = \sum_{i=0}^{n} \omega_i (y_i - f(x_i))^2$$
权越大表示该点越重要。通常取 $\omega_i = \frac{1}{\sigma_i^2}$，其中 $\sigma_i^2$ 是测量方差。
例题 5.1 给定数据：
| $x$ | 1 | 2 | 3 | 4 | 5 |
| —– | — | — | — | — | — |
| $y$ | 4 | 4.5 | 6 | 8 | 8.5 |
求最小二乘拟合直线。
解：$n+1 = 5$，计算：
法方程：
$$\begin{cases}
5a + 15b = 31
15a + 55b = 105.5
\end{cases}$$
解得：$a = 2.45$，$b = 1.25$。
拟合直线：$y = 2.45 + 1.25x$。
例题 5.2 对以下数据用 $y = ae^{bx}$ 拟合：
| $x$ | 1 | 2 | 3 | 4 | 5 |
| —– | — | — | — | — | — |
| $y$ | 2.3 | 3.5 | 5.4 | 8.1 | 12.3 |
解：取对数：$\ln y = \ln a + bx$。
令 $Y = \ln y$，求 $Y = A + bx$ 的线性拟合。
| $x$ | 1 | 2 | 3 | 4 | 5 |
| —– | — | — | — | — | — |
| $Y = \ln y$ | 0.833 | 1.253 | 1.686 | 2.092 | 2.509 |
计算：
法方程：
$$\begin{cases}
5A + 15b = 8.373
15A + 55b = 29.038
\end{cases}$$
解得：$A = 0.439$，$b = 0.412$。
因此 $a = e^{0.439} = 1.551$，拟合曲线：$y = 1.551e^{0.412x}$。
基础练习
1. 证明线性拟合的法方程有唯一解当且仅当 $x_i$ 不全相同。
2. 给定数据 $(1,2), (2,3), (3,5), (4,6)$，求最小二乘拟合直线。
3. 对以下数据求二次多项式拟合：
| $x$ | -2 | -1 | 0 | 1 | 2 | |-----|----|----|---|---|---| | $y$ | 5 | 1 | 0 | 2 | 4 |
进阶练习
4. 证明对于线性拟合，拟合直线必过点 $(\bar{x}, \bar{y})$。
5. 设 $H$ 是多项式拟合的Hankel矩阵，证明当 $x_i$ 互异时 $H$ 正定。
6. 证明多元线性回归的最小二乘解 $\hat{\beta} = (X^T X)^{-1} X^T y$ 满足法方程。
编程实践
7. 编写Python程序实现线性拟合，并与numpy.polyfit结果比较。
8. 用Python的scipy.optimize.curve_fit实现非线性最小二乘拟合。
1. 李庆扬等. 数值分析（第5版）. 清华大学出版社, 2008. 2. 何晓群. 应用回归分析. 中国人民大学出版社, 2019.