Aula 12 - Modelagem Estocástica.ipynb

Kernel: SageMath 9.7

Conceitos básicos de processos estocásticos em populações

Copyright (2013-15) Flávio Codeço Coelho, EMAp-FGV

Neste módulo revisitamos os elementos básicos da dinâmica de populações, como nascimento, morte, imigração e emigração. Mas agora vamos representá-los como realizações de eventos probabilísticos.

Processo Poisson simples

Neste processo consideramos apenas o nascimento, ou mais genericamente a entrada de indivíduos em uma população de forma aleatória, mas a uma taxa constante $\alpha$.

Seja $N(t, t+dt)$ o número de novos indivíduos que entraram na população no intervalo $(t, t+dt)$. Então podemos definir as probabilidades de transição da população como:

$(i)\, Pr\{N(t, t+dt)=0\}=1-\alpha dt + o(dt)$,

$(ii)\, Pr\{N(t, t+dt)=1\}= \alpha dt +o(dt)$

$(iii)\, Pr\{N(t, t+dt)>1\}= o(dt)$

$(iv)\, N(t, t+dt)$ é independente de todas as entradas no intervalo $(0, t)$

Aqui, o(dt) denota todos os termos  de order inferior a dt.

Tempo até o próximo evento

Suponha que $t_0 + Z$ seja o tempo do primeiro evento após $t_0$. Para determinar a distribuição de probabilidade de $Z$, podeŕiamos considerar a distribuição $F(t)=Pr(Z\leq t) (t \geq 0)$. Mas é muito mais fácil pensarmos sobre o evento complementar $G(t)=1-F(t)= Pr(Z>t)$, ou seja nenhum evento ocorreu antes de $t$.

Considere um intervalo pequeno $dt >0$, agora vamos dividir os eventos que ocorrem no intervalo $(t_0, t_0+t+dt)$ nos que ocorrem no primeiro intervalo $(t_0, t_{0+t})$ e nos que ocorrem no intervalo subsequente $(t_{0+t}, t_{0+t+dt})$. Podemos ver que

$G(t+dt)=Pr\{Z>t+dt\}$

$=Pr\{Z>t$ e nenhum evento ocorre no intervalo $(t_{0+t}, t_{0+t+dt})\}$ 

$=Pr\{Z>t\} Pr\{$nenhum evento ocorre no intervalo $(t_{0+t}, t_{0+t+dt})\mid Z>t\}$ 

aplicando a condição $(iv)$ acima, 

$G(t+dt) =Pr\{$nenhum evento ocorre no intervalo $(t_{0+t}, t_{0+t+dt})\}$

Aplicando a condição $(i)$ obtemos

