ODE拟合与Symfit for Python:如何获得初始值的估计？

Question

我想知道如何提取状态变量的初始值。

基本上，初始值也被认为是参数。

我提供初始值，因为它们是集成所需的。但是，这些值有时被视为附加参数(模型参数旁边)，因此提供的初始值也被视为初始猜测。事实上，在做生长实验时，我们有一些初始条件，但我们可能不知道所有的初始条件(取决于具体的实验条件和正在研究的系统)。

考虑一个简单的微生物生长模型，其生长速率µ由众所周知的Monod动力学(参数mumax和ks)、恒定的生物量-底物产量(参数yxs)和与生长相关的产物形成(参数ypx)控制。请参阅下面的代码。

对于进一步的实现，我想计算随时间变化的一阶灵敏度、参数相关性等。

from symfit import variables, parameters, ODEModel, Fit, D
import numpy as np
from matplotlib import pyplot as plt
# define ODE model
ks, mumax, ypx, yxs = parameters('ks, mumax, ypx, yxs')
p, s, x, t = variables('p, s, x, t')
model_dict = {
    D(x, t): mumax * s / (ks + s) * x,
    D(s, t): -1/yxs * mumax * s / (ks + s) * x,
    D(p, t): ypx * mumax * s / (ks + s) * x,
}
ode_model = ODEModel(model_dict, initial={t:0.0, x:0.1, s:20.0, p:0.0})

# Generate noisy measurement data
tSample = np.array([1, 2, 4, 6, 8, 12])
pSample, sSample, xSample = ode_model(t=tSample, ks=0.01, mumax=0.4, ypx=0.2, yxs=0.5)
xRelErr = 0.05
sRelErr = 0.10
pRelErr = 0.10
xNoise = xSample + xSample * xRelErr * np.random.randn(xSample.size)
sNoise = sSample + sSample * sRelErr * np.random.randn(sSample.size)
pNoise = pSample + pSample * pRelErr * np.random.randn(pSample.size)

# constraints for parameters
ks.min = 0.0
mumax.min = 0.0
ypx.min = 0.0
yxs.min = 0.0

# initial guesses for parameters
ks.value = 0.01
mumax.value = 0.4
ypx.value = 0.2
yxs.value = 0.5

# perform fitting
fit = Fit(ode_model, t=tSample, x=xNoise, s=sNoise, p=pNoise)
fit_result = fit.execute()

# show fit
print(fit_result)
print(fit_result.covariance_matrix)

Answer 1

问得好。简短的答案是否定的，因为这还没有实现。(尽管它出现在列表中已经有一段时间了，请参见Issue #53。)

但是，在某些情况下，可能已经有了解决方法。假设你想要

x0 = Parameter('x0')
ode_model = ODEModel(model_dict, initial={t: 0.0, x: x0, s: 20.0, p: 0.0})

在上面的例子中。symfit当前仅检查model_dict中的Parameter对象，因此它不会看到x0。因此，如果您显式地建模依赖于x0，那么它将会工作。也许你可以找到一些方法让x0出现在你的模型中。最坏的情况是，您可以将其中一个组件与一些简单的标识(如cos(x0)**2 + sin(x0)**2 )相乘

model_dict = {
    D(x, t): mumax * s / (ks + s) * x * (cos(x0)**2 + sin(x0)**2),
    D(s, t): -1/yxs * mumax * s / (ks + s) * x,
    D(p, t): ypx * mumax * s / (ks + s) * x,
}

但是，这是在没有保证的情况下提供的;)在我自己的研究中，当我有一个模型的model_dict确实明确地依赖于初始参数时，我注意到symfit很难估计这些参数的雅可比，所以拟合是不稳定的。然而，这是几个版本之前的事情，从那时起，许多事情都发生了变化，所以试一试吧！

您也可以使用基于非梯度的方法，如Nelder-Mead或希望很快使用Differential Evolution。

我会调查这一点，我为这个here打开了一个问题。如果你想参与其中，任何帮助都是绝对受欢迎的，例如，当涉及到提供一个我可以用作测试的最小工作示例时。

Answer 2

我可以确认this建议在Symfit 0.5.1中有效，但需要注意的是，“初始”值需要指定为，例如，x0*.value*。这就产生了一个奇怪的问题，当你调用"print(fit_result)“时，它会抱怨"'NoneType‘对象没有’sqrt‘属性”，但其他一切似乎都很好。

下面是一个(或多或少)演示这一点的最小示例：

from symfit.core.minimizers import DifferentialEvolution, BFGS, BasinHopping
import symfit as sf
import numpy as np
import matplotlib.pyplot as plt

signal = np.array([ 3.69798582e-04, -4.42968073e-04, -1.62901574e-04, -2.66074317e-03,
       -2.69579455e-03,  7.74354388e-04,  2.19358448e-03,  6.38607419e-03,
        3.78642692e-03, -3.27400548e-03, -4.34699571e-03,  2.41753615e-04,
        9.96158926e-03,  1.89990480e-02,  3.60202254e-02,  7.63199623e-02,
        1.43007912e-01,  2.48988030e-01,  3.95977066e-01,  5.65958061e-01,
        7.12988678e-01,  8.09083671e-01,  8.50628154e-01,  8.72826358e-01,
        8.83509851e-01,  8.81799211e-01,  8.81649458e-01,  8.88619439e-01,
        8.82045565e-01,  8.91226071e-01,  9.07541872e-01,  9.12674298e-01,
        9.06886981e-01,  8.98514353e-01,  9.02275457e-01,  9.07204574e-01,
        9.01857725e-01,  8.92016771e-01,  8.99608399e-01,  8.98652789e-01,
        8.98421894e-01,  8.94907751e-01,  9.09122059e-01,  9.13065656e-01,
        9.09254937e-01,  9.10644851e-01,  8.99605570e-01,  8.99581621e-01,
        9.12054065e-01,  9.06361431e-01,  9.03838466e-01,  9.09130042e-01,
        9.13443979e-01,  9.18176457e-01,  9.08160892e-01,  9.03154574e-01,
        9.03069088e-01,  9.02597396e-01,  8.98854582e-01,  9.07801262e-01])
cycles = np.arange(len(signal))

s,c = sf.variables('s,c')
r,K,s0 = sf.parameters('r,K,s0')
lg2_e = np.log2(np.exp(1))

s0.min = 1e-7
s0.value = 1e-6
s0.max = 1e-5
K.min = 0.1
K.value = 1
K.max = 2
r.min = 0.1/lg2_e
r.value = 1/lg2_e
r.max = 2/lg2_e


logistic_eqn = {
    sf.D(s,c): r*s*(1-s/K) * (sf.cos(s0)**2 + sf.sin(s0)**2)
}

logistic_model = sf.ODEModel(logistic_eqn, initial = {c:0.0,s:s0.value})

logistic_fit = sf.Fit(logistic_model,c=cycles,s=signal)#, minimizer=DifferentialEvolution)
fit_result = logistic_fit.execute()
#print(fit_result)
sig_fit, = logistic_model(c=cycles, **fit_result.params)
plt.plot(cycles, signal, 'o')
plt.plot(cycles, sig_fit)

Data points with fitted curve