1.4 第一个科学计算程序

1.4 第一个科学计算程序

本节摘要：把前三节的知识串成一条完整链路：读入一份带缺失值的实验数据，用插值补全缺口，用最小二乘拟合出物理规律，再用数值积分算出曲线下的总量。这段程序只有几十行，却覆盖了后续各章的核心模块，是你 SciPy 之旅的第一个里程碑。

阅读收获

阅读完本节，你应当能够：

1. 独立完成"数据读取 → 插值补缺 → 曲线拟合 → 数值积分"的完整流程；
2. 理解 curve_fit 的参数含义与协方差输出的作用；
3. 用代码组织意识（函数划分、分步输出）管理一个小的科学计算任务；
4. 说出每个环节如果不做会产生什么后果。

一、问题与直觉

假设你是一位生物实验员，连续 24 小时每小时测一次细菌培养皿的菌落量，得到一串随时间增长的观测点。但实验中途有两小时仪器故障，数据出现了两个缺口；而且你怀疑这个增长符合指数规律，想验证一下，还想知道 24 小时内的总增长量。

这就是一个典型的真实科学计算任务：数据有缺失（要插值）、有噪声（要拟合）、要总量（要积分）。三个问题分别对应 SciPy 的三个模块——interpolate、optimize、integrate。我们一次性把它们走通。动手优先的路线就是这样：先不管每个模块的边边角角，把整条流水线跑起来，再逐个深入研究。

二、核心原理

2.1 问题拆解

整个任务可以拆成四步：

• 插值：用已知点的规律估计缺失时段的值。这里用 interpolate.interp1d，它在相邻点之间按所选方式（线性、三次样条等）构造连续函数。
• 拟合：假设菌落量随时间服从指数增长，用 optimize.curve_fit 找出模型参数，并得到参数的协方差（不确定度）。
• 积分：对拟合出的连续模型求定积分，得到 0 到 24 小时的总增长量。

三步环环相扣：插值补全数据让拟合更稳定，拟合给出连续模型让积分有意义。

2.2 完整代码：一段能跑的最小程序

import numpy as np
from scipy import interpolate, optimize, integrate

# ---------- 1. 实验观测数据（t: 小时，N: 菌落量，含两个缺口） ----------
t_full = np.arange(0, 25, 1.0)
N_full = 100 * np.exp(0.05 * t_full) + 5 * np.random.default_rng(42).normal(size=len(t_full))
N_full[[7, 8]] = np.nan          # 第 7、8 小时数据缺失

# ---------- 2. 插值补全缺失值 ----------
mask = ~np.isnan(N_full)
interp = interpolate.interp1d(t_full[mask], N_full[mask], kind='cubic')
N_clean = N_full.copy()
N_clean[mask == False] = interp(t_full[mask == False])
print("补全后的第 7、8 小时估计值:", N_clean[[7, 8]].round(1))

# ---------- 3. 拟合指数模型 N = a * exp(b * t) ----------
model = lambda t, a, b: a * np.exp(b * t)
popt, pcov = optimize.curve_fit(model, t_full, N_clean, p0=[100, 0.05])
a, b = popt
print(f"拟合结果: a={a:.1f}, b={b:.4f}")
print(f"生长速率 b 的标准差: {np.sqrt(pcov[1, 1]):.4f}")

# ---------- 4. 数值积分求 24 小时总增长量 ----------
total, err = integrate.quad(model, 0, 24, args=(a, b))
print(f"24 小时总增长量: {total:.0f}（误差估计 {err:.1e}）")
​

这段代码的输出大致是：补全值落在 700~800 之间，拟合出的 b 接近真实值 0.05，总增长量约一万多——具体数字因随机种子而定，但数量级和趋势是稳定的。

2.3 每个环节的深意

为什么插值用三次样条而不是线性？线性插值在缺口两端会产生折角，拟合指数模型时这些不光滑点会拖累参数估计；三次样条保证一阶导数连续，拟合更稳。为什么拟合要有初值 p0？curve_fit 用局部优化算法（默认 Levenberg-Marquardt），从初值出发迭代，初值给得离谱会不收敛或收敛到局部最优。为什么积分用拟合后的模型而不是插值函数？拟合模型有明确的物理含义，积分结果才有解释意义；对插值函数积分只能得到"这段曲线下面积"，两者都能算，但回答的问题不同。

2.3 模型对不对：看残差，别只看拟合曲线

曲线贴着数据点走，不代表模型是对的。真正的检验是看残差——每个观测点与模型预测值的差。如果残差呈现随机散布，说明模型捕捉到了主要规律；如果残差有明显的弯曲或趋势，说明模型形式选错了。

