1231-Reed-Frost epidemic model.sagews

Reed-Frost epidemic model

In this stochastic model, there are three classes of individuals: Susceptible, Infected, and Recovered (or Removed). At each time step, susceptible individuals are infected by one of the infected individuals with probability p. Infected individuals recover in one time step.

This demo plots the number of infected (red) and susceptible (blue) individuals over time. A histogram of final epidemic sizes (total number infected) is also shown. By default, the code plots the results of 10 stochastic runs. You can set the number of runs, p, and the initial number of susceptible and infected individuals.

#Reed-Frost epidemic model
from scipy import stats

@interact
def reed_frost(S0 = input_box(100, label = 'S[0]: '), I0 = input_box(1, label = 'I[0]: '), p = input_box(0.05, label = 'p: '), n = input_box(10, label = 'trials: ')):
    q = 1-p
    iplot = line([0,0])
    splot = line([0,0])
    finalsizebinsize = 1
    while ((I0+S0)/finalsizebinsize>100):
        finalsizebinsize = finalsizebinsize*10
    while ((I0+S0)/finalsizebinsize>30):
        finalsizebinsize = finalsizebinsize*2
    finalsize = [0 for i in xrange(int((I0+S0)/finalsizebinsize)+1)]
    for i in xrange(n):
        I = I0
        S = S0
        Ilist = [I]
        Slist = [S]
        t = 1
        while (I>0):
            if (S>0):
                I = stats.binom.rvs(S,1-(q^I))
            else:
                I = 0
            S = S-I
            Ilist.append(I)
            Slist.append(S)
        bin = int(sum(Ilist)/finalsizebinsize)
        finalsize[bin] = finalsize[bin]+1
        iplot = iplot + list_plot(Ilist,plotjoined=True,rgbcolor=(1,0,0))
        splot = splot + list_plot(Slist,plotjoined=True,rgbcolor=(0,0,1))
    show(iplot+splot, axes_labels=['time', 'individuals'])
    fseq = IndexedSequence(finalsize,range(int((I0+S0)/finalsizebinsize)+1))
    show(fseq.plot_histogram(), axes_labels=['final size x %d' % finalsizebinsize, 'frequency'])

Interact: please open in CoCalc