非线性回归

一句话定位：当数据呈指数增长、饱和衰减、S 型爬升等曲线形态，且模型形式来自机理或先验知识时，用非线性最小二乘法估计曲线参数——这是数学建模竞赛中"机理建模 + 数据拟合"的标配工具。本文配套程序生成的全部图表与指标，均可直接复现（固定随机种子）。

一、算法含义

1.1 通俗理解：从"拟合直线"到"拟合曲线"

线性回归的模型是一条直线（或超平面）：

$y = \beta_0 + \beta_1 x + \varepsilon$

它假设 $y$ 与 $x$ 呈直线关系。但现实数据往往是曲线的：

• 传染病早期感染者人数随时间指数增长；
• 人口规模随时间呈 S 型（Logistic）增长；
• 药物浓度随时间指数衰减；
• 酶促反应速率随底物浓度升高而趋于饱和（Michaelis-Menten 曲线）；
• 学习成本随产量上升而边际递减（对数曲线）。

这些曲线的形状由"机理"决定（例如"增长率与现有规模成正比"解微分方程就得到指数模型），函数形式是明确的，但其中的参数（增长率、饱和值等）未知，需要从数据中估计——这就是非线性回归要解决的问题。

严格定义：非线性回归研究模型

$y_i = f(x_i; \boldsymbol{\theta}) + \varepsilon_i, \qquad i = 1, 2, \dots, n$

其中：

• $x_i$ 为第 $i$ 个自变量的观测值（可以是向量）；
• $\boldsymbol{\theta} = (\theta_1, \theta_2, \dots, \theta_p)^\top$ 为待估参数向量，$p$ 为参数个数；
• $f$ 是关于参数 $\boldsymbol{\theta}$ 非线性的已知函数（函数形式由机理给定）；
• $\varepsilon_i$ 为随机误差，通常假设 $\varepsilon_i \sim N(0, \sigma^2)$ 独立同分布。

与线性回归的本质区别：判断"线性还是非线性"，看的是模型对参数是否线性，而不是对自变量。例如 $y = a + b \ln x$ 对参数 $a$、$b$ 是线性的（把 $\ln x$ 看成一个新自变量即可），所以它属于线性回归；而 $y = a e^{bx}$ 对参数 $b$ 是非线性的（$b$ 在指数位置上），无论怎样变换自变量都无法写成参数的线性组合，必须用非线性回归方法求解。

1.2 常用非线性模型速查表

模型名称	函数形式	曲线形态	典型应用（竞赛场景）
指数增长模型	$y = a e^{bx}$	$b>0$ 单调加速增长；$b<0$ 指数衰减	传染病早期传播、人口/经济初期增长、放射性衰变、药物消除
幂函数模型	$y = a x^{b}$	幂律增长/下降，双对数坐标下呈直线	异速生长（体重—代谢率）、规模效应、城市位序—规模律
对数模型	$y = a + b \ln x$	增长越来越慢（边际收益递减）	学习曲线、广告投入—销量关系、经济增长与资本投入
Logistic 生长曲线	$y = \dfrac{K}{1 + e^{a - bx}}$	S 型：先加速、后减速，渐近趋于饱和值 $K$	人口预测、传染病累计感染人数、新产品市场渗透率
Michaelis-Menten	$y = \dfrac{a x}{b + x}$	过原点、单调上升、渐近饱和于 $a$	酶促反应速率、吸附等温线（Langmuir 型）

参数含义提示：Logistic 曲线中 $K$ 为环境容量（饱和值），$a/b$ 决定拐点位置，$b$ 控制曲线陡峭度；Michaelis-Menten 中 $a$ 为最大反应速率 $V_{\max}$，$b$ 为半饱和常数 $K_m$（$x = b$ 时 $y = a/2$）。参数有物理意义，是非线性回归区别于多项式回归的核心价值。

1.3 关键概念：可线性化 vs 不可线性化

(1) 可线性化（变量变换法）：通过取对数、倒数等变换，把模型改造成线性形式：

原模型	变换方式	变换后形式	新"变量"
$y = a e^{bx}$	两边取自然对数	$\ln y = \ln a + b x$	$Y = \ln y$
$y = a x^b$	两边取自然对数	$\ln y = \ln a + b \ln x$	$Y = \ln y,\ X = \ln x$
$y = \dfrac{ax}{b+x}$	两边取倒数	$\dfrac{1}{y} = \dfrac{1}{a} + \dfrac{b}{a} \cdot \dfrac{1}{x}$	$Y = 1/y,\ X = 1/x$（Lineweaver-Burk 变换）
$y = \dfrac{K}{1+e^{a-bx}}$（$K$ 已知）	移项取对数	$\ln\dfrac{K-y}{y} = a - bx$	$Y = \ln\dfrac{K-y}{y}$

变换后直接套用线性最小二乘（如 np.polyfit），再反变换还原参数。这是传统手算/表格软件时代的做法，代码简单，但有重要缺陷（见下）。

(2) 不可线性化：找不到任何变换能将其化为线性形式的模型，例如：

• $y = a + b e^{cx}$（带常数项的三参数指数，无法消去常数项）；
• $y = a e^{bx} + c e^{dx}$（双指数/混合指数）；
• 高斯峰 $y = a \exp\left(-\dfrac{(x-\mu)^2}{2\sigma^2}\right)$（参数 $\mu$ 在平方内，取对数后仍是非线性）；
• Logistic 曲线中 $K$ 也未知时（$K$ 无法通过变换消去）。

这类模型只能用迭代的非线性最小二乘法求解。

(3) 取对数线性化会改变误差结构——必须牢记的关键区别

原模型通常假设误差是加性的：$y = a e^{bx} + \varepsilon$，即观测值在真值附近上下波动，波动幅度与 $y$ 的大小无关。

取对数后模型变为 $\ln y = \ln a + bx + \ln\!\left(1 + \frac{\varepsilon}{ae^{bx}}\right)$。当 $\varepsilon$ 相对较小时，$\ln(1+\delta) \approx \delta$，于是：

$\ln y \approx \ln a + bx + \frac{\varepsilon}{a e^{bx}}$

此时对数尺度上的"误差"是 $\varepsilon/(a e^{bx})$——相对误差：$y$ 小处相对误差被放大，$y$ 大处相对误差被缩小。而线性化后的最小二乘最小化的是

$\sum_{i=1}^{n} \left( \ln y_i - \ln \hat{y}_i \right)^2 \approx \sum_{i=1}^{n} \left( \frac{y_i - \hat{y}_i}{\hat{y}_i} \right)^2$

即相对误差平方和，与原始尺度的 $\sum (y_i - \hat{y}_i)^2$ 是两个不同的优化目标。结论：

• 若真实误差是加性的（恒定波动），取对数线性化会给小 $y$ 数据点过大的权重，参数估计产生偏差，且残差方差结构被扭曲；
• 若真实误差是乘性的（$y = ae^{bx} e^{\varepsilon}$，波动幅度与 $y$ 成正比——增长类数据常如此），取对数后误差恰好变成加性的，此时对数线性化反而是"正确"的做法。

一句话总结：变换之前先想清楚误差长什么样。 拿不准时，直接用非线性最小二乘（curve_fit），它按原始尺度最小化误差平方和，并且标准误、置信区间一并给出，不用做任何变换。

1.4 求解方法：非线性最小二乘

目标：找到使残差平方和最小的参数

$S(\boldsymbol{\theta}) = \sum_{i=1}^{n} r_i^2 = \sum_{i=1}^{n} \left( y_i - f(x_i; \boldsymbol{\theta}) \right)^2 = \|\mathbf{r}(\boldsymbol{\theta})\|^2$

与线性回归不同，令 $\partial S/\partial \theta_j = 0$ 得到的正规方程一般没有解析解（因为 $f$ 非线性），只能用迭代法。两个最经典的迭代算法：

(1) Gauss-Newton 法（牛顿法在最小二乘问题上的近似）

核心思想：在迭代点 $\boldsymbol{\theta}^{(k)}$ 处对 $f$ 做一阶泰勒展开，把非线性问题在局部"线性化"：

$f(x_i; \boldsymbol{\theta}) \approx f(x_i; \boldsymbol{\theta}^{(k)}) + \mathbf{J}_i(\boldsymbol{\theta}^{(k)}) \left( \boldsymbol{\theta} - \boldsymbol{\theta}^{(k)} \right)$

其中雅可比矩阵 $\mathbf{J}$ 的元素为 $J_{ij} = \dfrac{\partial f(x_i; \boldsymbol{\theta})}{\partial \theta_j}$。代入后每步只需解一个线性最小二乘问题，得到迭代公式：

$\boxed{\ \boldsymbol{\theta}^{(k+1)} = \boldsymbol{\theta}^{(k)} + \left( \mathbf{J}^\top \mathbf{J} \right)^{-1} \mathbf{J}^\top \mathbf{r}\!\left(\boldsymbol{\theta}^{(k)}\right)\ }$

优点：真值附近收敛很快（近似二次收敛）。缺点：$\mathbf{J}^\top \mathbf{J}$ 接近奇异时步长爆炸，初值不好时容易发散。

(2) Levenberg-Marquardt（LM）法——实际最常用，scipy 的默认算法

在 Gauss-Newton 基础上加一个阻尼项：

$\boxed{\ \boldsymbol{\theta}^{(k+1)} = \boldsymbol{\theta}^{(k)} + \left( \mathbf{J}^\top \mathbf{J} + \lambda \mathbf{I} \right)^{-1} \mathbf{J}^\top \mathbf{r}\!\left(\boldsymbol{\theta}^{(k)}\right)\ }$

• $\lambda \to 0$：趋近 Gauss-Newton，步长大、收敛快；
• $\lambda \to \infty$：趋近梯度下降，步长小、方向稳、更稳健；
• 实际实现动态调整 $\lambda$：本次迭代使 $S(\boldsymbol{\theta})$ 下降就减小 $\lambda$（放大步长），否则增大 $\lambda$（更保守地重试）。

收敛判据：$\|\boldsymbol{\theta}^{(k+1)} - \boldsymbol{\theta}^{(k)}\|$ 或 $|S^{(k+1)} - S^{(k)}|$ 小于阈值（如 $10^{-8}$），或达到最大迭代次数。

初值问题：非线性最小二乘是"从一个起点出发，走到离它最近的那个谷底"，只能保证找到局部最优。初值必须合理，常用策略：① 利用机理含义（如增长率 $b$ 的数量级）；② 先用线性化方法粗估一组参数作为初值；③ 多起点试验，取 SSE 最小者。

1.5 优缺点

优点	缺点
参数有明确机理含义（增长率、饱和值、环境容量），便于解释、便于写进论文	依赖初值，可能收敛到局部最优甚至发散
模型受机理约束、参数少，外推比纯数据驱动模型更可信	必须预先知道（或敢假设）函数形式，选错模型一切白搭
可由雅可比矩阵直接给出参数标准误、置信区间、拟合曲线的置信带	迭代算法对病态数据（参数强相关、矩阵近奇异）敏感
能处理无法线性化的模型（多参数指数、未知饱和值等）	参数过多时易过拟合、可辨识性差（见 7.2 节）

二、何时使用（适用场景与条件）

2.1 核心判断：有模型形式的先验

使用非线性回归的前提是你知道（或敢假设）曲线的函数形式，且这个形式来自：

• 机理推导：由微分方程解出来的曲线。例："感染人数增长率 ∝ 当前感染人数" → 解微分方程得指数模型；"增长率 ∝ 当前规模 × 剩余空间" → Logistic 曲线；
• 公认经验规律：酶动力学 Michaelis-Menten 方程、放射性衰变指数律、生物异速生长的幂律（Kleiber 定律）；
• 数据形态的强先验：散点图呈明显 S 型或饱和型，且理论上就该如此（如总人口不可能无限增长）。

2.2 竞赛典型场景

场景	建议模型	参数含义（论文话术素材）
传染病早期传播（COVID-19、SIR 早期近似）	$y = a e^{bx}$	$b$ 为指数增长率，可进一步换算基本再生数 $R_0$
人口/城市规模预测	$y = \dfrac{K}{1+e^{a-bx}}$	$K$ 为环境容量/最大承载规模
经济增长、GDP 预测	$y = a e^{bx}$ 或 $y = ax^b$	$b$ 为增长率/弹性系数
化学反应速率、酶动力学	$y = \dfrac{ax}{b+x}$	$a = V_{\max}$（最大速率），$b = K_m$（半饱和常数）
广告投入与销量（边际递减）	$y = a + b \ln x$	$b$ 为对数边际效应
药物浓度随时间衰减	$y = a e^{-bx}$	$b$ 为消除速率常数，半衰期 $t_{1/2} = \ln 2 / b$

2.3 使用前提（四条，缺一不可）

1. 模型形式合理：有机理或经验依据，且散点图形状与模型形态一致（先画图，再选模型）；
2. 初值接近真值：给不出合理初值时，先用线性化方法粗估，或多起点搜索（见 7.2 节第 1 条）；
3. 样本量足够：经验上 $n \ge 10p$（$p$ 为参数个数），否则标准误、置信区间不可靠；
4. 误差假设近似成立：通常假设等方差、独立、正态；不成立时考虑加权非线性最小二乘或先做变量变换。

2.4 不适用情形

• 无机理先验、只想光滑拟合/预测：用多项式回归、样条（spline）或机器学习方法（随机森林、SVR），不要硬套某个机理曲线——机理模型选错，比"黑箱"方法更危险；
• 模型参数不可辨识：如 $y = a b\, e^{cx}$ 中 $a$ 与 $b$ 的乘积是一个整体，无法分开估计（见 7.2 节第 2 条）；
• 数据只覆盖曲线的局部：如只有 Logistic 曲线拐点左侧的数据，饱和值 $K$ 的估计极不稳定——机理模型同样需要数据支撑；
• 响应变量是计数、二分类、比例：应改用广义线性模型（泊松回归、Logistic 回归等，见第八节）。

2.5 与多项式回归、广义线性模型的对比选择

方法	模型形式	参数意义	外推能力	何时选择
非线性回归	$y = f(x;\theta)$，形式由机理给定	明确（增长率、容量等）	较好（受机理约束）	有模型先验，参数需要解释与外推
多项式回归	$y = \beta_0 + \beta_1 x + \cdots + \beta_d x^d$	单项系数无物理意义	差（高次项外推迅速发散）	无先验，仅需局部拟合/插值
广义线性模型（GLM）	$g(E[y]) = \mathbf{X}\boldsymbol{\beta}$	线性效应、OR/RR 等	一般	响应为计数/二分类/比例，或需处理异方差

注意：多项式回归在数学上仍是参数线性的模型（线性回归的特例），它通过提高次数逼近任意光滑函数，属于"万能近似"；非线性回归则是"机理约束下的精准描述"。竞赛中若被问"为什么不用多项式"，标准回答是：多项式参数无机理含义、无法外推、高阶易震荡；而本模型参数有明确物理意义，外推受机理约束。