# 拟合后计算残差并做简单诊断
residuals = N_clean - model(t_full, a, b)
print("残差均值（应接近 0）:", residuals.mean().round(3))
print("残差标准差（相对噪声水平）:", residuals.std().round(3))

# 检查残差是否随时间有趋势：分组看前段与后段的均值
half = len(residuals) // 2
print("前段残差均值:", residuals[:half].mean().round(3),
      "后段残差均值:", residuals[half:].mean().round(3))
​

如果前段与后段残差均值一正一负、符号相反，说明模型系统性偏差——比如真实增长不是纯指数，而是带饱和的逻辑斯蒂形式。这时候不要硬调参数，而是换模型形式重新拟合。残差分析是科学计算里最容易被跳过、却最有信息量的步骤，建议养成习惯：每次拟合完都看一眼残差。

2.4 给结论加上不确定性

科学报告不能只给"生长速率是 0.0498"，还要给"±多少"。curve_fit 返回的协方差矩阵正好用于此：对角线元素的平方根就是各参数的标准差，据此可以算出置信区间。

from scipy import stats

# 95% 置信区间：参数 ± t 值 × 标准差
se = np.sqrt(np.diag(pcov))
t_crit = stats.t.ppf(0.975, df=len(t_full) - 2)
ci_low, ci_high = b - t_crit * se[1], b + t_crit * se[1]
print(f"生长速率 b = {b:.4f}，95% 置信区间 [{ci_low:.4f}, {ci_high:.4f}]")
​

注意这里的自由度取"数据点数减参数个数"（两个参数，所以减 2）。有了置信区间，"这株细菌的繁殖速率约每小时增长百分之五，95% 置信区间从百分之四点几到百分之五点几"才是完整的科学结论。这也是为什么本教程反复强调：curve_fit 别只取 popt，把 pcov 一起存下来。

三、工程实践要点

3.1 程序的工程化组织

上面是最小程序，真实项目里建议这样组织：

• 数据与模型分离：数据读取、预处理写成一个函数，模型定义写成一个函数，拟合与输出再分开。
• 分步输出中间结果：每步打印关键数字（补全值、拟合参数、误差估计），方便定位问题出在哪个环节。
• 固定随机种子：示例里用了 default_rng(42)，保证结果可复现——科学计算里可复现性是底线。

把 2.2 的完整程序按函数重写一遍，就是这个样子：

def load_data():
    """生成或读取观测数据（此处为模拟）"""
    t = np.arange(0, 25, 1.0)
    y = 100 * np.exp(0.05 * t) + 5 * np.random.default_rng(42).normal(size=len(t))
    y[[7, 8]] = np.nan
    return t, y

def fill_missing(t, y):
    """三次样条插值补全缺失点"""
    mask = ~np.isnan(y)
    interp = interpolate.interp1d(t[mask], y[mask], kind='cubic')
    y_clean = y.copy()
    y_clean[~mask] = interp(t[~mask])
    return y_clean

def fit_exponential(t, y):
    """拟合指数模型并返回参数、协方差"""
    model = lambda t, a, b: a * np.exp(b * t)
    popt, pcov = optimize.curve_fit(model, t, y, p0=[100, 0.05])
    return model, popt, pcov

def main():
    t, y = load_data()
    y_clean = fill_missing(t, y)
    model, (a, b), pcov = fit_exponential(t, y_clean)
    total, err = integrate.quad(model, 0, 24, args=(a, b))
    print(f"b = {b:.4f} ± {np.sqrt(pcov[1, 1]):.4f}，总量 {total:.0f}")

if __name__ == "__main__":
    main()
​

这样每个环节都能单独测试、单独复用——换数据、换模型都不需要动其他函数。对一个小练习来说似乎小题大做，但养成这个习惯后，遇到几百行的大分析代码时你会感激自己。

3.2 常见失败模式对照

表里四条看着抽象，其实每条都能在五分钟内复现。我们把它们逐个触发一遍，以后见到同类报错就不会慌：

# 触发一：初值离谱，迭代直接飞掉
try:
    optimize.curve_fit(model, t_full, N_clean, p0=[1, 100])
except Exception as e:
    print("触发一，不收敛:", e)

# 触发二：参数比数据能支撑的多，方差爆炸
popt4, pcov4 = optimize.curve_fit(
    lambda t, a, b, c, d: a * np.exp(b * t) + c * t + d,
    t_full, N_clean)
print("触发二，参数标准差:", np.sqrt(np.diag(pcov4)).round(3))

# 触发三：插值函数在数据范围外求值
outside = interpolate.interp1d(t_full, N_clean, kind='cubic')
try:
    outside(25.5)