$G(t+dt)=G(t)(1-\alpha dt + o(dt)$

Quando $dt$ tende a $0$, temos

$\frac{dG(t)}{dt}= -\alpha G(t)$

Cuja solução é

$G(t)=G(0) e^{-\alpha t}$

Como $G(0)=Pr(Z>0)=1$ concluímos que a função de distribuição acumulada de $Z$ é

$F(t) = Pr(Z \leq t)=1-G(t)= 1- e^{-\alpha t}$

Logo a pdf do tempo até o próximo evento é:  

$f(t)=F'(t)=\alpha e^{-\alpha t} \,\, (\alpha > 0)$

que corresponde à distribuição exponencial com parâmetror $\alpha$.

Exercício: Simulando o processo de Poisson

A partir do exposto acima, implemente uma simulação de 100 eventos  de um processo de Poisson com taxa $\alpha$. Mostre graficamente que a distribuição do tempo até o centésimo evento segue uma distribuição Gamma(100,$\alpha$). 

import numpy as np
from numpy.random import exponential as expv
import pylab as P
from numpy import cumsum
import warnings
warnings.filterwarnings('always')
%display typeset
%matplotlib notebook

alpha = 1.5
times = cumsum(expv(alpha, size=100))
list_plot(list(zip(times, np.ones(100))))

from numpy.random import gamma as gammav
t100 = [cumsum(expv(alpha, size=100))[-1] for i in range(1000)]
gs = gammav(100, alpha, size=1000)
P.figure(figsize=(15,8))
P.hist(gs, bins=50, color='r', alpha=0.3, label='Gamma(100, %1.1f)'%alpha);
P.hist(t100,bins=50, alpha=.3, label='$T_{100}$');
P.title('1000 amostras')
P.legend();

<IPython.core.display.Javascript object>

Equações Diferenciais Estocásticas (EDEs)

 Uma equação diferencial estocástica é uma equação diferencial na qual um ou mais de seus termos é um processo estocástico. Por conseguinte, a sua solução também é um processo estocástico. Seguindo a notação do cálculo de Ito, um exemplo simples seria:

$$dX_t = a(X_t) \, d t + b(X_t) \, d W_t,$$

 

onde $W_t$ é um processo de Wiener, ou seja, um processo estocástico em tempo contínuo, também conhecido como movimento browniano.

A solução analítica destas equações foge ao escopo deste curso. Mas podemos explorar soluções numéricas.

Modelando com EDEs

Vamos agora desenvolver alguns modelos simples usando EDEs.

Dinâmica de uma população (Processo de nascimentos e mortes)

Seja  o seguinte modelo da dinâmica demográfica de uma população:

$$\frac{\mathrm{d}x}{\mathrm{d}t}=b x - d x; \,\,\,\,\,\, x(0)=x_0$$

onde $x(t)$ é a população no tempo $t$ e $b$ e $d$ são as taxas de nascimento e morte, respectivamente. A solução deste sistema é:

$$x(t)= x_0e^{(b-d)t}$$

 

var('t')
b = 0.3
d = 0.27
x0 = 25
x(t) = x0*exp((b-d)*t)
plot(x(t),(t,0,50))

Naturalmente, populações naturais não variam de forma tão bem comportada. Como existem heterogeneidades em uma população real, $x(t)$ como derivado acima corresponde ao valor esperado da população, e não ao valor exacto de um processo de nascimentos e mortes no tempo $t$. Além do mais $x(t)$ varia continuamente enquanto que populações variam de forma discreta.

Um modelo estocástico da demografia desta população, preservando as mesmas taxas de nascimento e morte pode ser construído, nos oferecendo uma série temporal mais realística.

 

Para nos ajudar a chegar na EDE correspondente, vamos construir um modelo estocástico direto primeiro. Vamos considerar os nascimentos e mortes aleatórios, em um curto intervalo de tempo $\Delta t$.

Neste intervalo, são possíveis três tipos de alterações na população:

1. Um nascimento, $\beta_1=1$,
2. uma morte, $\beta_2=-1$,
3. ou nada, $\beta_3=0$.

Estas alterações ocorrerão com  as  respectivas probabilidades:

1. $P_1 = bx\Delta t$,
2. $P_2 = dx\Delta t$ e
3. $P_3 = 1-(b+d)x\Delta t$.

Seja $p_k(t) = P(X(t) = x_k)$ a probabilidade de ser ter $x_k$ indivíduos na população no instante de tempo $t$.

Isto nos leva a:

$$p_k(t+\Delta t) = p_k(t)\left[1 -bx_k\Delta t -dx_k\Delta t\right]+p_{k-1}(t)\left[bx_{k-1}\Delta t\right]+p_{k+1}(t)\left[dx_{k+1}\Delta t\right]$$

Pode-se demonstrar que para $\Delta t$ muito pequeno, que a EDE abaixo satisfaz aproximadamente a mesma distribuição de probabilidade que o processo estocástico discreto acima.

$$\mathrm{d}x(t) = (b - d)x \mathrm{d}t + \sqrt{((b+d)x)}\mathrm{d}W(t)$$

onde $W(t)$ é um processo de Wiener.

Para resolver numericamente esta EDE precisamos entender o método de Euler-Maruyama.

Método de Euler-Maruyama

Considere a seguinte equação diferencial estocástica de  Itō :

$$\mathrm{d}X_t = a(X_t) \, \mathrm{d} t + b(X_t) \, \mathrm{d} W_t,$$

Com condições iniciais $X_0=x_0$, onde where $W_t$ é o processo de Wiener, e suponha que queremos resolver esta EDE em um dado intervalo $[0, T]$. Então a aproximação de Euler–Maruyama approximation à solução verdadeira de $X$ é a cadeia de Markov $Y$ definida abaixo:

• Divida o interval $[0,T]$ em $N$ sub-intervalos iguais de comprimento $\Delta T \gt 0$; $$0 = \tau_{0} < \tau_{1} < \cdots < \tau_{N} = T \mbox{ e } \Delta t = T/N;$$
• faça $Y_0=x_0;$
• Defina $Y_n$ recursivamente para $1 \leq n \leq N$: $$\, Y_{n + 1} = Y_n + a(Y_n) \Delta t + b(Y_n) \Delta W_n,$$ onde $$\Delta W_n=W_{\tau_{n+1}}-W_{\tau_n}.$$ $\Delta W \sim N(0,\Delta t)$, i.e. i.i.d. com média $0$ e variância $\Delta t$

Solução numérica do processo de nascimento e morte

Vamos agora definir as funções que correspondem às partes determinísticas e estocásticas da nossa EDE.

determ = lambda x,b,d:(b-d)*x
stoc = lambda x,b,d: sqrt((b+d)*x)

Agora definimos uma função que implementa o método de Euler-Maruyama, Junto com uma figura interativa aplicando-o ao nosso modelo.

def EulerMaruyama(tstart, ystart, tfinish, nsteps, f1, f2, params):
    sol = np.zeros(nsteps)
    tvals = np.empty(nsteps)
    sol[0] = ystart 
    tvals[0] = tstart 
    h = (tfinish-tstart)/nsteps
    for step in range(nsteps-1): 
        sol[step+1]=sol[step] + h*f1(sol[step],*params) +  f2(sol[step],*params)*normalvariate(0,h)
        tvals[step+1]=(tvals[step] + h) 
    return list(zip(tvals,sol))
    
out = Graphics()
save(out,'temp')

@interact
def EulerMaruyamaExample(b = slider(srange(0,1,.01),default=0.3),
                        d = slider(srange(0,1,.01),default=0.27),
                        plots_at_a_time = slider(range(1,100),default=10), 
                        number_of_steps = slider(range(1,1000),default=500), 
                        clear_plot = checkbox(True), 
                        auto_update=False):
    pretty_print(html(r'<center>Solução pelo método de Euler-Maruyama <br>da EDE:</center>'))
    pretty_print(html(r'<center>$\mathrm{d}x(t) = (b - d)x \mathrm{d}t + \sqrt{((b+d)x)}\mathrm{d}W(t)$</center>'))
    sol = EulerMaruyama(0,25,50,number_of_steps,determ,stoc,(b,d))
    emplot = list_plot(sol,plotjoined=True, color='cyan')
    exact_plot = plot(x(t),(t,0,50), color='red', legend_label=u'Solução determinística')
    for i in range(1,plots_at_a_time):
        sol2 = EulerMaruyama(0,25,50,number_of_steps,determ,stoc,(b,d))
        sol = [(v[0],sol[n][1]+v[1]) for n, v in enumerate(sol2)]
        emplot = emplot + list_plot(sol2,plotjoined=True, alpha=0.3)
    
    if clear_plot:
        avgplot = list_plot([(i,v/plots_at_a_time) for i,v in sol],plotjoined=True, color='green', legend_label=u"média")
        out = emplot + exact_plot + avgplot
        save(out,'temp')
    else:
        out = load('temp')
        out = out + emplot + exact_plot
        save(out,'temp')
    show(out, figsize = [8,5])

O Método de Milstein

O método de Milstein parte da mesma ideia de do método de EM, ou seja discretizar a dinâmica em n intervalos iguais. A forma como serão calculados os valores da função em função de seus valores anteriores é que muda:

$$Y_{{n+1}}=Y_{n}+a(Y_{n})\Delta t+b(Y_{n})\Delta W_{n}+{\frac {1}{2}}b(Y_{n})b'(Y_{n})\left((\Delta W_{n})^{2}-\Delta t\right)$$

onde $b'$ é a derivada de $b(x)$ com respeito a $x$, e $\Delta W_{n}=W_{{\tau _{{n+1}}}}-W_{{\tau _{n}}}$,  são variáveis aleatórias i.i.d. de distribuição Normal com média 0 e variância $\Delta t$.

Exercício:

Modifique o solver acima implementando o método de Milstein. 

def milstein(tstart, ystart, tfinish, h, f1, f2, params):
    nsteps = int((tfinish-tstart)/h)
    sol = np.zeros(nsteps)
    tvals = np.zeros(nsteps)
    sol[0] = ystart 
    tvals[0] = tstart 
    def deriv(f,x, params,h=0.0001):
        if x-h <0:
            h=x
        return (f(x+h,*params)-f(x-h, *params))/(2*h)
    
    for step in range(nsteps-1):
        R = normalvariate(0,h)
        sol[step+1]=sol[step] + h*f1(sol[step],*params) +  f2(sol[step],*params)*R + 0.5*f2(sol[step],*params)*deriv(f2,sol[step],params) * (R**2-h)
        tvals[step+1]=(tvals[step] + h)
    return list(zip(tvals,sol))

out = Graphics()
save(out,'temp')

@interact
def MilsteinExample(b = slider(srange(0,1,.01),default=0.3),
                    d = slider(srange(0,1,.01),default=0.27),
                    tfinish = slider(range(10,1000), default=50),
                    plots_at_a_time = slider(range(1,100),default=10), 
                    passo = slider(srange(0,.5,.001),default=0.25), 
                    clear_plot = checkbox(True), 
                    auto_update=True):
    pretty_print(html(r'<center>Solução pelo método de Milstein <br>da EDE:</center>'))
    pretty_print(html(r'<center>$\mathrm{d}x(t) = (b - d)x \mathrm{d}t + \sqrt{((b+d)x)}\mathrm{d}W(t)$</center>'))
    sol = milstein(0,25,tfinish,passo,determ,stoc,(b,d))
    mplot = list_plot(sol,plotjoined=False, color='cyan')
#     show(mplot)
    exact_plot = plot(25*exp((b-d)*t),(t,0,tfinish), color='red', legend_label=u'Solução determinística')
    for i in range(1,plots_at_a_time):
        sol2 = milstein(0,25,tfinish,passo,determ,stoc,(b,d))
        sol = [(v[0],sol[n][1]+v[1]) for n, v in enumerate(sol2)]
        mplot = mplot + list_plot(sol2,plotjoined=False, alpha=0.3)
    
    if clear_plot:
        avgplot = list_plot([(i,v/plots_at_a_time) for i,v in sol],plotjoined=False, color='green', legend_label=u"média")
        out = mplot + exact_plot + avgplot
        save(out,'temp')
    else:
        out = load('temp')
        out = out + mplot + exact_plot
        save(out,'temp')
    show(out, figsize = [8,5])

Derivando uma EDE para um modelo SIR

seja $\Delta X(t)=(\Delta S, \Delta I)^T$, onde $\Delta X_i = X_i(t+\Delta t) - X_i(t)$, $i=1,2$.

Então o valor esperado de $\Delta X(t)$ na escala detempo de $\Delta t$ é: $$E(\Delta X(t)) = \begin{pmatrix}-\frac{\beta}{N} S I\\ \frac{\beta}{N} S I -\gamma I\end{pmatrix} \Delta t.$$ A matriz de covariância de $\Delta X(t)$ é $V(\Delta X(t))\approx E(\Delta X(t)[\Delta X(t)]^T$, ou seja, $$V(\Delta x(t))=\begin{pmatrix}\frac{\beta}{N} S I & -\frac{\beta}{N} S I\\ -\frac{\beta}{N} S I & \frac{\beta}{N} S I + \gamma I \end{pmatrix}\Delta t.$$

A partir dos resultados acima podemos chegar ao seguinte sistema de equações diferenciais estocásticas para o modelo SIR:

$$\begin{align} \begin{matrix}\frac{dS}{dt}&=-\frac{\beta}{N} S I+B_{11}\frac{dW_1}{dt}+B_{12}\frac{dW_2}{dt}\\ \frac{dI}{dt} &= \frac{\beta}{N} S I -\gamma I+B_{21}\frac{dW_1}{dt}+B_{22}\frac{dW_2}{dt} \end{matrix} \end{align}$$

onde $B=\sqrt{V(\Delta X(t))}$.

Equações de Kolmogorov

Também conhecidas como "Markov Jump processes", estas equações descrevem um processo markoviano contínuo e diferenciável em um espaço de estados discreto e finito. A equações de Kolmogorov (anterógrada e retrógrada), modelam o problema de uma forma ligeiramente diferente. Nas EDE modelamos um sistema dinâmico estocástico como a soma de um processo médio, $a(X_t)dt$, e sua covariância, $b(X_t)dW$. Nas equações de Kolmogorov, modelamos a taxa de variação das probabilidaes de transição. Logo sua versão anterógrada é:

$$\frac{dP(t)}{dt} =QP(t),$$

onde $P(t)=(p_{ji}(t)$ é a matriz de probabilidades de transição e $Q=(q_{ji})$, é a matriz geradora. De forma similar, a versão retrógrada é dada por:

$$\frac{dP(t)}{dt} =P(t)Q$$

 

Os métodos para solução numérica deste tipo de processo estocástico foram originalmente desenvolvidos para modelar reações químicas, mas podem ser aplicados a quaisquer modelos compartimentais cujo espaço de estados seja discreto. Estes métodos se dividem em exatos e aproximados. Dentre os métodos exatos o método "direto" proposto por GIllespie em 1977 é o mais usado. Resumidamente, o método de gillespie é um método iterativo que faz uso de duas variáveis aleatórias: uma representando o tempo até o próximo evento, e a segunda. multinomial, para escolher qual dos evento ocorrerá.

Para descrever esta classe de modelos, é interessante pensar no modelo como um vetor de estados

$$X(t)=[X_1(t), X_2(t), \ldots, X_N(t)]$$

Onde $X_n(t)$ é a quantidade do $n$-ésimo tipo no tempo $t$.

Existem $M$ reacões, envolvendo os tipos existentes, que podem ocorrer alterando o estado do sistema. Cada reação é caracterizada por uma função de propensão $a_m$ e um vetor de alteração de estados $V_m=[V_{m1}, \ldots , V_{mN} ]$.

 

Implementando processos markovianos em tempo contínuo

Abaixo vamos aprender a simular um modelo SIR estocástico simples usando apenas as bibliotecas numpy e scipy. Vamos simular a versão determinística e também algumas realizações da versão estocástica.

import numpy as np
import pandas as pd
from scipy.integrate import odeint
from numpy.random import rand, gamma, exponential, poisson, multinomial
import pylab as P


def ode(y,t, *parms):
    S,I = y
    beta, gam = parms

    return[-beta/N*S*I,
            beta/N*S*I - gam*I
    ]




  
def run_sir(N, tf, nsims, *pars):
    """
    Runs simulation.

    :parameters:
    :param N: Population size
    :param tf: Stop time
    :param nsims:  Number of simulations
    :param pars: epidemic parameters 
    """
    beta, gam, I0 = pars

    sims = {}
    for k in range(nsims):
        t = [0]
        S = [N-I0]
        I = [I0]
        gts =[]
        while I[-1] > 0 and t[-1] < tf:
            U = rand()

            R = beta*S[-1]*I[-1]/N + gam*I[-1]
            pinf = ((beta/N)*S[-1]*I[-1])/R # Probability of next event being an infection
            gt = exponential(1/R)

            if U <= pinf: # Next event is an infection
                
                S.append(S[-1]-1)
                I.append(I[-1]+1)
                t.append(t[-1] + gt)
                gts.append(gt)
            else: # next event is a recovery
                S.append(S[-1])
                I.append(I[-1]-1)
                #print('removal')
                t.append(t[-1] + gt)  # -np.log(rand())/R)
        sims[k] = (np.array([t,S,I]).T, np.array(t), gts)
        P.plot(t,I,label='$O_t^{}$'.format(k+1), drawstyle='steps-post')
    return sims

beta=0.1
gam = 1/21
N = 500
constant = False
I0 = 1
tf = 365
ts = np.arange(tf)
nsims = 10
P.figure(figsize=(15,10))
sims = run_sir(N,tf,nsims, beta, gam, I0)
Is = [i[0][:,2] for i in sims.values()]
EI = odeint(ode, [N-2,2], ts, args=(beta,gam))
P.plot(ts, EI[:,1], 'k-',lw=2, label='$Y_t$')
P.legend(loc=0)
P.xlabel("t (days)")
P.ylabel("cases")
P.grid()
P.savefig('simulation.png')

<IPython.core.display.Javascript object>

WARNING: Some output was deleted.

Agora vamos simular uma versão extendida do modelo SIR com demografia, ou seja, com nascimentos e mortes.

def SIRd(y,t, *parms):
    S,I = y
    alpha, beta, gam, N = parms
    betat = lambda t: beta+(0.0*beta)*np.cos((2*np.pi*t)/365.)
    beta = betat(t)
    return[alpha*N-beta*S*I/N -alpha*S,
            beta*S*I/N - gam*I -alpha*I
    ]


def sir_dem(N, tf, nsims, *pars):
    """
    Runs simulation.

    :parameters:
    :param N: Population size
    :param tf: Stop time
    :param nsims:  Number of simulations
    :param gtd: Generation time distribution: 'e' for exponential anything else for gamma
    :param pars: epidemic parameters 
    """
    alpha, beta, gam, I0, Tmed, constant = pars
    R0=beta/(gam+alpha)
    sims = {}
    for k in range(nsims):
        t = [0]
        S = [N-I0]
        I = [I0]
        gts =[]
        while I[-1] > 0 and t[-1] < tf:
            T =t[-1]
            R = alpha*N+beta*S[-1]*I[-1]/N + gam*I[-1] + alpha*S[-1] + alpha*I[-1]
            pbirth = alpha*N/R  # Probability of the next event being a birth (S -> S+1)
            pinf = ((beta/N)*S[-1]*I[-1])/R # Probability of next event being an infection
            prec = gam*I[-1]/R  # Probability of the next event being a recovery (I -> I-1)
            pds = alpha*S[-1]/R
            pdi = alpha*I[-1]/R
            ev = multinomial(1, [pbirth, pinf, prec, pds, pdi]).nonzero()
            gt = exponential(1/R)
            if ev[0][0] == 0: #  event is a birth
                S.append(S[-1]+1)
                I.append(I[-1])
            elif ev[0][0] == 1:  #  event is an infection
                S.append(S[-1]-1)
                I.append(I[-1]+1)
            elif ev[0][0] == 2:  # event is a recovery
                S.append(S[-1])
                I.append(I[-1]-1)
            elif ev[0][0] == 3:  # event is a susceptible death
                S.append(S[-1]-1)
                I.append(I[-1])
            elif ev[0][0] == 4: # next event is a infectious death
                S.append(S[-1])
                I.append(I[-1]-1)
            t.append(t[-1] + gt)
            gts.append(gt)
        sims[k] = (np.array([t,S,I]).T, np.array(t), gts)
        P.plot(t,I,label='$O_t^{}$'.format(k+1), alpha=0.3, drawstyle='steps-post')
        #P.plot(t,R0*(np.array(S)/N),label='$R_t^{}$'.format(k+1), alpha=0.3, drawstyle='steps-post')
    return sims

alpha = 0.01
beta=0.1
gam = 1/21
R0 = beta/(alpha+gam)
N = 3000
I0 = 15
S0 = N-I0
print("R0={}".format(R0))
print("Endemic Equilibrium: S: {}: I:{}".format(S0*(alpha+gam)/beta, S0*alpha/beta*(R0-1)))
Tmed = 20
constant = False
P.figure(figsize=(15,10))
tf = 3650//2
ts = np.arange(tf)
nsims = 3
sims = sir_dem(N,tf,nsims, alpha, beta, gam, I0, Tmed, constant)
Is = [i[0][:,2] for i in sims.values()]
EI = odeint(SIRd, [N-I0,I0], ts, args=(alpha,beta,gam, N))

P.plot(ts, EI[:,1], 'k-',lw=2, label='$I(t)$')
P.plot(ts, EI[:,0], 'b-',lw=2, label='$S(t)$')
P.ylabel("cases")
P.legend(loc='center')
ax2 = P.gca().twinx()
ax2.plot(ts,EI[:,0]/N*R0, label='$R_t$')
for i in range(nsims):
    ax2.plot(sims[i][0][:,0], R0*sims[i][0][:,1]/N,label='$R_t^{}$'.format(1+1), alpha=0.3, linestyle=':')
ax2.plot()
ax2.set_ylabel('$R_t$')
P.legend(loc=0)
P.xlabel("t (days)")
P.grid()
P.savefig('simulation_dem.png')

R0=1.73553719008264
Endemic Equilibrium: S: 1719.92857142857: I:219.557851239669

<IPython.core.display.Javascript object>

WARNING: 1 intermediate output message was discarded.

Simulando modelos estocásticos em Python

A abordagem acima demonstra como é simples implementar uma simulação estocástica a partir das equações de Kolmogorov. Contudo podemos simplificar ainda mais a nossa vida, usando uma biblioteca de que implementa o algoritmo "direto" de Gillespie, assim como algumas de suas variações: o Stochpy. Para executar o código abaixo, é necessário instalar o stochpy no Sage, uma vez que este não faz parte da distribuição padrão do Sage.

!pip3 install --user stochpy

Requirement already satisfied: stochpy in /home/fccoelho/.sage/local/lib/python3.9/site-packages (2.4)
WARNING: You are using pip version 21.0.1; however, version 21.3 is available.
You should consider upgrading via the '/home/fccoelho/Downloads/SageMath/local/bin/python3 -m pip install --upgrade pip' command.

import stochpy;
import os
import pylab as P

Para escrever um modelo no Pysces, precisamos utilizar a Pysces Modeling language, vejamos como fica um simples modelo SIR

modelpath = '/home/fccoelho/Stochpy/pscmodels'
sir = """
Modelname: SIR
Description: PySCes Model Description Language Implementation of SIR model

# Set model to run with numbers of individuals
Species_In_Conc: False
Output_In_Conc: False

# Differential Equations as Reactions
Infeccao:
    S > I
    beta*S*I/(S+I+R)

Recuperacao:
    I > R
    gamma*I

# Parameter Values
S = 100
I = 5
R = 0
beta = 0.50
gamma = 0.1
"""
with open(os.path.join(modelpath,'SIR.psc'),'w') as f:
    f.write(sir)

Agora vamos definir o nosso método de solução numérica, e carregar o modelo:

smod = stochpy.SSA()
smod.Model('SIR.psc')

# os.system('rm -f '+DATA+'*.dat')
smod.DoStochSim(method="Tau-leap", end = int(100),mode = 'time',trajectories = int(10))

smod.PlotAverageSpeciesTimeSeries(species2plot=["S","I","R"])
P.grid()

*** WARNING ***: No regular grid is created yet. Use GetRegularGrid(n_samples) if averaged results are unsatisfactory (e.g. more or less 'samples')

smod.PlotSpeciesTimeSeries(species2plot=["I"])
P.grid()

smod.PlotSpeciesDistributions(linestyle= 'solid')

####################
just ran this line
just ran this line
####################
####################
just ran this line
just ran this line
####################
####################
just ran this line
just ran this line
####################
####################
just ran this line
just ran this line
####################
####################
just ran this line
just ran this line
####################
####################
just ran this line
just ran this line
####################
####################
just ran this line
just ran this line
####################
####################
just ran this line
just ran this line
####################
####################
just ran this line
just ran this line
####################
####################
just ran this line
just ran this line
####################
####################
just ran this line
just ran this line
####################
####################
just ran this line
just ran this line
####################
####################
just ran this line
just ran this line
####################
####################
just ran this line
just ran this line
####################
####################
just ran this line
just ran this line
####################
####################
just ran this line
just ran this line
####################
####################
just ran this line
just ran this line
####################
####################
just ran this line
just ran this line
####################
####################
just ran this line
just ran this line
####################
####################
just ran this line
just ran this line
####################
####################
just ran this line
just ran this line
####################
####################
just ran this line
just ran this line
####################
####################
just ran this line
just ran this line
####################
####################
just ran this line
just ran this line
####################
####################
just ran this line
just ran this line
####################
####################
just ran this line
just ran this line
####################
####################
just ran this line
just ran this line
####################
####################
just ran this line
just ran this line
####################
####################
just ran this line
just ran this line
####################
####################
just ran this line
just ran this line
####################

smod.PlotSpeciesTimeSeries(species2plot=["S","I","R"])