三、算法指标

本节所有指标中，$\hat{y}_i = f(x_i; \hat{\boldsymbol{\theta}})$ 为拟合值，$r_i = y_i - \hat{y}_i$ 为残差，$\bar{y}$ 为响应均值，$p$ 为参数个数，$n$ 为样本量。

3.1 拟合优度类

(1) SSE 残差平方和

$\text{SSE} = \sum_{i=1}^{n} r_i^2 = \sum_{i=1}^{n} \left( y_i - \hat{y}_i \right)^2$

• 取值范围 $[0, +\infty)$，越小越好；
• 量纲是 $y$ 量纲的平方，绝对值不能跨数据集比较；
• 它是优化目标本身——非线性最小二乘找的就是使 SSE 最小的参数。

(2) RMSE 均方根误差

$\text{RMSE} = \sqrt{\frac{\text{SSE}}{n}} = \sqrt{\frac{1}{n}\sum_{i=1}^{n}\left( y_i - \hat{y}_i \right)^2}$

• 取值范围 $[0, +\infty)$，越小越好；
• 与 $y$ 同量纲，可直接读作"平均误差约为多少"，适合在论文中与其他模型横向比较；
• 注意分母用 $n$ 还是 $n-p$ 各有说法，论文中写明口径即可（本文统一用 $n$）。

(3) R² 决定系数

$R^2 = 1 - \frac{\text{SSE}}{\text{SST}} = 1 - \frac{\sum_{i=1}^{n} (y_i - \hat{y}_i)^2}{\sum_{i=1}^{n} (y_i - \bar{y})^2}$

• 取值范围通常 $[0, 1]$（非线性回归中理论上可小于 0），越接近 1 拟合越好；
• 解释为"模型解释了响应变量总变异的百分之多少"；
• 警示：非线性回归的 R² 是类比定义，不能像线性回归那样做 F 检验；且 R² 高 ≠ 模型正确——残差有明显结构、参数不显著时 R² 依然可以很高。

3.2 参数推断类

(4) 参数估计值 ± 标准误

非线性最小二乘的参数协方差矩阵由收敛点处的雅可比矩阵给出（$J_{ij} = \partial f(x_i; \boldsymbol{\theta}) / \partial \theta_j$）：

$\widehat{\mathrm{Cov}}(\hat{\boldsymbol{\theta}}) \approx \hat{\sigma}^2 \left( \mathbf{J}^\top \mathbf{J} \right)^{-1}, \qquad \hat{\sigma}^2 = \frac{\text{SSE}}{n - p}$

参数标准误取对角线开方：

$\text{SE}(\hat{\theta}_j) = \sqrt{\left[ \widehat{\mathrm{Cov}}(\hat{\boldsymbol{\theta}}) \right]_{jj}}$

• 解读：SE 越小，参数估计越精确；SE 相对估计值过大说明参数不可靠（可辨识性差，见 7.2 节）；
• 这是大样本近似，$n$ 小时偏乐观；scipy 的 curve_fit 返回的 pcov 正是这个矩阵（默认 absolute_sigma=False）。

(5) 参数 95% 置信区间

$\hat{\theta}_j \pm t_{1-\alpha/2}(n-p) \cdot \text{SE}(\hat{\theta}_j), \qquad \alpha = 0.05$

• $t_{1-\alpha/2}(n-p)$ 是自由度 $n-p$ 的 t 分布上分位数（$n$ 大时约 1.96）；
• 解读：区间不包含 0 ⇒ 该参数在 5% 显著性水平下显著（即该效应"真实存在"）；区间很宽 ⇒ 参数估计不确定；
• 竞赛话术："参数 $b$ 的 95% 置信区间为 [0.385, 0.412]，不包含 0，说明增长趋势显著。"

(6) AIC 赤池信息准则

$\text{AIC} = n \ln\!\left( \frac{\text{SSE}}{n} \right) + 2p \qquad \left(\text{即 } n \ln \hat{\sigma}^2 + 2p + \text{常数}\right)$

• 取值范围为全体实数，越小越好；比较模型时只看差值：$\Delta\text{AIC} > 2$ 即有实质性差异；
• 第一项奖励拟合好，第二项惩罚参数多——用于不同非线性模型之间选型（指数 vs Logistic vs Michaelis-Menten）；
• 小样本（$n/p < 40$ 左右）用修正版 $\text{AICc} = \text{AIC} + \dfrac{2p(p+1)}{n-p-1}$；
• 特别注意：只能在同一尺度（同一响应变量）下比较 AIC——对 $y$ 拟合的模型与对 $\ln y$ 拟合的模型，AIC 不可直接比较（本文程序里两者都换算到原始尺度后再比）。

(7) 残差随机性检验

直观判断（论文中必须有图）：残差图应呈水平带状随机散布：无弯曲趋势（否则模型形式错误）、无喇叭形（否则异方差）、无连续同号段（否则自相关）。

定量辅助（无需额外库）：符号游程检验。记残差符号 "+"、"−" 交替的段数为游程数 $R$，$n_+$、$n_-$ 分别为正、负残差个数，残差随机时：

$Z = \frac{R - \mu_R}{\sigma_R} \sim N(0,1), \qquad \mu_R = \frac{2 n_+ n_-}{n} + 1, \qquad \sigma_R^2 = \frac{2 n_+ n_- (2 n_+ n_- - n)}{n^2 (n-1)}$

$|Z| < 1.96$ 表示在 5% 水平无法拒绝"残差随机"（游程过少表示残差成串同号，过多表示符号频繁交替）。也可用 Durbin-Watson 统计量 $DW = \dfrac{\sum_{i=2}^{n} (r_i - r_{i-1})^2}{\sum_{i=1}^{n} r_i^2}$，接近 2 表示无自相关。

3.3 指标汇总表

指标（中文）	公式	取值范围	解读要点
残差平方和 SSE	$\sum (y_i - \hat{y}_i)^2$	$[0, +\infty)$	优化目标，越小越好，量纲为 $y^2$
均方根误差 RMSE	$\sqrt{\text{SSE}/n}$	$[0, +\infty)$	平均误差水平，与 $y$ 同量纲，可跨模型比较
决定系数 R²	$1 - \text{SSE}/\text{SST}$	通常 $[0, 1]$	拟合优度；高 R² 不代表模型正确
参数估计 ± 标准误	$\hat{\theta}_j \pm \text{SE}(\hat{\theta}_j)$，SE 取自 $\hat{\sigma}^2(\mathbf{J}^\top\mathbf{J})^{-1}$ 对角线	SE ≥ 0	SE 小才可靠；SE 过大提示可辨识性差
参数 95% 置信区间	$\hat{\theta}_j \pm t_{0.975}(n-p)\cdot\text{SE}(\hat{\theta}_j)$	区间	不含 0 ⇒ 参数显著
AIC / AICc	$n\ln(\text{SSE}/n) + 2p$；AICc 再加 $\frac{2p(p+1)}{n-p-1}$	实数	越小越好；用于模型选型；须同尺度比较
残差随机性	游程检验 $Z$、$DW \approx 2$	$Z \sim N(0,1)$	$\|Z\| < 1.96$ 或 DW≈2 ⇒ 残差随机、模型合适

四、可视化图表

4.1 图表总览

