张叶安的小站

本页面只读。您可以查看源文件，但不能更改它。如果您觉得这是系统错误，请联系管理员。
====== 第五章 曲线拟合的最小二乘法 ======

===== 5.1 引言 =====

在科学实验和统计分析中，常常需要根据一组离散数据点 $(x_i, y_i)$（$i = 0, 1, \ldots, n$）来确定变量 $x$ 和 $y$ 之间的函数关系 $y = f(x)$。由于数据通常存在测量误差，我们不要求曲线经过所有数据点，而是寻求一条能反映数据总体趋势的曲线，使得各数据点到曲线的偏差（在某种意义下）最小。这就是**曲线拟合**问题。

最常用的方法是**最小二乘法**（Least Squares Method），它使各点偏差的平方和最小。

===== 5.2 线性拟合 =====

==== 5.2.1 问题的提法 ====

给定数据 $(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$$

最小。

==== 5.2.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}$$

==== 5.2.3 解的表达式 ====

引入记号：
  - $\bar{x} = \frac{1}{n+1}\sum x_i$，$\bar{y} = \frac{1}{n+1}\sum y_i$
  - $l_{xx} = \sum (x_i - \bar{x})^2 = \sum x_i^2 - (n+1)\bar{x}^2$
  - $l_{xy} = \sum (x_i - \bar{x})(y_i - \bar{y}) = \sum x_i y_i - (n+1)\bar{x}\bar{y}$
  - $l_{yy} = \sum (y_i - \bar{y})^2$

则：
$$b = \frac{l_{xy}}{l_{xx}}, \quad a = \bar{y} - b\bar{x}$$

==== 5.2.4 拟合效果评价 =====

**残差平方和**：$Q = \sum (y_i - a - bx_i)^2$

**相关系数**：
$$r = \frac{l_{xy}}{\sqrt{l_{xx} l_{yy}}} \in [-1, 1]$$

  - $|r|$ 越接近1，线性相关性越强
  - $r > 0$ 表示正相关，$r < 0$ 表示负相关
  - $r = 0$ 表示无线性相关

**决定系数**：$R^2 = r^2$，表示拟合优度。

===== 5.3 多项式拟合 =====

==== 5.3.1 一般形式 ====

求 $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$$

最小。

==== 5.3.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$ 较大时病态严重。

==== 5.3.3 用正交多项式拟合 ====

为避免病态，可用正交多项式 $\{\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)}$$

===== 5.4 非线性拟合 =====

==== 5.4.1 可线性化的非线性模型 ====

许多非线性模型可通过变量替换化为线性形式：

| 原模型 | 变换 | 线性形式 |
|--------|------|----------|
| $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$ |

==== 5.4.2 直接非线性最小二乘 ====

对于一般形式 $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法。

===== 5.5 多元线性拟合 =====

==== 5.5.1 二元线性拟合 ====

求 $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}$$

==== 5.5.2 一般多元线性回归 ====

求 $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}$$

===== 5.6 加权最小二乘 =====

当各数据点的可靠性不同时，引入权因子 $\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.7 典型例题 =====

**例题 5.1** 给定数据：
| $x$ | 1 | 2 | 3 | 4 | 5 |
|-----|---|---|---|---|---|
| $y$ | 4 | 4.5 | 6 | 8 | 8.5 |

求最小二乘拟合直线。

**解**：$n+1 = 5$，计算：
  - $\sum x_i = 15$，$\sum y_i = 31$
  - $\sum x_i^2 = 55$，$\sum x_i y_i = 105.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 |

计算：
  - $\sum x_i = 15$，$\sum Y_i = 8.373$
  - $\sum x_i^2 = 55$，$\sum x_i Y_i = 29.038$

法方程：
$$\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}$。

===== 5.8 习题 =====

**基础练习**

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.