Make the SIR model less weird

[ben18785]

At the moment, it's a really bizarre model because:

• It only outputs infected and recovered populations (even though we don't observe the susceptible population, we should still output it). I think since most people will expect SIR to be a three output model.
• A parameter is the initial susceptible population size (which is fine), but this overwrites whatever the value of S0 the model is instantiated with. This could be solved by making the y0 instantiation vector only two dimensional.

[MichaelClerx]

It only outputs infected and recovered populations (even though we don't observe the susceptible population, we should still output it). I think since most people will expect SIR to be a three output model.

Nooooooo! I want it to be realistic! I've seen too many papers where people fit to unobservable quantities in toy models. If you really want you can make some kind of switch to output all three, but there's no way of outputting all 3 without making the fitting problem unrealistic

[MichaelClerx]

Example from the Beeler-Reuter model:

def simulate(self, parameters, times):
    """ See :meth:`pints.ForwardModel.simulate()`. """
    y0 = [self._v0,
          self._cai0,
          self._m0,
          self._h0,
          self._j0,
          self._d0,
          self._f0,
          self._x10]

    solved_states = scipy.integrate.odeint(
        self._rhs, y0, times, args=(parameters,), hmax=self._I_Stim_length,
        rtol=self._rtol, atol=self._atol)

    # Only return the observable (V, Cai)
    return solved_states[:, 0:2]

def simulate_all_states(self, parameters, times):
    """
    Returns all state variables that ``simulate()`` does not return.
    """
    y0 = [self._v0,
          self._cai0,
          self._m0,
          self._h0,
          self._j0,
          self._d0,
          self._f0,
          self._x10]

    solved_states = scipy.integrate.odeint(
        self._rhs, y0, times, args=(parameters,), hmax=self._I_Stim_length)

    # Return all states
    return solved_states

[ben18785]

Cool. Thanks. To be clear, I like the example — real data and hence more persuasive. The only things I was wondering about were: a) don’t we have a way of fitting to only a subset of outputs via multi-output likelihoods? and b) change the way the S0 parameter is set as currently the parameters overwrite the instantiation values.

On 26 Sep 2019, at 11:00, Michael Clerx [email protected] wrote:

Example from the Beeler-Reuter model:

def simulate(self, parameters, times): """ See :meth:pints.ForwardModel.simulate(). """ y0 = [self._v0, self._cai0, self._m0, self._h0, self._j0, self._d0, self._f0, self._x10]

solved_states = scipy.integrate.odeint(
    self._rhs, y0, times, args=(parameters,), hmax=self._I_Stim_length,
    rtol=self._rtol, atol=self._atol)

# Only return the observable (V, Cai)
return solved_states[:, 0:2]

def simulate_all_states(self, parameters, times): """ Returns all state variables that simulate() does not return. """ y0 = [self._v0, self._cai0, self._m0, self._h0, self._j0, self._d0, self._f0, self._x10]

solved_states = scipy.integrate.odeint(
    self._rhs, y0, times, args=(parameters,), hmax=self._I_Stim_length)

# Return all states
return solved_states

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub, or mute the thread.

[MichaelClerx]

don’t we have a way of fitting to only a subset of outputs via multi-output likelihoods

I can't really think of a real-life situation where you'd do that, other than testing/benchmarking

Happy for change (b) though! Or for adding a 2nd simulation method like we have for BR