图名	用途	关键解读点
① 原始数据 + 非线性拟合曲线 + 95% 置信带	展示拟合效果与拟合的不确定性	曲线应穿过数据点云中心；置信带在数据密集处窄、两端（尤其外推区）变宽；带越窄参数越确定
② 直接非线性拟合 vs 取对数线性化对比图	揭示两种估计方式的差异来源	原始坐标下两条曲线整体接近（低值区差异略明显）；半对数（ln y–x）坐标下两者都是直线，斜率截距的差异被放大显示；差异源于"原始尺度 vs 对数尺度"的优化目标不同
③ 残差图	检验模型假设（最重要的一张图）	残差围绕 0 水平带状随机散布＝好；弯曲＝模型形式错；喇叭形＝异方差；成串同号＝自相关
④ 拟合值—实际值散点图	直观展示整体预测能力	点沿对角线 $y=\hat{y}$ 密集分布＝好；系统性偏离对角线＝模型有偏；方差随 $y$ 增大属正常方差结构（或提示异方差）

4.2 好图与异常图的特征

好图特征：

• 拟合曲线穿行于数据点云中心，两侧点数大致均衡；
• 置信带窄而平滑，宽度在数据范围内变化自然；
• 残差图呈"随机噪声"外观：围绕 0 水平线均匀散布、无形状、无趋势；
• 拟合值—实际值散点紧贴对角线，无弯弓形偏离。

异常图特征（论文中出现必须解释或处理）：

• 拟合曲线与数据形态不匹配（如 S 型数据用指数模型硬拟合）→ 换模型；
• 置信带在数据端部急剧张开 → 外推区不可信，慎用模型做远端预测；
• 残差呈 U 形/倒 U 形 → 模型缺项或选错形式；
• 残差喇叭形（方差随 $\hat{y}$ 增大）→ 异方差：考虑加权最小二乘、取对数（乘性误差）或改用相对误差目标；
• 残差连续同号成串 → 自相关：时间序列数据考虑加入自回归结构或差分后再拟合。

五、符号说明

符号	含义	示例/单位
$x_i$	第 $i$ 个自变量的观测值	时间 $t = 3$（天）
$y_i$	第 $i$ 个响应变量的观测值	感染人数 128（人）
$f(x; \boldsymbol{\theta})$	非线性模型函数（形式由机理给定）	$f = a e^{bx}$
$\boldsymbol{\theta}$ / $p$	待估参数向量 / 参数个数	$\boldsymbol{\theta} = (a, b)^\top$，$p = 2$
$\hat{\boldsymbol{\theta}}$	参数的非线性最小二乘估计	$\hat{a} = 5.03$
$\varepsilon_i$	随机误差项	$\varepsilon_i \sim N(0, \sigma^2)$
$\sigma^2$ / $\hat{\sigma}^2$	误差方差 / 其无偏估计	$\hat{\sigma}^2 = \text{SSE}/(n-p)$
$r_i$	第 $i$ 个残差	$r_i = y_i - \hat{y}_i$
$\hat{y}_i$	拟合值（预测值）	$\hat{y}_i = f(x_i; \hat{\boldsymbol{\theta}})$
$\mathbf{J}$	雅可比矩阵，$J_{ij} = \partial f(x_i; \boldsymbol{\theta})/\partial \theta_j$	$n \times p$ 矩阵
$S(\boldsymbol{\theta})$	残差平方和目标函数	$S = \sum r_i^2$
$\lambda$	LM 算法的阻尼因子	$\lambda = 0.01$
SE / CI	标准误 / 置信区间	$\text{SE}(\hat{b}) = 0.006$
$t_{1-\alpha/2}(\nu)$	自由度 $\nu$ 的 t 分布上分位数	$t_{0.975}(78) \approx 1.99$
$K$	Logistic 曲线的饱和值（环境容量）	$K = 10000$（万人）
$a, b$	具体模型参数（随模型解释）	$a$：初始规模；$b$：增长率
SSE / RMSE / R²	残差平方和 / 均方根误差 / 决定系数	见第三节
AIC / AICc	赤池信息准则及其小样本修正	数值越小越好
$n$	样本量	$n = 80$

六、可运行程序（完整代码）

运行环境：Python 3.12，依赖 numpy、scipy、scikit-learn、matplotlib、pandas（无其他第三方库）。将所有代码块按顺序拼接保存为一个 .py 文件（如 nlr_demo.py），在命令行执行 python nlr_demo.py 即可；无图形界面的服务器环境可用 MPLBACKEND=Agg python nlr_demo.py 静默运行（仅保存图片、不弹窗）。

程序功能：① 生成合成数据（真实模型 $y = 5\,e^{0.4x}$，乘性噪声，$x \in [0, 10]$，$n = 80$，随机种子 42）；② 用 curve_fit 做非线性最小二乘并给出参数标准误、95% 置信区间；③ 用取对数线性化的传统做法拟合同一数据；④ 附手写 Gauss-Newton 迭代片段体现原理；⑤ 打印第三节全部指标；⑥ 绘制第四节全部 4 张图并保存到 figures/ 子目录（文件名前缀 nlr_）。

# -*- coding: utf-8 -*-
"""
《非线性回归》配套完整可运行程序（Python 3.12）
数据：合成指数增长数据  y = 5·e^{0.4x}·e^{ε}（乘性噪声），x∈[0,10]，n=80
方法一：scipy.optimize.curve_fit 非线性最小二乘（Levenberg-Marquardt 算法）
方法二：取对数线性化（np.polyfit 对 ln y 做一元线性回归）
输出：第三节全部指标（SSE/RMSE/R²/标准误/95%CI/AIC/残差检验）+ 第四节 4 张图
依赖：numpy, scipy, scikit-learn, matplotlib, pandas（无其他第三方库）
"""

# ========== 0. 导入库与全局设置 ==========
import os
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy import stats
from scipy.optimize import curve_fit
from sklearn.metrics import r2_score, mean_squared_error

# 中文字体设置（防止图片中中文显示为方块）
plt.rcParams["font.sans-serif"] = ["PingFang SC", "Arial Unicode MS", "SimHei"]
plt.rcParams["axes.unicode_minus"] = False   # 让负号正常显示

# 图片保存目录：脚本运行目录下的 figures/ 子目录
os.makedirs("figures", exist_ok=True)

第 1 步：生成模拟数据。 固定随机种子，保证任何人运行得到完全一致的结果；噪声取乘性（对数尺度加性、约 15% 相对误差），这是增长类数据最常见的误差结构，也便于后面讲解"误差结构决定方法选择"。

# ========== 1. 生成模拟数据（真实模型：指数增长 + 乘性噪声） ==========
np.random.seed(42)                            # 固定随机种子：运行结果完全可复现
n = 80                                        # 样本量
x = np.linspace(0.0, 10.0, n)                 # 自变量均匀取 80 个点
TRUE_A, TRUE_B = 5.0, 0.4                     # 真实参数：y = 5·e^{0.4x}
SIGMA_LN = 0.15                               # 对数尺度噪声标准差 ≈ 15% 相对误差
eps = np.random.normal(0.0, SIGMA_LN, n)      # 正态随机扰动（对数尺度）
y = TRUE_A * np.exp(TRUE_B * x) * np.exp(eps) # 乘性噪声：波动幅度与 y 大小成正比

print("=" * 66)
print("《非线性回归》演示：指数增长模型 y = a·e^(b·x)")
print(f"真实参数: a={TRUE_A}, b={TRUE_B}；乘性噪声 σ(ln 尺度)={SIGMA_LN}；n={n}")
print("-" * 66)
print("生成数据预览（前 5 行）：")
print(pd.DataFrame({"x": x, "y": y}).head().round(4).to_string(index=False))

