1

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

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

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

考虑一个简单的微生物生长模型，其生长速率 mu 由著名的 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)
4

2 回答 2

1

我可以确认这个建议在 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)

具有拟合曲线的数据点

于 2019-11-04T01:19:36.977 回答
1

好问题。简短的回答是否定的，因为这还没有实施。（虽然它已经在列表上一段时间了，请参阅问题 #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或希望很快就会使用差分进化。

我会调查一下，我在这里打开了一个问题。如果您想参与其中，任何帮助肯定都会受到欢迎，例如在提供我可以用作测试的最小工作示例时。

于 2018-05-21T10:49:35.107 回答