Kernel: SageMath 9.5
In [1]:
%auto typeset_mode(True, display=False)
UsageError: Line magic function `%auto` not found.
In [0]:
%html <h1>Conceitos básicos de processos estocásticos em populações</h1> <p>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.</p> <h2>Processo Poisson simples</h2> <p>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$.</p> <p>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:</p> <p>$(i)\, Pr\{N(t, t+dt)=0\}=1-\alpha dt + o(dt)$,</p> <p>$(ii)\, Pr\{N(t, t+dt)=1\}= \alpha dt +o(dt)$</p> <p>$(iii)\, Pr\{N(t, t+dt)>1\}= o(dt)$</p> <p>$(iv)\, N(t, t+dt)$ é independente de todas as entradas no intervalo $(0, t)$</p> <p>Aqui, o(dt) denota todos os termos de order inferior a dt.</p> <h3>Tempo até o próximo evento</h3> <p>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$.</p> <p>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_o+t, t_0+t+dt)$. Podemos ver que</p> <p>$G(t+dt)=Pr\{Z>t+dt\}$</p> <p>$=Pr\{Z>t$ e nenhum evento ocorre no intervalo $(t_0+t, t_0+t+dt)\}$ </p> <p>$=Pr\{Z>t\} Pr\{$nenhum evento ocorre no intervalo $(t_0+t, t_0+t+dt)\mid Z>t\}$ </p> <p>aplicando a condição $(iv)$ acima, </p> <p>$G(t+dt) =Pr\{$nenhum evento ocorre no intervalo $(t_0+t, t_0+t+dt)}$</p> <p>Aplicando a condição $(i)$ obtemos</p> <p>$G(t+dt)=G(t)(1-\alpha dt + o(dt)$</p> <p>Quando $dt$ tende a $0$, temos</p> <p>$\frac{dG(t)}{dt}= -\alpha G(t)$</p> <p>Cuja solução é</p> <p>$G(t)=G(0) e^{-\alpha t}$</p> <p>Como $G(0)=Pr(Z>0)=1$ concluímos que a função de distribuição acumulada de $Z$ é</p> <p>$F(t) = Pr(Z \leq t)=1-G(t)= 1- e^{-\alpha t}$</p> <p>Logo a pdf do tempo até o próximo evento é: </p> <p>$f(t)=F'(t)=\alpha e^{-\alpha t} \,\, (\alpha > 0)$</p> <p>que corresponde à distribuição exponencial com parâmetror $\alpha$.</p> <h2>Exercício: Simulando o processo de Poisson</h2> <p>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$). </p>
In [0]:
import numpy as np from numpy.random import exponential as expv import pylab as P from numpy import cumsum alpha = 1.5 times = cumsum(expv(alpha, size=100)) list_plot(zip(times, np.ones(100)))
In [0]:
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.clf() 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() P.savefig('hist.png')
In [0]:
%html <h1>Equações Diferenciais Estocásticas (EDEs)</h1> <h4>Copyright (2013-15) Flávio Codeço Coelho, EMAp-FGV</h4> <p> 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:</p> <p>$$\mathrm{d}X_t = a(X_t) \, \mathrm{d} t + b(X_t) \, \mathrm{d} W_t,$$</p> <p> </p> <p>onde $W_t$ é um processo de Wiener, ou seja, um processo estocástico em tempo contínuo, também conhecido como movimento browniano.</p> <p>A solução analítica destas equações foge ao escopo deste curso. Mas podemos explorar soluções numéricas.</p> <h2>Modelando com EDEs</h2> <p>Vamos agora desenvolver alguns modelos simples usando EDEs.</p> <h3>Dinâmica de uma população (Processo de nascimentos e mortes)</h3> <p>Seja o seguinte modelo da dinâmica demográfica de uma população:</p> <p>$$\frac{\mathrm{d}x}{\mathrm{d}t}=b x - d x; \,\,\,\,\,\, x(0)=x_0$$</p> <p>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 é:</p> <p>$$x(t)= x_0e^{(b-d)t}$$</p> <p> </p>
In [0]:
%auto var('t') b = 0.3 d = 0.27 x0 = 25 x(t) = x0*exp((b-d)*t) plot(x(t),(t,0,50))
In [0]:
%html <p>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.</p> <p>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.</p> <p> </p> <p>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$.</p> <p>Neste intervalo, são possíveis três tipos de alterações na população:</p> <ol> <li>Um nascimento, $\beta_1=1$,</li> <li>uma morte, $\beta_2=-1$,</li> <li>ou nada, $\beta_3=0$.</li> </ol> <p>Estas alterações ocorrerão com as respectivas probabilidades:</p> <ol> <li>$P_1 = bx\Delta t$,</li> <li>$P_2 = dx\Delta t$ e</li> <li>$P_3 = 1-(b+d)x\Delta t$.</li> </ol> <p>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$.</p> <p>Isto nos leva a:</p> <p>$$p_k(t+\Delta t) = p_k(t)[1 - bx\Delta t - dx\Delta t]+p_{k-1}(t)[bx_{k-1}\Delta t]+p_{k+1}(t)[dx_{k+1}\Delta t]$$</p> <p>Que pode ser reescrita como:</p> <p>$$\frac{p_k(t+\Delta t)-p_k(t)}{\Delta t} = \frac{(b-d)(x_{k+1}p_{k+1}(t) - x_{k-1}p_{k-1}(t))}{2\Delta x} + \frac{1}{2}\frac{(b+d)(x_{k+1}p_{k+1}(t) - 2x_kp_k(t) + x_{k-1}p_{k-1}(t))}{(\Delta x)^2}$$</p> <p>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.</p> <p>$$\mathrm{d}x(t) = (b - d)x \mathrm{d}t + \sqrt{((b+d)x)}\mathrm{d}W(t)$$</p> <p>onde $W(t)$ é um processo de Wiener.</p> <p>Para resolver numericamente esta EDE precisamos entender o método de Euler-Maruyama.</p> <h2>Método de Euler-Maruyama</h2> <p>Considere a seguinte equação diferencial estocástica de <a title="Itō calculus" href="http://en.wikipedia.org/wiki/It%C5%8D_calculus">Itō</a> :</p> <p>$$\mathrm{d}X_t = a(X_t) \, \mathrm{d} t + b(X_t) \, \mathrm{d} W_t,$$</p> <p>Com <a href="http://en.wikipedia.org/wiki/Initial_condition" target="_blank">condições iniciais</a> $X_0=x_0$, onde where $W_t$<em> é o <a href="http://en.wikipedia.org/wiki/Wiener_process" target="_blank">processo de Wiener</a></em>, e suponha que queremos resolver esta EDE em um dado intervalo $[0, T]$. Então a aproximação de <strong>Euler–Maruyama approximation</strong> à solução verdadeira de $X$ é a cadeia de Markov $Y$ definida abaixo:</p> <ul> <li>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;$$</li> <li>faça $Y_0=x_0;$</li> <li>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$</li> </ul> <h3>Solução numérica do processo de nascimento e morte</h3> <p>Vamos agora definir as funções que correspondem às partes determinísticas e estocásticas da nossa EDE.</p>
In [0]:
%auto det = lambda x,b,d:(b-d)*x stoc = lambda x,b,d: sqrt((b+d)*x)
In [0]:
%html <p>Agora definimos uma função que implementa o método de Euler-Maruyama, Junto com uma figura interativa aplicando-o ao nosso modelo.</p>
In [0]:
%auto def EulerMaruyama(tstart, ystart, tfinish, nsteps, f1, f2, params): sol = [ystart] tvals = [tstart] h = N((tfinish-tstart)/nsteps) for step in range(nsteps): sol.append(sol[-1] + h*f1(sol[-1],*params) + f2(sol[-1],*params)*normalvariate(0,h)) tvals.append(tvals[-1] + 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,det,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 Não estocástica') for i in range(1,plots_at_a_time): sol2 = EulerMaruyama(0,25,50,100,det,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])
In [0]:
%html <h2>O Método de Milstein</h2> <p>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:</p> <p>$$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)$$</p> <p>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$.</p> <h4>Exercício:</h4> <p>Modifique o solver acima implementando o <a href="https://en.wikipedia.org/wiki/Milstein_method">método de Milstein</a>. </p>
In [0]:
%html <table border="0" rules="cols"> <tbody> <tr> <td> </td> <td> </td> </tr> <tr> <td> </td> <td> </td> </tr> </tbody> </table> <h2>Derivando uma EDE para duas populações:</h2> <p>seja $\Delta X=(\Delta X_1, \Delta X_2)^T$, onde $\Delta X_i = X_i(t+\Delta t) - X_i(t)$, $i=1,2$.</p> <p>Seja as possíveis transições nas duas populações de finidas na tabela abaixo:</p> <table border="0" align="center"> <thead> <tr> <td>$i$ </td> <td>Mudança $(\Delta X)_i$ </td> <td>Probabilidade, $p_i$ </td> </tr> </thead> <tbody> <tr> <td>1</td> <td> $(1,0)^T$</td> <td> $b_1 \Delta t$</td> </tr> <tr> <td>2 </td> <td> $(0,1)^T$ </td> <td> $b_2 \Delta t$</td> </tr> <tr> <td>3 </td> <td> </td> <td> </td> </tr> <tr> <td>4 </td> <td> </td> <td> </td> </tr> <tr> <td>5 </td> <td> </td> <td> </td> </tr> <tr> <td>6 </td> <td> </td> <td> </td> </tr> </tbody> </table> <h1>Equações de Kolmogorov</h1> <p>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 é:</p> <p>$$\frac{dP(t)}{dt} QP(t),$$</p> <p>onde $P(t)=(p_{ji}(t)$ é a <em>matriz de probabilidades de transição</em> e $Q=(q_{ji})$, é a <em>matriz geradora</em>. De forma similar, a versão retrógrada é dada por:</p> <p>$$\frac{dP(t)}{dt} P(t)Q$$</p> <p> </p> <p>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 "<a href="http://en.wikipedia.org/wiki/Gillespie_algorithm" target="_blank">direto</a>" 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á.</p> <p>Para descrever esta classe de modelos, é interessante pensar no modelo como um vetor de estados</p> <p>$$X(t)=[X_1(t), X_2(t), \ldots, X_N(t)]$$</p> <p>Onde $X_n(t)$ é a quantidade do $n$-ésimo tipo no tempo $t$.</p> <p>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} ]$.</p> <p> </p> <h2>Simulando modelos estocásticos em Python</h2> <p>para simular tais modelos faremos uso de 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.</p>
In [0]:
import stochpy; import pylab as plt import os
In [0]:
%html <p>Para escrever um modelo no Pysces, precisamos utilizar a <a href="http://stochpy.sourceforge.net/html/inputfile_doc.html#pysces-inputfile" target="_blank">Pysces Modeling language</a>, vejamos como fica um simples modelo SIR</p>
In [0]:
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(DATA+'SIR.psc','w') as f: f.write(sir)
In [0]:
%html <p>Agora vamos definir o nosso método de solução numérica, e carregar o modelo:</p>
In [0]:
smod = stochpy.SSA() smod.Model(DATA+'SIR.psc')
In [0]:
%python os.system('rm -f '+DATA+'*.dat') smod.DoStochSim(method="Direct", end = 50,mode = 'time',trajectories = 10)
In [0]:
%python plt.hold(True) plt.clf() plt.grid() smod.PlotAverageSpeciesTimeSeries(species2plot=["S","I","R"]) plt.grid() plt.hold(False) plt.savefig('sir.png')
In [0]:
smod.PlotSpeciesTimeSeries(species2plot=["I"]) plt.grid() plt.savefig('i.png')
In [0]:
plt.cla() smod.PlotSpeciesDistributions(linestyle= 'solid') plt.savefig('sir_dist.png')
In [0]:
plt.cla() smod.PlotWaitingtimesDistributions() plt.savefig('sir_wait.png')
In [0]:
plt.cla() smod.PlotSpeciesTimeSeries(species2plot=["S","I","R"]) plt.savefig('sir_ts.png')