第 2 步：方法一——curve_fit 非线性最小二乘。 注意 p0 初值必须合理（按机理猜：初始规模约 3、增长率约 0.2），给错初值可能收敛到局部最优或发散（见 7.2 节）。返回的 pcov 即协方差矩阵 $\hat{\sigma}^2(\mathbf{J}^\top\mathbf{J})^{-1}$。

# ========== 2. 方法一：scipy.optimize.curve_fit（非线性最小二乘） ==========
def exp_model(x, a, b):
    """指数增长模型 y = a·e^(b·x)：a 为初始规模，b 为增长率"""
    return a * np.exp(b * x)

# curve_fit 默认 method='lm'，即 Levenberg-Marquardt 算法；p0 必须给合理初值
p0 = [3.0, 0.2]                               # 初值：按机理猜测（初始规模 3，增长率 0.2）
popt, pcov = curve_fit(exp_model, x, y, p0=p0)
a_nls, b_nls = popt

# 参数标准误：pcov ≈ (JᵀJ)^{-1}·σ̂²（curve_fit 默认 absolute_sigma=False）
se_nls = np.sqrt(np.diag(pcov))
se_a_nls, se_b_nls = se_nls

# 参数 95% 置信区间：θ̂ ± t_{0.975}(n-p)·SE
p = 2                                         # 参数个数
t_crit = stats.t.ppf(1 - 0.05 / 2, df=n - p)  # t 分布临界值（n 大时 ≈ 1.96）
ci_a_nls = (a_nls - t_crit * se_a_nls, a_nls + t_crit * se_a_nls)
ci_b_nls = (b_nls - t_crit * se_b_nls, b_nls + t_crit * se_b_nls)

y_hat_nls = exp_model(x, *popt)               # 拟合值
resid_nls = y - y_hat_nls                     # 残差（原始尺度）

第 3 步：计算拟合优度指标，并用 sklearn 交叉验证手算结果（两者应完全一致）。

# ========== 3. 指标计算（方法一） ==========
SSE_nls = float(np.sum(resid_nls ** 2))                     # 残差平方和
RMSE_nls = float(np.sqrt(SSE_nls / n))                      # 均方根误差
SST = float(np.sum((y - y.mean()) ** 2))                    # 总离差平方和
R2_nls = 1.0 - SSE_nls / SST                                # 决定系数
AIC_nls = n * np.log(SSE_nls / n) + 2 * p                   # AIC（略去公共常数）
AICc_nls = AIC_nls + 2 * p * (p + 1) / (n - p - 1)          # 小样本修正 AICc

# 用 sklearn 交叉验证上面手算的指标（结果应完全一致）
R2_sk = float(r2_score(y, y_hat_nls))
RMSE_sk = float(np.sqrt(mean_squared_error(y, y_hat_nls)))

第 4 步：方法二——取对数线性化。 对 $y = ae^{bx}e^{\varepsilon}$ 两边取对数得 $\ln y = \ln a + bx + \varepsilon$，乘性噪声变成加性噪声，于是可在 $\ln y$ 上直接做一元线性回归（np.polyfit），再反变换还原参数。注意：反变换 exp(ln ŷ) 预测的是 $y$ 的中位数，预测均值需乘偏置修正因子 $\exp(\hat{\sigma}^2/2)$。

# ========== 4. 方法二：取对数线性化（传统做法） ==========
# 两边取对数：ln y = ln a + b·x + ε_ln。原乘性噪声变成对数尺度的加性噪声，
# 于是可在 ln y 上直接做一元线性回归（np.polyfit），再反变换还原参数。
Y = np.log(y)
coef, cov_log = np.polyfit(x, Y, 1, cov=True)   # 一元线性回归 + 协方差矩阵
b_log, ln_a_log = coef[0], coef[1]              # 斜率 = b，截距 = ln a
a_log = float(np.exp(ln_a_log))                 # 反变换还原 a

# 对数尺度参数标准误；注意 np.polyfit 的协方差按 [斜率, 截距] 顺序排列
se_b_log = float(np.sqrt(cov_log[0, 0]))   # 斜率 b 的标准误
se_ln_a = float(np.sqrt(cov_log[1, 1]))    # 截距 ln a 的标准误
se_a_log = float(a_log * se_ln_a)          # a 的标准误用 delta 方法近似：SE(a) ≈ a·SE(ln a)
ci_a_log = (a_log - t_crit * se_a_log, a_log + t_crit * se_a_log)
ci_b_log = (b_log - t_crit * se_b_log, b_log + t_crit * se_b_log)

# 拟合值与残差：原始尺度 + 对数尺度（两种残差都要看）
y_hat_log = a_log * np.exp(b_log * x)
resid_log_orig = y - y_hat_log                  # 原始尺度残差
resid_log_ln = Y - (ln_a_log + b_log * x)       # 对数尺度残差

# 注意：反变换 exp(ln ŷ) 预测的是 y 的中位数；预测均值需乘偏置修正因子 exp(σ̂²/2)
sigma2_ln = float(np.sum(resid_log_ln ** 2) / (n - p))   # 对数尺度残差方差估计
bias_corr = float(np.exp(sigma2_ln / 2))                 # 均值预测偏置修正因子

# 对数线性化模型的指标（统一在"原始尺度"上算，才能与 curve_fit 公平比较）
SSE_log = float(np.sum(resid_log_orig ** 2))
RMSE_log = float(np.sqrt(SSE_log / n))
R2_log = 1.0 - SSE_log / SST
AIC_log = n * np.log(SSE_log / n) + 2 * p
AICc_log = AIC_log + 2 * p * (p + 1) / (n - p - 1)

第 5 步：手写 Gauss-Newton 迭代核心（约 10 行），体现"逐次线性化 + 解正规方程"的原理，收敛结果应与 curve_fit 完全一致（LM 在 $\lambda \to 0$ 时退化为 Gauss-Newton）。注意纯 Gauss-Newton 对初值敏感：把初值改成 [3.0, 0.2] 它就会发散（读者可自行试验），这正是 LM 需要阻尼项 $\lambda$ 的原因。

# ========== 5. 手写 Gauss-Newton 迭代（约 10 行核心，体现算法原理） ==========
# 迭代公式：θ ← θ + (JᵀJ)^{-1}·Jᵀ·r，其中 r 为残差向量，J 为雅可比矩阵
# 注意：纯 Gauss-Newton 对初值敏感——若初值取 [3.0, 0.2] 会发散（读者可自行试验），
# 这正是 LM 算法要加阻尼 λ 的原因；这里取较接近真值的初值 [4.0, 0.3]。
theta = np.array([4.0, 0.3])                    # 初值：需较接近真值（否则 GN 发散）
for it in range(50):                            # 最多迭代 50 次
    r = y - exp_model(x, *theta)                # 残差向量 r(θ)
    J = np.column_stack([np.exp(theta[1] * x),  #  ∂f/∂a = e^{b x}
                         theta[0] * x * np.exp(theta[1] * x)])  # ∂f/∂b = a x e^{b x}
    delta = np.linalg.solve(J.T @ J, J.T @ r)   # 解正规方程 (JᵀJ)·Δθ = Jᵀr
    theta = theta + delta                       # 参数更新
    if np.linalg.norm(delta) < 1e-12:           # 收敛判断：步长足够小即停
        break
print("-" * 66)
print(f"手写 Gauss-Newton：迭代 {it + 1} 次收敛 → a = {theta[0]:.4f}, b = {theta[1]:.4f}")
print("（与 curve_fit 结果一致，验证了 LM 算法在 λ→0 时退化为 Gauss-Newton）")

