<jupyter_start><jupyter_text>Introduction: Forecast for SEAIRQ model with stochastic parameters In this notebook, we consider the SEAIRQ model. We assume that the parameters * $\beta$ (probability of infection on contact),* $\gamma_{E}$ (rate of progression for exposed individual to class A), * $\gamma_{AA}$ (rate of progression from class A to asymptomatic infective class), * $\gamma_{AS}$ (rate of progression from class A to symptomatic infective class), * $\gamma_{I_a}$ (rate of removal for asymptomatic infected individuals), and* $\gamma_{I_s}$ (rate of removal for symptomatic infected individuals) * $ \tau_S$ (quarantining rate for susceptibles)* $ \tau_E$ (quarantining rate for exposed)* $ \tau_A$ (quarantining rate for A)* $ \tau_{I_a}$ (quarantining rate for asymptomatic infectives)* $ \tau_{I_s}$ (quarantining rate for symptomatic infectives)are not known exactly, but rather are characterized by a 11D Gaussian distribution with known mean and covariance matrix. The Gaussian distribution function is trunacted, i.e. set to zero if any parameter is $< 0$.**We now illustrate how uncertainties in the parameters affect the predictions of the SEAIRQ model.**For this we simulate the SEIR model $N_s = 500$ times; for each simulation the above parameters are sampled from a given 11D Gaussian distribution. The resulting 500 trajectories are shown together with their mean, standard deviation, median, and 5 as well as 95 percentiles.We perform this analysis for the deterministic SEAIRQ model.  Define model parameters and initialise pyross.forecast.SEAIRQ<jupyter_code>M = 1                  # the SEAIRQ model we consider has no age structure
Ni = 50000*np.ones(M)   # so there is only one age group 
N = np.sum(Ni)         # and the total population is the size of this age group


E0  = np.array([0])
A0  = np.array([1])
Ia0 = np.array([0])     # the SEAIRQ model we consider has only one kind of infective 
Is0 = np.array([20])     # we take these to be symptomatic 
Q0  = np.array([0])
R0  = np.array([0])     # and assume there are no recovered individuals initially 
S0  = N-(Ia0+Is0+R0+E0)    # The initial susceptibles are obtained from S + E + A + Ia + Is + R = N


# there is no contact structure
def contactMatrix(t):   
    return np.identity(M) 

# duration of simulation and output datapoints
Tf = 500;  Nt=Tf+1


# These parameters we consider exact
fsa   = 1              # the self-isolation parameter 
tE    = 0.00           # rate E  -> Q
tA    = 0.00           # rate A  -> Q
tIa   = 0.00           # rate Ia -> Q
tIs   = 0.05           # rate Is -> Q


# These are the parameters that we sample stochastically
# means
alpha = 0.0            # fraction of asymptomatic infectives 
beta  = 0.2            # infection rate
gIa   = 0.1            # removal rate of asymptomatic infectives 
gIs   = 0.1            # removal rate of symptomatic infectives 
gE    = 0.04           # removal rate of E
gA    = 0.2            # rate to go from A to Ia



# order in covariance matrix:
# beta, gE, gAA, gAS, gIa, gIs, tS, tE, tA, tIa, tIs
#
cov = np.zeros([6,6],dtype=float)
cov[0,0] = 0*alpha**2 # cov(alpha, alpha) = Var(alpha)
cov[1,1] = 0.1*beta**2 # cov(beta, beta) = Var(beta)
cov[2,2] = 0.01*gIa**2 # cov(gIa,gIa) = Var(gIa)
cov[3,3] = 0.01*gIs**2 # cov(gIs,gIs) = Var(gIs)
cov[4,4] = 0.01*gA**2 # cov(gA, gA) = Var(gA)
cov[5,5] = 0.01*gE**2 # cov(gE, gE) = Var(gE)

# 
cov[1,5] = 0.01*beta*gE # cov(beta, gE)
cov[5,1] = cov[1,5] # covariance matrix is symmetric
#
cov[2,3] = cov[2,2] # cov(gIa, gIs)
cov[3,2] = cov[2,3] 