except ValueError as e:
    print("触发三，越界求值:", e)

# 触发四：被积函数带奇异点，静默失败
total_bad, err_bad = integrate.quad(lambda x: 1 / x, -1, 1)
print("触发四，奇异积分结果:", total_bad, err_bad)
​

触发一的教训是初值不能拍脑袋：指数项在 t 取 24 时已经放大到惊人的量级，初值差一个数量级，迭代就再也回不来。触发二更隐蔽——模型照样"拟合成功"，但参数标准差变成 inf 或天文数字，这是数据信息量撑不起参数个数的信号，解决办法是砍参数或加数据。触发三的教训是外推要显式声明：interp1d 默认不允许范围外求值，想控制行为就传 bounds_error 和 fill_value 两个参数。触发四最危险：奇异点不一定会抛异常，而是返回 inf 加一条警告——看到警告里出现 divergent 字样，先检查被积函数有没有分母为零的位置，再谈积分方法。这四条凑齐了，你在真实项目里遇到报错时，基本都能在表里对号入座。

⚠️ 常见坑：interp1d 默认不允许在数据范围外求值（会抛异常）。如果你的模型预测需要超出观测范围，要么显式处理，要么意识到外推本身不可靠——样条外推的结果通常没有物理意义。
💡 关键直觉：这条"插值 → 拟合 → 积分"的流水线是科学计算中最常见的组合拳。往后你会在信号处理、气象、金融、生物等几乎一切领域看到它的变体。跑通它，你就拥有了第一块完整的 SciPy 拼图。

3.3 扩展练习：把本章方法用在别处

下面的练习都不难，但能检验你是否真的理解了这条流水线：

1. 把指数模型换成逻辑斯蒂模型 N = K / (1 + exp(-r * (t - t0)))，重跑拟合与积分。你会看到拟合曲线在数据末端变得平缓——思考一下：当观测数据只覆盖增长早期时，逻辑斯蒂模型的参数还能被可靠估计吗？
2. 把插值方法从 cubic 换成 linear，重新拟合，比较两次的 b 和置信区间。缺口只有两个点时，方法选择的影响有多大？
3. 自己造一份带三个缺口的数据（随机挑三个时间点置为缺失），走完整流程。多缺口会不会让拟合的置信区间明显变宽？

这三个练习加起来半小时左右，做完后你对"数据质量如何影响结论"会有非常直观的认识——这是任何 API 文档都给不了的体验。

3.4 何时该找别的工具

这条流水线覆盖了很多问题，但它不是万能钥匙。数据规模、结论性质、数学形式越过某个边界后，硬用 SciPy 反而不划算——知道什么时候换工具，也是能力的一部分。

第一个边界是数据规模。示例只有 25 个点，插值拟合都是瞬时完成；当数据变成几百万行、还带分组和缺失标记时，pandas 的清洗与聚合语法比手写数组索引高效得多。正确姿势不是二选一，而是分工：pandas 管数据整理，把整理好的列转成数组交给 SciPy 算数值，再把结果接回表格。

第二个边界是结论的性质。curve_fit 给出的协方差适合简单场景，但正式的回归分析——多个解释变量、分类变量、异方差检验、模型选择——应该交给 statsmodels 这类专门做统计建模的工具，它输出的结果本身就带完整的推断信息。我们手动算置信区间的思路没有错，只是问题变大后，别自己手搓标准误。

第三个边界是数学形式。如果问题要的是符号推导，比如"这个积分的原函数长什么样"，数值积分永远给不了答案，这时该找 sympy 做符号计算。SciPy 给数字，sympy 给公式，两者是互补关系而不是竞争关系。

判断标准其实很朴素：当"数据整理"的篇幅开始超过"数值计算"时，引入表格工具；当你要的是公式而非数字时，引入符号工具。SciPy 的定位始终是数值算法库，认识它的边界，和认识它的能力一样重要。

本节速览

• 要点一：插值补缺用三次样条比线性更平滑，利于后续拟合。
• 要点二：curve_fit 需要合理初值 p0，否则可能不收敛或陷入局部最优。
• 要点三：拟合结果中的协方差矩阵给出参数的不确定度，是科学报告必需的。
• 要点四：对拟合模型积分才能得到有物理含义的总量。
• 要点五：程序组织上做到数据、模型、流程分离，分步打印中间结果。
• 要点六：固定随机种子保证可复现，是科学计算的底线习惯。

至此，你已经用 SciPy 走完一个完整的小项目。第 2 章我们深入这三件武器本身：积分模块的自适应算法、优化模块的各类求解器、插值的样条家族，从"会用"走向"会选、会调"。