第 6 步：残差随机性检验（符号游程检验，纯 numpy 实现）+ 打印第三节全部指标。

# ========== 6. 残差随机性检验：符号游程检验（纯 numpy 实现） ==========
def runs_test(resid):
    """H0: 残差符号随机排列。返回 (游程数 R, 标准化统计量 Z)，|Z|<1.96 接受 H0"""
    s = np.sign(resid)
    s = s[s != 0]                               # 去掉恰好为 0 的残差
    n1, n2 = int(np.sum(s > 0)), int(np.sum(s < 0))
    if n1 == 0 or n2 == 0:
        return None, np.nan
    runs = 1 + int(np.sum(s[1:] != s[:-1]))     # 数游程：符号切换次数 + 1
    mu = 2.0 * n1 * n2 / (n1 + n2) + 1          # 随机排列下游程数的期望
    var = 2.0 * n1 * n2 * (2.0 * n1 * n2 - n1 - n2) / ((n1 + n2) ** 2 * (n1 + n2 - 1))
    z = (runs - mu) / np.sqrt(var)
    return runs, z

runs_nls, z_nls = runs_test(resid_nls)
runs_log, z_log = runs_test(resid_log_ln)

# ========== 7. 打印第三节全部指标 ==========
print("=" * 66)
print("方法一：curve_fit 非线性最小二乘（原始尺度最小化 SSE）")
print("-" * 66)
print(f"  参数估计      : a = {a_nls:.4f} ± {se_a_nls:.4f}    b = {b_nls:.4f} ± {se_b_nls:.4f}")
print(f"  95% 置信区间  : a ∈ [{ci_a_nls[0]:.4f}, {ci_a_nls[1]:.4f}]    "
      f"b ∈ [{ci_b_nls[0]:.4f}, {ci_b_nls[1]:.4f}]")
print(f"  SSE = {SSE_nls:.2f}    RMSE = {RMSE_nls:.4f}    R² = {R2_nls:.5f}")
print(f"  AIC = {AIC_nls:.2f}    AICc = {AICc_nls:.2f}")
print(f"  残差游程检验  : 游程数 R = {runs_nls}, Z = {z_nls:+.3f}  (|Z|<1.96 → 残差随机)")
print(f"  sklearn 交叉验证: R² = {R2_sk:.5f}, RMSE = {RMSE_sk:.4f}（与手算一致）")
print("=" * 66)
print("方法二：取对数线性化（ln y 尺度最小化 SSE）")
print("-" * 66)
print(f"  ln 尺度回归   : ln y = {ln_a_log:.4f} + {b_log:.4f}·x")
print(f"  反变换参数    : a = e^({ln_a_log:.4f}) = {a_log:.4f} ± {se_a_log:.4f}    "
      f"b = {b_log:.4f} ± {se_b_log:.4f}")
print(f"  95% 置信区间  : a ∈ [{ci_a_log[0]:.4f}, {ci_a_log[1]:.4f}]    "
      f"b ∈ [{ci_b_log[0]:.4f}, {ci_b_log[1]:.4f}]")
print(f"  原始尺度指标  : SSE = {SSE_log:.2f}    RMSE = {RMSE_log:.4f}    "
      f"R² = {R2_log:.5f}    AIC = {AIC_log:.2f}    AICc = {AICc_log:.2f}")
print(f"  对数尺度残差  : 游程数 R = {runs_log}, Z = {z_log:+.3f}")
print(f"  对数尺度 σ̂² = {sigma2_ln:.4f}，均值预测偏置修正因子 exp(σ̂²/2) = {bias_corr:.4f}")
print("=" * 66)
print("两种方法对比（指标均在原始尺度计算，才可比较）：")
compare = pd.DataFrame({
    "方法": ["curve_fit 非线性最小二乘", "取对数线性化"],
    "a 估计": [a_nls, a_log],
    "b 估计": [b_nls, b_log],
    "RMSE": [RMSE_nls, RMSE_log],
    "R²": [R2_nls, R2_log],
    "AIC": [AIC_nls, AIC_log],
})
print(compare.round(4).to_string(index=False))
print("差异原因：curve_fit 在原始尺度最小化 SSE，原始尺度 RMSE 必然更小（优化目标所致）；")
print("对数线性化在对数尺度最小化，相当于最小化相对误差，其对数尺度残差更均匀。")
print("若数据误差本质是乘性的（本数据正是），对数线性化的误差假设更贴合生成机制；")
print("若误差是加性的，则应信任 curve_fit 的结果。")

第 7 步：图一——原始数据 + 非线性拟合曲线 + 95% 置信带。 置信带由协方差传播公式 $\mathrm{Cov}(\hat{y}) = \mathbf{J}\,\mathrm{pcov}\,\mathbf{J}^\top$ 在细网格上计算，展示"数据密集处带窄、外推区带宽"的典型特征。

# ========== 8. 图一：原始数据 + 非线性拟合曲线 + 95% 置信带 ==========
x_grid = np.linspace(x.min(), x.max(), 300)     # 细网格：画光滑曲线用
y_grid = exp_model(x_grid, *popt)               # 细网格上的拟合值
# 置信带：由协方差传播公式 Cov(ŷ) = J·pcov·Jᵀ 得到每个网格点的 SE(ŷ)
J_grid = np.column_stack([np.exp(b_nls * x_grid),
                          a_nls * x_grid * np.exp(b_nls * x_grid)])
se_grid = np.sqrt(np.diag(J_grid @ pcov @ J_grid.T))
band = t_crit * se_grid                         # 95% 置信带半宽

fig1, ax1 = plt.subplots(figsize=(9, 6))
ax1.scatter(x, y, s=22, color="#4C72B0", alpha=0.8, label="观测数据 (n=80)")
ax1.plot(x_grid, y_grid, color="#C44E52", lw=2.5, label="拟合曲线 y = a·e^(b·x)")
ax1.fill_between(x_grid, y_grid - band, y_grid + band,
                 color="#C44E52", alpha=0.15, label="95% 置信带（均值曲线）")
ax1.set_xlabel("x"); ax1.set_ylabel("y")
ax1.set_title("图1 非线性最小二乘拟合曲线与 95% 置信带")
ax1.legend()
fig1.tight_layout()
fig1.savefig(os.path.join("figures", "nlr_fit_band.png"), dpi=150)
print("已保存图片: figures/nlr_fit_band.png")

第 8 步：图二——两种方法的对比。 左图原始坐标下两条曲线整体接近、低值区略有差异；右图半对数坐标下两者都是直线，斜率截距的差异被放大显示——这就是"线性化与不线性化"的全部差异所在。

# ========== 9. 图二：直接非线性拟合 vs 取对数线性化（对比） ==========
fig2, (ax2a, ax2b) = plt.subplots(1, 2, figsize=(13, 5.5))
# 左图：原始坐标——两条曲线整体接近，低值区略有差异
ax2a.scatter(x, y, s=18, color="#4C72B0", alpha=0.75, label="观测数据")
ax2a.plot(x_grid, y_grid, color="#C44E52", lw=2, label="curve_fit（原始尺度 LS）")
ax2a.plot(x_grid, a_log * np.exp(b_log * x_grid), color="#55A868", lw=2,
          ls="--", label="对数线性化（ln y 尺度 LS）")