# Define parameters for simulations
parameters = {'alpha':alpha, 'beta':beta,
              'gE':gE,'gA':gA,
              'gIa':gIa, 'gIs':gIs, 'gE':gE, 'fsa':fsa,
             'tE':tE,'tA':tA,'tIa':tIa,'tIs':tIs,
             'cov':cov
                }

# Initialise pyross forecast module
model_forecast = pyross.forecast.SEAIRQ(parameters, M, Ni)


# Number of simulations over which we average, use 500
Ns = 10 
# Define a function which we use below to plot simulation results
def plot_trajectories(result,
                percentile=-1,
                plot_index = 4, # which time series should be plotted? 
               filename='None'): # set filename for saving figures
    # plot_index   class
    #   0          susceptibles
    #   1          exposed
    #   2          asymptomatic and infectious
    #   3          asymptomatic infectives
    #   4          symptomatic infectives
    #   5          quarantined
    if plot_index == 0:
        title='Susceptibles'
        ylabel = r'$N_S$'
    elif plot_index == 1:
        title='Exposed'
        ylabel = r'$N_{E}$'
    elif plot_index == 2:
        title=r'Asymptomatic, infectious (A)'
        ylabel = r'$N_{A}$'
    elif plot_index == 3:
        title='Asymptomatic infectives'
        ylabel = r'$N_{I,a}$'
    elif plot_index == 4:
        title='Symptomatic infectives'
        ylabel = r'$N_{I,s}$'
    elif plot_index == 5:
        title='Quarantined'
        ylabel = r'$N_{Q}$'
    else:
        raise RuntimeError("plot_index should be 0, 1, 2, or 3.")
    #
    fontsize=25
    #
    #
    trajectories = result['X']
    t_arr = result['t']
    traj_mean = result['X_mean']
    traj_std = result['X_std']
    #
    #
    # Plot trajectories
    #
    fig, ax = plt.subplots(1,1,figsize=(7,5))
    ax.set_title(title,
                 y=1.05,
                fontsize=fontsize)
    for i,e in enumerate(trajectories):
        ax.plot(t_arr,e[plot_index],
                alpha=0.15,
               )
    ax.fill_between(t_arr,traj_mean[plot_index] - traj_std[plot_index],
                    traj_mean[plot_index] + traj_std[plot_index],
                    alpha=0.7,
                    color='limegreen',
                           label='Std deviation')
    ax.plot(t_arr,traj_mean[plot_index] - traj_std[plot_index],
                    alpha=1,
            label='Std deviation',
            lw=1.5,
            ls='--',
                    color='black')
    ax.plot(t_arr,traj_mean[plot_index] + traj_std[plot_index],
                    alpha=1,
            #label='Std deviation',
            lw=1.5,
            ls='--',
                    color='black')
    ax.plot(t_arr,traj_mean[plot_index],
                    alpha=1,
            lw=2,
            color='black',
           label='Mean')
    ax.set_xlim(np.min(t_arr),np.max(t_arr))
    ax.set_ylabel(ylabel,fontsize=fontsize)
    ax.set_xlabel(r'$t$ [days]',fontsize=fontsize)
    ax.legend(loc='upper right',fontsize=18)
    plt.show()
    if filename != 'None':
        fig.savefig(filename + '_trajs.png', bbox_inches='tight',dpi=100)
    plt.close()
    #
    #
    #
    # Plot percentiles
    #
    if percentile > 0:
        percentiles_lower = np.percentile(trajectories[:,plot_index],percentile,axis=0)
        percentiles_upper = np.percentile(trajectories[:,plot_index],100-percentile,axis=0)
        percentiles_median = np.percentile(trajectories[:,plot_index],50,axis=0)
        print("In the following plot, red dashed lines denote {0} and {1} percentiles of the numerical data:".format(percentile,
                                                                                             100-percentile))
        fig, ax = plt.subplots(1,1,figsize=(7,5))
        ax.set_title(title,
                     y=1.05,
                    fontsize=fontsize)
        for i,e in enumerate(trajectories):
            ax.plot(t_arr,e[plot_index],
                    alpha=0.15,
                   )
        ax.fill_between(t_arr,percentiles_lower,
                        percentiles_upper,
                        alpha=0.1,
                        color='red',
                               label='Percentiles')
        ax.plot(t_arr,percentiles_lower,
                        alpha=1,
                lw=2,
                label='Percentiles',
                ls='--',
                color='red',
               )
        ax.plot(t_arr,percentiles_upper,
                        alpha=1,
                lw=2,
                color='red',
                ls='--',
               )
        ax.plot(t_arr,percentiles_median,
                        alpha=1,
                lw=2,
                color='red',
               label='Median')
        ax.plot(t_arr,traj_mean[plot_index],
                        alpha=1,
                lw=2,
                color='black',
               label='Mean')
        ax.set_xlim(np.min(t_arr),np.max(t_arr))
        ax.set_ylabel(ylabel,fontsize=fontsize)
        ax.set_xlabel(r'$t$ [days]',fontsize=fontsize)
        ax.legend(loc='upper right',fontsize=18)
        plt.show()
        if filename != 'None':
            fig.savefig(filename + '_trajs2.png', bbox_inches='tight',dpi=100)
        plt.close()
        
        
