1

I would like to know how to extract the intial values for the state variables.

Basically, the intial values are also considered as parameters.

I'm providing the initial values because they are needed for integration. However, these values are sometimes considered as additional parameters (next to the model parameters), so that the initial values provided are also considered are initial guesses. In fact, when doing growth experiments, we have some initial conditions, but we may not know all of them (depending on the specific experimental conditions and system under investigation).

Consider a simple microbial growth model with the growth rate mu governed by the well-known Monod kinetic (parameters mumax and ks) and a constant biomass-substrate yield (parameter yxs) and a growth-associated product formation (parameter ypx). See the code below.

For further implementations I would like to calculate first-order sensitivities over time, parameter correlations, etc.

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)
Ethan Field
  • 4,646
  • 3
  • 22
  • 43
kossikater
  • 41
  • 5
  • Could you add an example? I'd love to help you with this (and thanks for using symfit) but I'm not sure what you need. You already provide the initial values to the `ODEModel` in order for it to be able to integrate, right? – tBuLi May 10 '18 at 17:43

2 Answers2

1

Great question. The short answer is no, because this has not been implemented yet. (Though it has been on the list for a while, see Issue #53.)

However, in some cases there might already be a workaround. Lets say you want

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

in your example above. symfit will currently only check model_dict for Parameter objects, so it will not see x0. So if you model explicitly depends on x0 it will work. Perhaps you can find some way to make x0 appear in your model. Worst case scenario, you could multiply one of the components with some trivial identity like 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,
}

However, this is provided without warranty ;). Something I noticed in my own research when I had a model who's model_dict did explicitly depend on the initial parameter was that symfit had difficulty estimating the Jacobian for such parameters, so fitting was unstable. However, this was a couple of versions ago and many things have changed since then so give it a try!

You can also use a non-gradient based method such as Nelder-Mead or hopefully very soon Differential Evolution.

I will look into this, I opened an issue for this here. If you'd like to get involved in this any help would definitely be welcome, for example when it comes to providing a minimal working example that I could use as a test.

tBuLi
  • 2,295
  • 2
  • 16
  • 16
  • Thx for the workaround. For sure I can provide a minimal working example for testing. Just pm me. – kossikater May 22 '18 at 06:23
  • I had also some thoughts about initial values as parameters to optimize: From my experience, we do have quite often some model states that we can't (or don't) measure, so at least for these states we need the initial values to see as parameters, although we are primarily interested in the model parameters. Also, I think one has always to carefully evaluate parameter estimates, e.g. wrt correlations, etc. So there will always be some tradeoff between usability (important for acceptance) and complexity of output. But that's some kind of philosopical question... Your work is highly appreciated! – kossikater May 22 '18 at 06:44
  • Unfortunately I can't PM you on stackoverflow. Perhaps you can reply to this issue and provide a minimal example? Your help is much appriciated :). https://github.com/tBuLi/symfit/issues/160 – tBuLi Jun 05 '18 at 11:05
1

I can confirm that this suggestion works in Symfit 0.5.1, with the caveat that the "initial" value needs to be specified as, e.g., x0*.value*. This creates an odd issue when you call "print(fit_result)", it complains that "'NoneType' object has no attribute 'sqrt'", but everything else seems fine.

Here's a (more or less) minimal example to demonstrate this:

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

JinjerJohn
  • 403
  • 4
  • 6
  • I think the print(fit_result) error happens when the variance is undefined - perhaps too much error in the estimate. I solved this in the above example by defining initial s as "10**s0.value" and changing the s0 bounds to their respective log10s. This is probably better in any case, for a variable that should really be hopping around in log-space rather than linear space. – JinjerJohn Nov 04 '19 at 13:35
  • Also, it's worth noting that _truly_ trivial identities like `s0-s0+1` don't work, sympy's too clever and just ignores them. Looking forward to incorporation of [this update](https://github.com/tBuLi/symfit/pull/279)! – JinjerJohn Nov 04 '19 at 13:42