ax2a.set_xlabel("x"); ax2a.set_ylabel("y")
ax2a.set_title("原始坐标：两曲线接近，低值区略有差异")
ax2a.legend()
# 右图：半对数坐标——两者都是直线，斜率/截距的微小差异被放大显示
ax2b.scatter(x, Y, s=18, color="#4C72B0", alpha=0.75, label="ln y 观测数据")
ax2b.plot(x_grid, np.log(exp_model(x_grid, *popt)), color="#C44E52", lw=2,
          label="curve_fit 对应的 ln ŷ")
ax2b.plot(x_grid, ln_a_log + b_log * x_grid, color="#55A868", lw=2,
          ls="--", label="对数线性化的直线")
ax2b.set_xlabel("x"); ax2b.set_ylabel("ln y")
ax2b.set_title("半对数坐标：斜率截距差异被放大显示")
ax2b.legend()
fig2.tight_layout()
fig2.savefig(os.path.join("figures", "nlr_log_compare.png"), dpi=150)
print("已保存图片: figures/nlr_log_compare.png")

第 9 步：图三——残差图（2×2 诊断面板）。 上半部是 curve_fit 在原始尺度的残差（喇叭形 → 乘性噪声导致的异方差）；下半部是对数线性化在对数尺度的残差（随机均匀 → 同方差）。同一份数据，在不同尺度下的残差结构完全不同——这正是"取对数会改变误差结构"的直观证据。

# ========== 10. 图三：残差图（2×2 诊断面板） ==========
fig3, axes = plt.subplots(2, 2, figsize=(13, 9))
# (a) curve_fit：残差 vs x（原始尺度）
axes[0, 0].scatter(x, resid_nls, s=18, color="#4C72B0")
axes[0, 0].axhline(0, color="k", lw=1)
axes[0, 0].set_xlabel("x"); axes[0, 0].set_ylabel("残差 e_i")
axes[0, 0].set_title("(a) curve_fit：残差 vs x（喇叭形 → 乘性噪声/异方差）")
# (b) curve_fit：残差 vs 拟合值（原始尺度）
axes[0, 1].scatter(y_hat_nls, resid_nls, s=18, color="#4C72B0")
axes[0, 1].axhline(0, color="k", lw=1)
axes[0, 1].set_xlabel("拟合值 ŷ"); axes[0, 1].set_ylabel("残差 e_i")
axes[0, 1].set_title("(b) curve_fit：残差 vs 拟合值")
# (c) 对数线性化：对数尺度残差 vs x
axes[1, 0].scatter(x, resid_log_ln, s=18, color="#55A868")
axes[1, 0].axhline(0, color="k", lw=1)
axes[1, 0].set_xlabel("x"); axes[1, 0].set_ylabel("ln 尺度残差")
axes[1, 0].set_title("(c) 对数线性化：ln 尺度残差 vs x（随机均匀 → 同方差）")
# (d) 对数线性化：对数尺度残差 vs 拟合值
axes[1, 1].scatter(ln_a_log + b_log * x, resid_log_ln, s=18, color="#55A868")
axes[1, 1].axhline(0, color="k", lw=1)
axes[1, 1].set_xlabel("ln 尺度拟合值"); axes[1, 1].set_ylabel("ln 尺度残差")
axes[1, 1].set_title("(d) 对数线性化：ln 尺度残差 vs 拟合值")
fig3.suptitle("图3 残差诊断：两种方法在不同尺度下的残差结构", fontsize=14)
fig3.tight_layout()
fig3.savefig(os.path.join("figures", "nlr_residuals.png"), dpi=150)
print("已保存图片: figures/nlr_residuals.png")

第 10 步：图四——拟合值—实际值散点图。 点越贴近对角线 $y = \hat{y}$，预测越准；两种方法的点几乎重叠，说明两者预测能力接近。

# ========== 11. 图四：拟合值-实际值散点图 ==========
fig4, ax4 = plt.subplots(figsize=(8, 7))
ax4.scatter(y, y_hat_nls, s=22, color="#4C72B0", alpha=0.8, label="curve_fit")
ax4.scatter(y, y_hat_log, s=22, color="#55A868", alpha=0.6, marker="s",
            label="对数线性化")
lims = [0, float(y.max()) * 1.05]
ax4.plot(lims, lims, color="#C44E52", lw=2, ls="--", label="y = ŷ 理想线")
ax4.set_xlim(lims); ax4.set_ylim(lims)
ax4.set_xlabel("实际值 y"); ax4.set_ylabel("拟合值 ŷ")
ax4.set_title("图4 拟合值-实际值散点（越贴近对角线越好）")
ax4.legend()
ax4.set_aspect("equal")
fig4.tight_layout()
fig4.savefig(os.path.join("figures", "nlr_pred_actual.png"), dpi=150)
print("已保存图片: figures/nlr_pred_actual.png")

第 11 步：展示图片。 本地交互环境（命令行/Jupyter）会弹出 4 个图形窗口；无图形界面的服务器环境下仅产生一条警告，不影响运行与保存。

# ========== 12. 展示图片 ==========
# 本地交互环境（命令行/Jupyter）下弹出 4 个图形窗口；
# 无图形界面的服务器环境下该行仅产生一条警告，图片已保存为文件
print("=" * 66)
print("运行结束：4 张图已保存至 figures/ 目录（前缀 nlr_）")
try:
    plt.show()
except Exception as e:
    print("当前环境无图形界面，跳过弹窗：", e)

七、结果解读与注意事项

7.1 以程序输出为例的完整解读

运行第六节程序，得到如下结果（随机种子固定为 42，数值完全可复现）：

（1）参数估计与统计推断

• 方法一（curve_fit）：$\hat{a} = 5.71 \pm 0.72$（真值 5.00）、$\hat{b} = 0.384 \pm 0.014$（真值 0.400）。$b$ 的 95% 置信区间 $[0.356, 0.413]$ 不包含 0，说明增长效应统计显著；$\hat{a}$ 的区间 $[4.28, 7.14]$ 覆盖真值 5.0；
• 方法二（对数线性化）：$\hat{a} = 4.87 \pm 0.16$、$\hat{b} = 0.401 \pm 0.006$——更接近真值，且标准误明显更小（$b$ 的标准误 0.006 对比 0.014）。原因是本例误差是乘性的，对数尺度回归恰好在"正确尺度"上做最小二乘，效率更高；
• 解读思路：先看标准误与估计值的相对大小（变异系数 $\text{SE}/\hat{\theta}$：$b$ 约 3.7%，$a$ 约 13%——$a$、$b$ 强相关，噪声在两个参数间分配），再看置信区间是否跨 0。

（2）拟合优度

• R² = 0.952：模型解释了响应变量约 95% 的变异；
• RMSE = 15.78（y 的范围约 5~250，平均相对误差约 15% 量级，与设定的噪声水平一致）；
• AIC = 445.43、AICc = 445.58；
• 提醒：R² 高只能说明"曲线贴数据"，不能说明"模型选得对"——必须结合残差图（图3）与机理判断。

（3）两种方法对比

• 参数有可见差异：$\hat{b}$ 差 0.017（0.384 vs 0.401），$\hat{a}$ 差约 17%（5.71 vs 4.87）。本例噪声为乘性，原始尺度最小二乘被大 $y$ 点主导，采样波动被放大（$b$ 的标准误 0.014 vs 0.006）；对数线性化的估计更接近真值；
• 原始尺度 RMSE：curve_fit（15.78）略小于对数线性化（15.95）——必然结果，因为它的优化目标正是原始尺度 SSE；
• AIC：445.43 vs 447.15，$\Delta\text{AIC} = 1.7 < 2$，无实质差异——两种方法在"预测能力"上接近；
• 真正的区别在残差结构（图3）：curve_fit 的原始尺度残差呈喇叭形（方差随拟合值增大），对数线性化的对数尺度残差随机均匀。由于本数据误差本质是乘性的，对数线性化的误差假设更贴合生成机制，参数估计应更信任它（或对 curve_fit 做加权/改用对数尺度目标）——这正是第三节强调"先想清楚误差结构"的原因。