# Define a function which we use below to plot parameters used for simulations
def plot_sample_parameters(result,
               filename='None'): # set filename for saving figures
    #
    fontsize=25
    #
    # Scatterplot of used parameters
    #
    sample_parameters = result['sample_parameters'].T
    beta = result['beta']
    gE = result['gE']
    gIa = result['gIa']
    gIs = result['gIs']
    #
    title = r'Samples for stochastic $\beta$, $\gamma_{E}$'
    labelx = r'$\beta $'
    labely = r'$\gamma_{E}$'
    x_mean = beta
    y_mean = gE
    labelx_mean = r'$\langle \beta \rangle$'
    labely_mean = r'$\langle \gamma_{E} \rangle$'
    data_index_x = 1
    data_index_y = 4
    fig, ax = plt.subplots(1,1,figsize=(7,5))
    ax.set_title(title,y=1.05,fontsize=fontsize)
    ax.axvline(x_mean,color='limegreen',ls='--',lw=2,label=labelx_mean)
    ax.axhline(y_mean,color='dodgerblue',ls='--',lw=2,label=labely_mean)
    ax.scatter(sample_parameters[data_index_x], sample_parameters[data_index_y] ,
               label='sampled data',
               color='black',s=10) #, c = truth)
    ax.set_xlabel(labelx,fontsize=fontsize)
    ax.set_ylabel(labely,fontsize=fontsize)
    ax.set_xlim(0,1.05*np.max(sample_parameters[data_index_x]))
    ax.set_ylim(0,1.05*np.max(sample_parameters[data_index_y]))
    ax.legend(loc='best',fontsize=15)
    plt.show()
    if filename != 'None':
        fig.savefig(filename + '_samples1.png', bbox_inches='tight',dpi=100)
    plt.close()
    #<jupyter_output><empty_output><jupyter_text>Forecast based on deterministic model<jupyter_code>result = model_forecast.simulate(S0, E0, A0, Ia0, Is0, Q0,
                                 contactMatrix, Tf, Nt,
                                      verbose=True,
                                     Ns=Ns) 

plot_trajectories(result,
            plot_index = 2,
            percentile=5,
           )

plot_trajectories(result,
#            filename='forecast_deterministic',
            percentile=5,
           )

plot_trajectories(result,
            plot_index = 5,
            percentile=5,
           )

plot_sample_parameters(result)<jupyter_output>Finished. Time needed for evaluation: 00:00:05
StarcoderdataJupyterStructuredCleanDedup