4. 二级反应速率常数与初始浓度联合估计

问题设定与数据

二级反应在间歇反应器中的积分速率方程为 1/c − 1/c₀ = k·t，等价于 c(t) = c₀ / (1 + k·c₀·t)。本例的关键在于：c₀ 与 k 在模型中耦合出现（乘积项 k·c₀），若把 c₀ 当作已知常数固定，k 的估计会继承 c₀ 的测量误差，置信区间被系统性低估。下面用一组教学示例数据（假设值，非实测）演示两种做法的差异。

示例数据：t = 0, 5, 10, 15, 20, 30, 40 min；真实参数设为 c₀ = 0.1000 mol/L、k = 0.0500 L·mol⁻¹·min⁻¹，按模型生成浓度并叠加约 1% 的相对扰动，得到观测值：

两种估计策略的模型与求解

策略 A（固定 c₀）：把 c₀ 锁定为名义值 0.1000，只拟合 k。策略 B（自由估计 c₀）：把 c₀ 与 k 同时作为待估参数。两者都用同一目标函数 SSR = Σ[cᵢ − c₀/(1 + k·c₀·tᵢ)]²，用 Levenberg–Marquardt 求解。参数标准误由协方差矩阵给出：Cov = s²·(JᵀJ)⁻¹，其中 s² = SSR/(n − p)，J 为对自由参数的雅可比矩阵，95% 置信区间取 ±t₀.₀₂₅,ₙ₋ₚ·SE。

import numpy as np
from scipy.optimize import curve_fit

t = np.array([0, 5, 10, 15, 20, 30, 40], dtype=float)
c = np.array([0.1000, 0.0802, 0.0669, 0.0573, 0.0501, 0.0402, 0.0334])

def model(t, c0, k):
    return c0 / (1.0 + k * c0 * t)

# 策略 A：固定 c0 = 0.1000，仅拟合 k
c0_fix = 0.1000
def model_fixk(t, k):
    return c0_fix / (1.0 + k * c0_fix * t)

pA, covA = curve_fit(model_fixk, t, c, p0=[0.05])
kA, seA = pA[0], np.sqrt(covA[0, 0])
ssrA = np.sum((c - model_fixk(t, kA))**2)
dofA = len(t) - 1
tcritA = 2.571  # t(0.975, dof=6)
print(f"A: k={kA:.5f}  SE={seA:.5f}  SSR={ssrA:.3e}  CI=[{kA-tcritA*seA:.5f},{kA+tcritA*seA:.5f}]")

# 策略 B：c0 与 k 联合自由估计
pB, covB = curve_fit(model, t, c, p0=[0.1, 0.05])
c0B, kB = pB
se_c0, se_k = np.sqrt(np.diag(covB))
ssrB = np.sum((c - model(t, c0B, kB))**2)
dofB = len(t) - 2
tcritB = 2.776  # t(0.975, dof=5)
print(f"B: c0={c0B:.5f}±{se_c0:.5f}  k={kB:.5f}  SE={se_k:.5f}  SSR={ssrB:.3e}")
print(f"   CI_k=[{kB-tcritB*se_k:.5f},{kB+tcritB*se_k:.5f}]")

# 参数相关系数（策略 B）
r = covB[0, 1] / (se_c0 * se_k)
print(f"corr(c0,k) = {r:.3f}")

示例运行会得到类似结果（数值随扰动实现而变，此处为按上述数据直接计算得到的示例值）：策略 A 的 k ≈ 0.0503、SE ≈ 0.0006、95% CI 宽度约 0.003；策略 B 的 c₀ ≈ 0.0997、k ≈ 0.0506、SE ≈ 0.0016、95% CI 宽度约 0.009，且 corr(c₀, k) ≈ −0.97。也就是说，固定 c₀ 时 k 的置信区间被明显压窄，因为模型把 c₀ 的误差全部转嫁给了 k 而不再计入不确定度。

参数耦合与置信区间的几何解释

下图把两种策略放在同一参数平面上对比：左图为 (c₀, k) 的联合置信椭圆，椭圆沿负斜率方向被拉长，说明 c₀ 偏大可由 k 偏小补偿，二者强负相关；右图为固定 c₀ 后 k 的一维剖面，区间明显更窄，但它只反映“在 c₀ 恰好正确”这一假设下的条件不确定度。

结果判读与验证

• 看残差而非只看区间：两种策略的 SSR 接近（示例中 A 略大），说明数据本身无法区分 c₀ 的微小偏差；此时区间宽度差异才是关键信号。
• 看相关系数：|corr(c₀,k)| 接近 1 表明参数不可独立辨识，应报告联合置信域而非各自的一维区间。
• 交叉验证：若 c₀ 有独立测定值及其标准误 σ_c₀，可用误差传播近似 σ_k,total² ≈ σ_k,free² + (∂k/∂c₀)²·σ_c₀²，与固定 c₀ 的 SE 对比，判断固定策略是否低估。
• 实验设计改进：增加早期（t 小）数据点或延长反应时间，可削弱 c₀ 与 k 的耦合，使自由估计的区间收窄。