如何使用scipy将两组带共享参数的实验数据拟合至微分方程（已解决）

云上老周

2026-9-1

问题背景

我需要同时拟合两组存在共享参数的实验数据，数据对应酶促反应Gl→Gm伴随HP衰减的化学反应过程，期望得到参考的拟合效果。
此前我已经通过symfit包完成数据拟合，但为了后续开展蒙特卡洛模拟等数据处理工作，需要改用scipy/numpy实现拟合逻辑。

我尝试编写的scipy代码如下：

import matplotlib.pyplot as plt
import numpy as np
import scipy as sp

# 从预处理后的txt文件读取数据集
with open("ydata.txt", "r") as csv_file:
    ydata = np.loadtxt(csv_file, delimiter = ',')

with open("ydata2.txt", "r") as csv_file:
    ydata2 = np.loadtxt(csv_file, delimiter = ',')

xdata = np.arange(0, len(ydata))
fulldata = np.column_stack([ydata,ydata2])


# 定义Gl -> Gm酶促反应伴随HP衰减的方程
def f(C, t, k, a, b):
    GL = ydata
    GM = ydata2
    
    dGLdt = -k*GL - GL/a
    dGMdt = k*GL - GM/b
    
    return [dGLdt, dGMdt] 

guess = (1e-3, 10, 10,1 )

popt, pcov = sp.optimize.curve_fit(f, xdata, fulldata, guess)

运行代码时出现如下报错：

File "/Users/karensantos/Desktop/Codes/Stack_question.py", line 52, in <module>
    popt, pcov = sp.optimize.curve_fit(f, xdata, fulldata, guess)

  File "/opt/anaconda3/lib/python3.8/site-packages/scipy/optimize/minpack.py", line 784, in curve_fit
    res = leastsq(func, p0, Dfun=jac, full_output=1, **kwargs)

  File "/opt/anaconda3/lib/python3.8/site-packages/scipy/optimize/minpack.py", line 410, in leastsq
    shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)

  File "/opt/anaconda3/lib/python3.8/site-packages/scipy/optimize/minpack.py", line 24, in _check_func
    res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))

  File "/opt/anaconda3/lib/python3.8/site-packages/scipy/optimize/minpack.py", line 484, in func_wrapped
    return func(xdata, *params) - ydata

ValueError: operands could not be broadcast together with shapes (2,98) (98,2)

目前单独拟合其中一个方程可以正常运行，但两组数据存在共享参数k，且产物GM浓度依赖底物GL浓度，必须同时拟合两组数据才能得到准确参数，需要明确使用scipy优化模块完成两组实验数据联合拟合的实现方法。

解决方案

你的代码存在三个核心错误：

定义的f只是微分方程组的导数形式，curve_fit要求模型函数输入自变量、参数后，直接返回和观测数据形状一致的预测值，不能直接返回微分表达式。
输出形状不匹配：观测数据fulldata形状为(时间点数, 2)，原函数返回值形状为(2, 时间点数)，触发广播报错。
微分方程定义逻辑错误：直接把全局观测值代入导数计算，ODE求解需要从初始浓度出发递推所有时间点的浓度，不能在递推步骤中直接代入观测值，同时初始参数猜测漏了反应初始浓度项。

联合拟合带共享参数的常微分方程数据，正确实现步骤为：

用scipy.integrate.odeint求解给定参数下的微分方程数值解，得到对应时间点的GL、GM预测值
编写适配curve_fit的模型函数，保证输出形状和观测数据完全一致
传入拼接后的观测数据做拟合，参数k、a、b会在拟合过程中同时约束两组数据的残差，实现共享参数的联合估计

修正后的可运行代码：

import numpy as np
import scipy as sp
from scipy.integrate import odeint

# 读取数据
with open("ydata.txt", "r") as csv_file:
    ydata = np.loadtxt(csv_file, delimiter = ',')
with open("ydata2.txt", "r") as csv_file:
    ydata2 = np.loadtxt(csv_file, delimiter = ',')

xdata = np.arange(0, len(ydata))
fulldata = np.column_stack([ydata, ydata2])  # 形状 (n_t, 2)，列对应GL、GM观测值

# 定义微分方程组，传入初始状态C=[GL0, GM0]和待拟合参数
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适配的模型函数：输入自变量x、待拟合参数，返回和fulldata同形状的预测值
def model(x, k, a, b, GL0, GM0):
    # 求解ODE，得到每个时间点的GL、GM预测值，形状(n_t, 2)
    pred = odeint(ode_system, [GL0, GM0], x, args=(k, a, b))
    return pred

# 初始参数猜测：k, a, b, GL初始值, GM初始值
guess = (1e-3, 10, 10, ydata[0], ydata2[0])

# 联合拟合
popt, pcov = sp.optimize.curve_fit(model, xdata, fulldata, p0=guess)

# 提取拟合结果
k_fit, a_fit, b_fit, GL0_fit, GM0_fit = popt
print(f"拟合参数：k={k_fit:.4e}, a={a_fit:.2f}, b={b_fit:.2f}, GL0={GL0_fit:.2f}, GM0={GM0_fit:.2f}")

# 生成拟合曲线用于后续绘图或分析
y_pred = model(xdata, *popt)

关键说明

初始浓度GL0、GM0如果实验测量存在误差，可以作为待拟合参数和k、a、b一起估计，这是动力学拟合的常规操作；如果初始浓度是准确已知的固定值，可以用functools.partial把固定参数传给model函数，不用放到待拟合参数列表里。
输出形状必须严格和输入的观测数据对齐：观测数据是(时间点数, 2)，模型返回值也必须是这个形状，否则会触发广播类报错。
后续开展蒙特卡洛模拟时，直接在拟合得到的最优参数基础上，给观测值加对应噪声水平的随机扰动，重复上述拟合流程即可，拟合返回的参数协方差矩阵pcov也可以直接用来生成参数的模拟分布。

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

如何使用scipy将两组带共享参数的实验数据拟合至微分方程（已解决）

关键说明 #

关键说明