如何将2个实验数据拟合成带枕的导数方程？

Question

我试图同时拟合2个实验数据，因为它有一些共享的参数。这是一个化学反应，我希望得到拟合，如所附图片所示。我已经成功地使用symfit包来拟合我的数据，但是，为了进一步处理数据(用monte模拟)，我需要使用scipy/numpy --我尝试过的代码是：

GL转换为GM及拟合

供下载的数据集

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

# Open dataset from txt file after extraction from brute data:
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])


# Define the equation considering the enzymatic reaction Gl -> Gm with the HP decay.
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) 

我可以用curve_fit一次求解一个方程，但我需要在一起找到所有正确的共享参数(k)，因为GM依赖GL (产品和衬底)。

如何使用枕优化来拟合这两个实验数据？

提前谢谢你，

Answer 1

您可以将数组连接在一维数组中，以便与curve_fit一起运行。

我不能举你的例子，所以我会做一个例子

import numpy as np
from scipy.optimize import curve_fit
def cost(x, a, b):
    return np.hstack(f(a, b, x))
def f(a,b,x):
    return a * x**3, a**2*np.exp(-(x-b/a)**2/a)
x = np.linspace(-2, 2)
y1, y2 = f(4.5, 2.3, x)
initial_guess = (1,1)
params, _ = curve_fit(cost, x, np.hstack([y1, y2]),initial_guess)
print(params)

在这个例子中，函数f接受两个参数和x数据，我使用它来计算(y1，y2)，然后使用curve_fit来确定生成(y1，y2)的参数。

编辑1

使用OP提供的数据，合适的情况可能是

def f(params, xdata, ydata, ydata2):
    C = xdata
    t, k, a, b = params
    GL = ydata
    GM = ydata2
    
    dGLdt = -k*GL - GL/a
    dGMdt = k*GL - GM/b
    
    return np.hstack([dGLdt, dGMdt])

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

popt, pcov = scipy.optimize.leastsq(f, guess, args=(xdata, ydata, ydata2))

参数的选择是[ 1.00000000e-03, -1.71943255e-69, 1.60693865e+61, 1.60694078e+60]，这是一个平凡的解决方案，(a，b)增长到无穷大，k变为零。t不改变，因为目标函数不随t而变化。

我觉得你应该重新考虑你的模型方程。