Defining observed model states and observation process in model declaration

[twallema]

Current implementation

# Import the ODEModel class
from pySODM.models.base import ODEModel

# Define the model equations
class ODE_SIR(ODEModel):
    """
    Simple SIR model without dimensions
    """
    
    state_names = ['S','I','R']
    parameter_names = ['beta','gamma']

    @staticmethod
    def integrate(t, S, I, R, beta, gamma):
        
        # Calculate total population
        N = S+I+R
        # Calculate differentials
        dS = -beta*S*I/N
        dI = beta*S*I/N - 1/gamma*I
        dR = 1/gamma*I

        return dS, dI, dR

# Define parameters and initial condition
params={'beta':0.35, 'gamma':5}
init_states = {'S': 1000, 'I': 1}

# Initialize model
model = ODE_SIR(states=init_states, parameters=params)

# Define dataset
data=[d, ]
states = ["I",]
log_likelihood_fnc = [ll_negative_binomial,]
log_likelihood_fnc_args = [alpha,]

# Calibated parameters and bounds
pars = ['beta',]
labels = ['$\\beta$',]
bounds = [(1e-6,1),]

# Setup objective function (no priors --> uniform priors based on bounds)
objective_function = log_posterior_probability(model, pars, bounds, data, states, log_likelihood_fnc, log_likelihood_fnc_args, labels=labels)

Proposed implementation

Define what states are observed, and what statistical process underpins the observations in the model declaration. Then provide the parameters of the observation process (if any) upon model declaration. When simulating the model, automatically add the observation noise to the observed states, thus avoiding the use of add_gaussian_noise() or add_poisson_noise(). This further simplifies the call to log_posterior_probability.

# Import the ODEModel class
from pySODM.models.base import ODEModel

# Define the model equations
class ODE_SIR(ODEModel):
    """
    Simple SIR model without dimensions
    """
    
    state_names = ['S','I','R']
    observed_state_names = ['I']
    observation_process = ['negative_binomial']
    parameter_names = ['beta','gamma']

    @staticmethod
    def integrate(t, S, I, R, beta, gamma):
        
        # Calculate total population
        N = S+I+R
        # Calculate differentials
        dS = -beta*S*I/N
        dI = beta*S*I/N - 1/gamma*I
        dR = 1/gamma*I

        return dS, dI, dR

# Define parameters and initial condition
params={'beta':0.35, 'gamma':5}
init_states = {'S': 1000, 'I': 1}

# Define the parameters of the observation process
observation_process_params = {'I': 0.04} # overdispersion parameter `alpha` 

# Initialize model
model = ODE_SIR(states=init_states, parameters=params, observation_parameters=observation_process_params)

# Define dataset
data=[d, ]
states=[I,]

# Calibated parameters and bounds
pars = ['beta',]
labels = ['$\\beta$',]
bounds = [(1e-6,1),]

# Setup objective function (no priors --> uniform priors based on bounds)
objective_function = log_posterior_probability(model, pars, bounds, data, states, labels=labels)