（4）残差检验

• 两种方法的游程检验分别为 $Z = -1.15$ 与 $Z = +0.23$，均满足 $|Z| < 1.96$：残差符号随机，无自相关迹象；
• 若你的数据出现 $|Z| > 1.96$ 或残差图有明显形状，先检查模型形式，再检查数据是否有序（时间序列）或存在异方差。

7.2 常见坑（竞赛中高频失分点）

坑 1：初值敏感性 → 收敛到局部最优甚至发散。

curve_fit 是从初值出发向"最近的谷底"走。若初值离真值太远（如把增长率初值给成 -0.5），可能收敛到错误参数组合，或直接报"迭代不收敛/协方差矩阵奇异"的错误。对策：① 用机理常识定初值数量级；② 先用线性化方法粗估一组参数当初值；③ 多试几组初值，取 SSE 最小且参数合理的那组；④ 必要时用 bounds 参数限定参数范围（如 $K > 0$）。

坑 2：参数可辨识性（过参数化）。

若模型参数之间存在函数关系（如 $y = ab\,e^{cx}$ 中 $a$ 与 $b$ 只能以乘积 $ab$ 出现），或数据只覆盖曲线的一段（如只有 Logistic 拐点左侧数据，$K$ 无法确定），则参数"识别不出来"：表现为 SE 巨大、置信区间极宽、不同初值收敛到完全不同但 SSE 几乎一样的解。对策：合并参数（令 $c = ab$）、固定部分参数、或补充覆盖关键区域的数据。

坑 3：线性化改变误差结构导致估计偏差。

对加性噪声数据盲目取对数拟合，相当于给小 $y$ 数据点过大的权重，参数估计有偏（小值区残差被过度放大）；对乘性噪声数据硬用原始尺度最小二乘，残差会呈喇叭形（异方差），虽然参数通常仍近似无偏，但标准误不再可信。对策：观察残差图判断误差结构，再选择方法——加性误差用 curve_fit，乘性误差用对数线性化（或加权最小二乘）。

坑 4：外推风险。

指数模型外推时 $y$ 以 $e^{b}$ 倍速膨胀，10 个时间单位后就是天文数字；Logistic 模型在数据不足拐点时 $K$ 的置信区间极宽。图1 的置信带在数据右端快速张开，直观展示了这一点。论文对策：只做短期外推；外推时同时给出预测值的置信区间；对指数模型，最好说明其"只适用于早期/无资源约束阶段"的适用边界。

坑 5：对数变换后的偏置与反变换错误。

对数尺度上拟合的是 $\ln y$ 的均值，反变换 $\exp(\ln \hat{y})$ 给出的是 $y$ 的中位数而非均值：$\mathbb{E}[y] = \exp(\ln\hat{y})\cdot \exp(\sigma^2/2)$。忽略修正因子会系统性低估预测值（本程序打印的修正因子约为 1.01，噪声小时影响不大，噪声大时不可忽略）。另外，忘记反变换（直接拿对数尺度的拟合值当预测值）是竞赛中的低级错误。

坑 6：混淆置信带与预测带。

图1 的 95% 置信带描述均值曲线的不确定性，$n \to \infty$ 时趋于 0；而单个新观测的 95% 预测带还要额外加上误差方差 $\sigma^2$，永远不会消失。论文中注明画的是哪一种，别用置信带冒充预测带。

7.3 竞赛论文写作话术模板

模型选择段：

"由传染病传播机理可知，早期感染人数近似服从指数增长规律，故采用指数增长模型 $y = a e^{bx}$ 描述……其中 $b$ 为增长率，具有明确的流行病学意义。相较于缺乏机理约束的多项式拟合，该模型参数可解释、外推受机理约束，更适合本问题。"

拟合结果段：

"基于机理，采用指数增长模型 $y = a e^{bx}$ 对数据（$n = 80$）进行非线性最小二乘拟合（Levenberg-Marquardt 算法），得 $\hat{a} = 5.71$（95% CI [4.28, 7.14]）、$\hat{b} = 0.384$（95% CI [0.356, 0.413]）。参数 95% 置信区间均不包含 0，说明初始规模与增长效应均显著；决定系数 $R^2 = 0.952$，RMSE = 15.78，模型解释了约 95% 的变异。"

模型比较段：

"为验证模型选择的合理性，将指数模型与 Logistic 模型、幂函数模型进行 AIC 比较，指数模型 AIC 最小（$\Delta\text{AIC} > 2$），且其参数均有明确物理意义，故选择指数模型。"

残差诊断段：

"拟合残差围绕零线随机散布，符号游程检验 $Z = -1.15$（$|Z| < 1.96$），无法拒绝残差随机性的原假设；残差方差随拟合值增大呈喇叭形，提示误差为乘性结构，与数据生成机制一致，模型假设合理。"

八、延伸阅读

• 广义线性模型（GLM）：当响应变量是计数（传染病日新增病例 → 泊松回归）、二分类（是否患病 → Logistic 回归）或存在方差—均值关系时，非线性回归的正态误差假设失效，应改用 GLM：$g(\mathbb{E}[y]) = \mathbf{X}\boldsymbol{\beta}$，通过链接函数 $g$（log、logit 等）把期望值映射到线性预测子。R 的 glm() 是标准工具；Python 下 sklearn 提供 PoissonRegressor、LogisticRegression（statsmodels 需另行安装）。注意：GLM 对线性预测子是线性的，与本文"参数非线性"是两码事，二者可组合（非线性链接 + 非线性预测子）。
• 广义加性模型（GAM）：把线性项换成光滑函数之和：$y = f_1(x_1) + f_2(x_2) + \varepsilon$，$f_j$ 用样条基表示，属于"半参数"方法——既不需要预设具体非线性形式，又比纯黑箱模型可解释。适合探索性建模；Python 可用 pyGAM 包（需另行安装）。竞赛中若机理不清但需要可解释的非线性效应，GAM 是比多项式回归更好的选择。
• 分段回归（piecewise regression / 变点模型）：机理可能在某个阈值处切换（如疫情管控前后增长率突变、酶促反应在不同温度区间的速率机制不同），此时用分段模型：$y = a_1 + b_1 x\ (x \le c),\ y = a_2 + b_2 x\ (x > c)$，其中变点 $c$ 也是待估参数（可用网格搜索 + 分段拟合，或专门的 piecewise-regression 包）。竞赛中"政策拐点""阶段划分"类问题常用。
• 贝叶斯非线性建模：给参数设先验分布（如 $K \sim U(0, 10^6)$、$b \sim N(0, 1)$），用 MCMC 采样（PyMC、Stan）得到参数的后验分布，直接给出任意量（含外推预测）的不确定性，且天然规避局部最优问题。代价是计算量大、需要设定先验。竞赛中用于小样本、强不确定性或需要传播不确定性的场景。
• 经典参考书：Bates & Watts《Nonlinear Regression Analysis and Its Applications》（1988）、Seber & Wild《Nonlinear Regression》（2003）；实践文档：scipy 官方文档 scipy.optimize.curve_fit 页面（含 LM 算法说明与协方差矩阵解释）。