python - 使用 Symfit for Python 进行 ODE 拟合 : How to get estimations for intial values?

Question

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

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

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

考虑一个简单的微生物生长模型，其生长速率 mu 由众所周知的莫诺动力学(参数 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 。如果您想参与其中，任何帮助肯定都会受到欢迎，例如，在提供一个我可以用作测试的最小工作示例时。