SharedMM4 - Aula 11 - Modelagem Estocástica.ipynbOpen in CoCalc
Author: Flávio Coelho
Views : 40
License: Apache License 2.0

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)N(t, t+dt) o número de novos indivíduos que entraram na população no intervalo (t,t+dt)(t, t+dt). Então podemos definir as probabilidades de transição da população como:

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

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

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

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

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

Tempo até o próximo evento

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

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

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

=Pr{Z>t=Pr\{Z>t e nenhum evento ocorre no intervalo (t0+t,t0+t+dt)}(t_0+t, t_0+t+dt)\} 

=Pr{Z>t}Pr{=Pr\{Z>t\} Pr\{nenhum evento ocorre no intervalo (t0+t,t0+t+dt)Z>t}(t_0+t, t_0+t+dt)\mid Z>t\} 

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

G(t+dt)=Pr{G(t+dt) =Pr\{nenhum evento ocorre no intervalo

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

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

Quando dtdt tende a 00, temos

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

Cuja solução é

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

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

F(t)=Pr(Zt)=1G(t)=1eαtF(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)=αeαt(α>0)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). 

In [6]:
import numpy as np from numpy.random import exponential as expv import pylab as P from numpy import cumsum %matplotlib inline
In [2]:
alpha = 1.5 times = cumsum(expv(alpha, size=100)) list_plot(zip(times, np.ones(100)))
In [3]:
from numpy.random import gamma as gammav t100 = [cumsum(expv(alpha, size=100))[-1] for i in xrange(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();

%html

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:

dXt= a(Xt)dt+b(Xt)dWt,\mathrm{d}X_t = a(X_t) \, \mathrm{d} t + b(X_t) \, \mathrm{d} W_t,

 

onde WtW_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:

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

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

x(t)=x0e(bd)tx(t)= x_0e^{(b-d)t}

 

In [3]:
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)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 tt. Além do mais x(t)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 Δt\Delta t.

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

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

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

  1. P1=bxΔtP_1 = bx\Delta t,
  2. P2=dxΔtP_2 = dx\Delta t e
  3. P3=1(b+d)xΔtP_3 = 1-(b+d)x\Delta t.

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

Isto nos leva a:

pk(t+Δt)=pk(t)[1 bxkΔt dxkΔt]+pk1(t)[bxk1Δt]+pk+1(t)[dxk+1Δt]p_k(t+\Delta t) = p_k(t)[1 - bx_k\Delta t - dx_k\Delta t]+p_{k-1}(t)[bx_{k-1}\Delta t]+p_{k+1}(t)[dx_{k+1}\Delta t]

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

dx(t)=(bd)xdt+((b+d)x)dW(t)\mathrm{d}x(t) = (b - d)x \mathrm{d}t + \sqrt{((b+d)x)}\mathrm{d}W(t)

onde W(t)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ō :

dXt= a(Xt)dt+b(Xt)dWt,\mathrm{d}X_t = a(X_t) \, \mathrm{d} t + b(X_t) \, \mathrm{d} W_t,

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

  • Divida o interval [0,T][0,T] em NN sub-intervalos iguais de comprimento ΔT>0\Delta T \gt 0;
  • faça Y0=x0;Y_0=x_0;
  • Defina YnY_n recursivamente para 1nN1 \leq n \leq N: Yn+1=Yn+a(Yn)Δt+b(Yn)ΔWn,\, Y_{n + 1} = Y_n + a(Y_n) \Delta t + b(Y_n) \Delta W_n, onde ΔWn=Wτn+1Wτn.\Delta W_n=W_{\tau_{n+1}}-W_{\tau_n}. ΔWN(0,Δt)\Delta W \sim N(0,\Delta t), i.e. i.i.d. com média 00 e variância Δt\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.

In [4]:
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.

In [7]:
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 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('<center>Solução pelo método de Euler-Maruyama <br>da EDE:</center>')) pretty_print(html('<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,100,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:

Yn+1=Yn+a(Yn)Δt+b(Yn)ΔWn+12b(Yn)b(Yn)((ΔWn)2Δt)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 bb' é a derivada de b(x)b(x) com respeito a xx, e ΔWn=Wτn+1Wτn\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 Δt\Delta t.

Exercício:

Modifique o solver acima implementando o método de Milstein

In [8]:
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 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('<center>Solução pelo método de Milstein <br>da EDE:</center>')) pretty_print(html('<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,0.5,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 ΔX(t)=(ΔS,ΔI)T\Delta X(t)=(\Delta S, \Delta I)^T, onde ΔXi=Xi(t+Δt)Xi(t)\Delta X_i = X_i(t+\Delta t) - X_i(t), i=1,2i=1,2.

Então o valor esperado de ΔX(t)\Delta X(t) na escala detempo de Δt\Delta t é: E(ΔX(t))=(βNSIβNSIγI)Δ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 ΔX(t)\Delta X(t) é V(ΔX(t))E(ΔX(t)[ΔX(t)]TV(\Delta X(t))\approx E(\Delta X(t)[\Delta X(t)]^T, ou seja, V(Δx(t))=(βNSIβNSIβNSIβNSI+γI)Δt.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: {dSdt=βNSI+B11dW1dt+B12dW2dtdIdt=βNSIγI+B21dW1dt+B22dW2dt\{\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} onde B=V(ΔX(t))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(Xt)dta(X_t)dt, e sua covariância, b(Xt)dWb(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 é:

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

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

dP(t)dt=P(t)Q\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)=[X1(t),X2(t),,XN(t)]X(t)=[X_1(t), X_2(t), \ldots, X_N(t)]

Onde Xn(t)X_n(t) é a quantidade do nn-ésimo tipo no tempo tt.

Existem MM 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 ama_m e um vetor de alteração de estados Vm=[Vm1,,VmN]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.

In [9]:
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 # Divided by N (500) because in the stochastic model, # beta is a probability, not a rate 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, Tmed, constant = pars betat = lambda t: beta+(0.5*beta)*np.cos((2*np.pi*t)/365.) 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 if U <= pinf: # Next event is an infection gt = exponential(1/R) 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] + exponential(1/R)) # -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
In [21]:
beta=0.1 gam = 1/21 N = 500 Tmed = 20 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, Tmed, constant) 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')

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

In [54]:
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) betat = lambda t: beta+(0.0*beta)*np.cos((2*np.pi*t)/365.) 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+betat(T)*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 = ((betat(T)/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
In [55]:
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

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.

In [ ]:
multi
In [8]:
!pip3 install --user stochpy
Requirement already satisfied: stochpy in ./.local/lib/python3.6/site-packages (2.3)
In [9]:
import stochpy; import os
####################################################################### # # # Welcome to the interactive StochPy environment # # # ####################################################################### # StochPy: Stochastic modeling in Python # # http://stochpy.sourceforge.net # # Copyright(C) T.R Maarleveld, B.G. Olivier, F.J Bruggeman 2010-2015 # # DOI: 10.1371/journal.pone.0079345 # # Email: [email protected] # # VU University, Amsterdam, Netherlands # # Centrum Wiskunde Informatica, Amsterdam, Netherlands # # StochPy is distributed under the BSD licence. # ####################################################################### Version 2.3.0 Output Directory: /home/user/Stochpy Model Directory: /home/user/Stochpy/pscmodels

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

In [10]:
modelpath = '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 = 999 I = 500 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:

In [11]:
smod = stochpy.SSA() smod.Model('SIR.psc')
In [12]:
# os.system('rm -f '+DATA+'*.dat') smod.DoStochSim(method="Direct", end = int(50),mode = 'time',trajectories = int(10))
[|| ] 10%, Runtime: 0.30 sec [|||| ] 20%, Runtime: 0.42 sec [|||||| ] 30%, Runtime: 0.54 sec [|||||||| ] 40%, Runtime: 0.65 sec [|||||||||| ] 50%, Runtime: 0.75 sec [|||||||||||| ] 60%, Runtime: 0.85 sec [|||||||||||||| ] 70%, Runtime: 0.97 sec [|||||||||||||||| ] 80%, Runtime: 1.07 sec [|||||||||||||||||| ] 90%, Runtime: 1.17 sec
In [13]:
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')
In [14]:
smod.PlotSpeciesTimeSeries(species2plot=["I"]) P.grid()
In [15]:
smod.PlotSpeciesDistributions(linestyle= 'solid')
/home/user/.local/lib/python2.7/site-packages/stochpy/modules/Analysis.py:577: DeprecationWarning: object of type <type 'numpy.float64'> cannot be safely interpreted as an integer. L_bin_edges = np.linspace(dat_min-bin_size/2.0,dat_max+bin_size/2.0,n_bins+1) /ext/sage/sage-8.8_1804/local/lib/python2.7/site-packages/matplotlib/axes/_axes.py:6571: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg. warnings.warn("The 'normed' kwarg is deprecated, and has been "
In [ ]:
In [16]:
smod.PlotSpeciesTimeSeries(species2plot=["S","I","R"])
In [ ]:
In [ ]: