使用scipy拟合共享参数的两组化学反应ODE数据集求解

开发者阿泽

2026-9-3

问题说明

长期潜水，首次发帖求助。
当前研究的化学体系仅能在特定时间段内被检测，需要同时表征反应进程与信号衰减过程，对应的微分方程如下：

Derivative(GL, t): (-k*GL) - GL/a,
Derivative(GM, t): (k*GL) - GM/b,

此前已经使用symfit包完成了数据拟合（系统拟合效果见下图），但后续需要开展Monte Carlo模拟，需改用scipy实现相同拟合功能。
Chemical reaction and fitting using symfit

最初尝试按如下方式定义方程：

def f(C, xdata):
    GL = ydataScaled
    GM = ydataScaled2
    
    dGLdt = -k*GL - GL/a
    dGMdt = k*GL - GM/b
    
    return [dGLdt, dGMdt]

但调用optimize.minimize或odeint都无法完成拟合，需要scipy环境下针对存在共享参数的双y数据集的正确拟合实现方案。

运行代码时抛出如下错误：

runfile('/Users/karensantos/Desktop/Codes/Stack_question.py', wdir='/Users/karensantos/Desktop/Codes')
2
(512, 32768)
float64
/opt/anaconda3/lib/python3.8/site-packages/nmrglue/fileio/convert.py:68: UserWarning: Incompatible dtypes, conversion not recommended
  warn("Incompatible dtypes, conversion not recommended")
Traceback (most recent call last):

  File "/Users/karensantos/Desktop/Codes/Stack_question.py", line 112, in <module>
    popt, pcov = sp.optimize.minimize(f, xdata, args = (ydataScaled, ydataScaled2))

  File "/opt/anaconda3/lib/python3.8/site-packages/scipy/optimize/_minimize.py", line 612, in minimize
    return _minimize_bfgs(fun, x0, args, jac, callback, **options)

  File "/opt/anaconda3/lib/python3.8/site-packages/scipy/optimize/optimize.py", line 1101, in _minimize_bfgs
    sf = _prepare_scalar_function(fun, x0, jac, args=args, epsilon=eps,

  File "/opt/anaconda3/lib/python3.8/site-packages/scipy/optimize/optimize.py", line 261, in _prepare_scalar_function
    sf = ScalarFunction(fun, x0, args, grad, hess,

  File "/opt/anaconda3/lib/python3.8/site-packages/scipy/optimize/_differentiable_functions.py", line 76, in __init__
    self._update_fun()

  File "/opt/anaconda3/lib/python3.8/site-packages/scipy/optimize/_differentiable_functions.py", line 166, in _update_fun
    self._update_fun_impl()

  File "/opt/anaconda3/lib/python3.8/site-packages/scipy/optimize/_differentiable_functions.py", line 73, in update_fun
    self.f = fun_wrapped(self.x)

  File "/opt/anaconda3/lib/python3.8/site-packages/scipy/optimize/_differentiable_functions.py", line 70, in fun_wrapped
    return fun(x, *args)

TypeError: f() takes 2 positional arguments but 3 were given

解决方案

代码存在三个核心问题，逐一修正即可正常拟合：

ODE定义逻辑错误：微分方程求解时，状态变量GL、GM是求解过程中动态计算的值，不能直接把观测到的ydataScaled、ydataScaled2赋值给状态变量，观测值仅用于计算拟合残差。
函数传参不符合接口规范：odeint要求传入的ODE函数参数顺序为(状态变量数组, 时间序列, 待拟合参数)，原函数的参数顺序、参数个数都不匹配，直接触发了参数数量报错。
优化器调用逻辑错误：minimize是标量优化函数，需要手动传入返回单值损失的目标函数；拟合双输出ODE更简便的方式是配合curve_fit实现，只需把两个通道的预测值拼接为一维数组，和拼接后的观测值一一对应即可。

注意：用scipy拟合时不需要保留symfit相关的变量、参数定义代码，删掉对应行即可。

替换原代码中拟合部分的内容为如下代码即可（前面NMR数据读入、积分处理的部分不需要改动）：

from scipy.integrate import odeint
from scipy.optimize import curve_fit
import numpy as np

# 定义ODE系统
def ode_system(C, t, k, a, b):
    GL, GM = C
    dGLdt = -k * GL - GL / a
    dGMdt = k * GL - GM / b
    return [dGLdt, dGMdt]

# 定义curve_fit调用的拟合函数，返回拼接后的一维预测值
def fit_func(t, k, a, b, GL0, GM0):
    C_init = [GL0, GM0]
    sol = odeint(ode_system, C_init, t, args=(k, a, b))
    return np.hstack([sol[:, 0], sol[:, 1]])

# 拼接两个通道的观测值为一维数组
y_merged = np.hstack([ydataScaled, ydataScaled2])
# 待拟合参数初始值，根据数据的实际时间尺度调整即可
init_params = [0.01, 10, 10, 1.0, 0.0]

# 执行拟合
popt, pcov = curve_fit(fit_func, xdata, y_merged, p0=init_params)
k_fit, a_fit, b_fit, GL0_fit, GM0_fit = popt

# 生成拟合曲线用于后续绘图或分析
fit_result = odeint(ode_system, [GL0_fit, GM0_fit], xdata, args=(k_fit, a_fit, b_fit))
GL_curve = fit_result[:, 0]
GM_curve = fit_result[:, 1]

补充说明：

初始参数init_params不要随意设置，否则容易出现拟合不收敛的问题，可以先手动调整参数让预测曲线大致贴合观测值，再作为初始值传入。
如果两个通道的噪声水平差异较大，可以给curve_fit传入sigma参数做加权拟合，结果稳定性会更好。
如果需要用minimize实现，只需自行编写损失函数：输入待拟合参数，调用ODE求解得到预测值，返回预测值和观测值的均方误差即可，核心逻辑和上述代码一致。

内容的提问来源于stack exchange，提问作者Karen DOS SANTOS