Simulação Estocástica

Simulação Estocástica Thiago Rodrigo Ramos 2 Sumário 1 Elementos básicos de probabilidade 9 1.1 Axiomas da probabilidade .................................. 9 1.1.1 Probabilidade condicional e independência .................... 9 1.2 Variáveis aleatórias ...................................... 10 1.3 Valor esperado ......................................... 11 1.4 Variância ............................................ 12 1.4.1 Covariância ...................................... 12 1.5 Desigualdades básicas de concentração .......................... 14 1.6 Teoremas assintóticos ..................................... 14 2 Variáveis discretas e como simulá-las 17 2.1 Variáveis com suporte finito ................................. 19 2.2 Bernoulli ............................................ 20 2.3 Distribuição binomial ..................................... 21 2.3.1 Simulando via Bernoullis .............................. 21 2.3.2 Simulando via identidade recursiva ........................ 22 2.3.3 Aspectos computacionais .............................. 23 2.3.4 Número médio de passos em algoritmos de inversão recursiva ........ 24 2.4 Distribuição geométrica ................................... 25 2.4.1 Simulando via Bernoullis .............................. 25 2.4.2 Simulando a geométrica via inversão ....................... 26 2.5 Distribuição de Poisson ................................... 27 2.5.1 Simulação a Poisson via inversão e recursão ................... 28 2.5.2 Algoritmo melhorado ................................ 29 2.5.3 Relação com a binomial ............................... 30 2.6 Distribuição binomial negativa ............................... 31 2.6.1 Simulando via Bernoullis .............................. 32 2.6.2 Simulando via soma de geométricas ........................ 32 2.6.3 Simulando via inversão recursiva ......................... 33 2.6.4 Por que o nome “Binomial Negativa”? ...................... 34 2.7 Distribuição hipergeométrica ................................ 34 2.7.1 Simulando a Hipergeométrica ........................... 35 3 4SUMÁRIO 3 Variáveis contínuas e como simulá-las 37 3.1 Método da Inversão ...................................... 37 3.1.1 Distribuição exponencial .............................. 38 3.2 Método da rejeição-aceitação ................................ 41 3.2.1 Distribuição normal ................................. 44 3.3 Distribuição Gamma ..................................... 47 3.3.1 Simulando quando αé inteiro ........................... 48 3.3.2 Simulando quando α>1 via aceitação–rejeição com Exponencial ....... 48 3.3.3 Simulando quando α<1.............................. 50 3.4 Distribuição Beta ....................................... 51 3.4.1 Simulando a Beta via aceitação–rejeição com proposta uniforme ....... 52 3.4.2 Simulando a Beta via Gammas independentes .................. 53 3.5 Transformações de Variáveis Aleatórias .......................... 53 3.5.1 Geração de Normais via Método de Box–Muller ................. 54 3.5.2 Geração da normal bivariada ............................ 56 3.5.3 Distribuição Qui-quadrado ............................. 59 3.5.4 Simulando a distribuição t de Student ....................... 61 4 Simulação via Monte Carlo 63 4.1 Estimando médias ...................................... 63 4.1.1 Exemplos ....................................... 63 4.2 Intervalos de Confiança ................................... 66 5 Redução de variância 69 5.1 Uso de variáveis antitéticas ................................. 69 5.2 O uso de variáveis de controle ............................... 74 5.3 Redução de Variância por Condicionamento ....................... 77 6 Amostragem por importância 81 6.1 Densidades Inclinadas (Tilted Densities) .......................... 82 6.2 Desigualdade de Chernoff .................................. 85 6.3 Variância sob Inclinação Exponencial ........................... 86 6.4 Variáveis Sub-Gaussianas e Desigualdade de Hoeffding ................ 88 6.5 Por que a Inclinação Exponencial? ............................. 89 7 Cadeias de Markov e MCMC 93 7.1 Cadeias de Markov (Resumo) ................................ 94 7.1.1 Classificação dos estados .............................. 97 7.1.2 Distribuição estacionária ..............................100 7.1.3 Reversibilidade ....................................104 7.2 Markov Chain Monte Carlo .................................107 7.3 Algoritmo de Metropolis–Hastings .............................108 7.4 Amostragem de Gibbs ....................................113 SUMÁRIO 5 8 Processos de Difusão 119 8.1 Movimento Browniano e SDE ................................119 8.2 Equação de Fokker–Planck .................................122 8.3 Amostragem de Langevin ..................................123 8.4 Denoising Score Matching ..................................125 8.4.1 Estimando o score na prática ............................127 8.4.2 Etapa de Inferência via Langevin ..........................129 9 Bootstrap 133 9.1 Uma visão pragmática de Bootstrap ............................133 9.2 Uma visão teórica de Bootstrap ...............................136 9.2.1 A desigualdade de Dvoretzky–Kiefer–Wolfowitz ................136 9.2.2 Bootstrap e pontes brownianas ...........................139 10 Estratégias para acelerar códigos em Python 145 10.1 Profiling com cProfile ...................................145 10.2 Paralelização com joblib.Parallel ............................148 10.3 Compilação Just-In-Time com Numba ...........................150 10.4 Paralelismo simples em Bash ................................153 6SUMÁRIO IMPORTANTE Estas notas de aula ainda estão em construção. Diversas partes do texto encontram-se em revisão e, em particular, as referências bibliográficas aos artigos e livros utilizados em sua elaboração ainda serão incluídas nas próximas versões. O conteúdo atual deve, portanto, ser considerado preliminar. Última atualização: 2025-11-12 SUMÁRIO 7 Um conselho: a importância de ser ruim antes de ser bom É natural que, quando começamos a fazer algo, a gente faça essa coisa muito malfeita ou cheia de defeitos. Isso é comum em qualquer processo de aprendizagem, e sempre foi assim, desde o início dos tempos. Quando comecei a programar em Python, muita coisa sobre a linguagem eu aprendi por conta própria, apesar de já ter feito alguns cursos básicos em C. Programei de forma amadora em Python por muitos anos, até que, no doutorado, precisei aprender a programar de forma mais organizada e profissional. Lembro que, nessa época, um amigo da pós-graduação me apresentou ao "submundo da programação". Foi aí que aprendi muito do que sei hoje sobre terminal do Linux, Git, e foi também quando comecei a usar o Vim. Uma das coisas que esse amigo me mostrou foi o Pylint, que nada mais é do que um verificador de bugs e qualidade de código para Python. O Pylint é bem rigoroso na análise, e ainda te dá, ao final, uma nota que vai até 10. Nessa fase, apesar de já ter evoluído bastante, meus códigos ainda recebiam notas por volta de 6 ou 7. Resolvi então rodar o Pylint nos meus códigos antigos pra ter uma noção de quão ruins eles eram — e a nota final foi -900. Pois é, existe um limite superior para o quão bem você consegue fazer algo, mas aparentemente o fundo do poço é infinito. O que eu queria mostrar com essa história é que faz parte do processo de aprendizado ser ruim no começo e melhorar com o tempo. Falo isso porque, hoje em dia, com o crescimento dos LLMs, a gente fica tentado a pular essa etapa de errar muito até acertar, e ir direto pra fase em que escrevemos códigos limpos, bem comentados, identados e organizados. Mas não se enganem: apesar da aparência profissional, depender de LLMs pra escrever tudo atrapalha justamente essa parte essencial de aprender errando. Neste curso, vários exercícios envolvem escrever códigos em Python. Meu conselho é: não tenham vergonha de errar, de escrever soluções ruins ou confusas. Isso é absolutamente normal. Vocês estão aqui para evoluir — e errar faz parte do processo. 8SUMÁRIO Capítulo 1 Elementos básicos de probabilidade 1.1 Axiomas da probabilidade Um espaço de probabilidade é uma tupla composta por três elementos: o espaço amostral, o conjunto de eventos e uma distribuição de probabilidade: •Espaço amostral Ω:Ωé o conjunto de todos os eventos elementares ou resultados possíveis de um experimento. Por exemplo, ao lançar um dado, Ω={1,2,3,4,5,6}. •Conjunto de eventos F:Fé uma σ-álgebra, ou seja, um conjunto de subconjuntos de Ω que contém Ωe é fechado sob complementação e união enumerável (e, consequentemente, também sob interseção enumerável). Um exemplo de evento é: “o dado mostra um número ímpar”. •Distribuição de probabilidade P:Pé uma função que associa a cada evento de Fum número em [0,1], tal que P[Ω] = 1, P[∅] = 0 e, para eventos mutuamente exclusivos A1, . . . , An, temos: P[A1∪···∪ An]= n ∑ i=1 P[Ai]. A distribuição de probabilidade discreta associada ao lançamento de um dado justo pode ser definida como P[Ai] = 1/6 para i∈ {1, . . . , 6}, onde Aié o evento “o dado mostra o valor i”. 1.1.1 Probabilidade condicional e independência A probabilidade condicional do evento Adado o evento Bé definida como a razão entre a probabilidade da interseção A∩Be a probabilidade de B, desde que P[B]=0: P[A|B] = P[A∩B] P[B]. Dois eventos AeBsão ditos independentes quando a probabilidade conjunta P[A∩B]pode ser fatorada como o produto P[A]P[B]: P[A∩B] = P[A]P[B]. 9 16 CAPÍTULO 1. ELEMENTOS BÁSICOS DE PROBABILIDADE Portanto, ∆k=EhfD(k) n+Xk √n−fD(k) n+Nk √ni. Passo 3: Expansão de Taylor condicional. Fixe D(k) n=d. Aplicando Taylor em torno de d, temos: fd+Xk √n=f(d) + Xk √nf′(d) + X2 k 2nf′′(d) + X3 k 6n3/2 f(3)(d+ξX), fd+Nk √n=f(d) + Nk √nf′(d) + N2 k 2nf′′(d) + N3 k 6n3/2 f(3)(d+ξN), para alguns ξX,ξNentre 0 e Xk/√nou Nk/√n. Subtraindo, fd+Xk √n−fd+Nk √n=1 √nf′(d)(Xk−Nk) + 1 2nf′′(d)(X2 k−N2 k) + Rk(d), onde Rk(d) = 1 6n3/2 X3 kf(3)(d+ξX)−N3 kf(3)(d+ξN). Passo 4: Tomando esperança condicional. Voltamos para ∆k=EhfD(k) n+Xk √n−fD(k) n+Nk √ni. Usando a decomposição anterior e condicionando em D(k) n, temos: ∆k=Eh1 √nf′(D(k) n)(Xk−Nk)i +Eh1 2nf′′(D(k) n)(X2 k−N2 k)i +E[Rk]. Agora, como XkeNksão independentes de D(k) n, obtemos: E[f′(D(k) n)(Xk−Nk)] = E[f′(D(k) n)] ·(E[Xk]−E[Nk]) = 0, E[f′′(D(k) n)(X2 k−N2 k)] = E[f′′(D(k) n)] ·(E[X2 k]−E[N2 k]) = 0. Portanto, só resta ∆k=E[Rk]. Passo 5: Controle do resto. Do termo Rk, temos |Rk| ≤ 1 6n3/2 |Xk|3sup |f(3)|+|Nk|3sup |f(3)|. Tomando esperança, |E[Rk]| ≤ C n3/2 E[|X1|3] + E[|N1|3], onde C=1 6sup |f(3)|. Somando sobre k,  n ∑ k=1 ∆k≤n·C n3/2 E[|X1|3] + E[|N1|3]=O1 √n→0. Logo, E[f(An)] −E[f(Bn)] →0, e como Bn∼N(0,1)para todo n, segue que And →N(0,1). Capítulo 2 Variáveis discretas e como simulá-las O ponto de partida do nosso curso será sempre o mesmo: só podemos utilizar variáveis uniformes para gerar todas as demais distribuições. Ou seja, assumimos que temos disponível uma variável aleatória U∼Uniforme(0,1), e a partir dela construiremos algoritmos para simular outras variáveis. A propriedade fundamental dessa variável é: P(a<U<b) = b−a, 0 ≤a<b≤1. Isto é, a probabilidade de Ucair em um subintervalo do intervalo (0,1)é igual ao comprimento desse subintervalo. Exercício 9. Seja U ∼Uniforme(0,1). Mostre que, para quaisquer números 0≤a<b≤1, P(a<U<b) = b−a. Para variáveis discretas, essa ideia pode ser usada da seguinte forma: suponha que Xassuma valores x1,x2, . . . , xmcom probabilidades p1,p2, . . . , pm, onde pk=P(X=xk),pk≥0, m ∑ k=1 pk=1. Definimos as probabilidades acumuladas Fk= k ∑ i=1 pi,k=1, . . . , m. Então, o algoritmo de simulação é: 1. Gerar U∼Uniforme(0,1); 2. Encontrar o menor índice ktal que U≤Fk; 3. Retornar X=xk. 17 18 CAPÍTULO 2. VARIÁVEIS DISCRETAS E COMO SIMULÁ-LAS A propriedade P(a<U<b) = b−agarante que P(X=xk) = pk. De forma intuitiva, dividimos o intervalo (0,1)em subintervalos consecutivos de comprimentos pk. Ao sortearmos U∼Uniforme(0,1), o valor de Xserá aquele correspondente ao subintervalo no qual Ucair. Esse procedimento é conhecido como método da inversão para variáveis discretas. A Figura 2.1 ilustra esse processo para uma variável Bernoulli. Figura 2.1: Particionamento do intervalo (0,1)para simular uma variável Bernoulli com p=0.7. Sorteia-se U∼Uniforme(0,1); se Ucair na região azul, definimos X=0, e caso contrário, X=1. A mesma ideia se aplica quando o conjunto de valores possíveis de Xé infinito (ou muito grande). Nesse caso, o intervalo (0,1)é particionado em uma sequência de subintervalos, cada um correspondente a um valor de X, como ilustrado na Figura 2.2. Figura 2.2: Particionamento do intervalo (0,1)para simular uma variável discreta com suporte infinito. O nome método da inversão vem do fato de que a simulação utiliza a função de distribuição acumulada (CDF) e sua inversa generalizada. Seja Xuma variável aleatória com CDF F(x). Então, se U∼Uniforme(0,1), vale que X=F−1(U), onde a inversa generalizada é definida por F−1(u) = min{x:F(x)≥u}, 0 <u<1. No caso discreto, isto corresponde exatamente ao passo do algoritmo em que escolhemos o menor ktal que U≤Fk. Ou seja, sorteamos U, e depois “invertemos” a CDF para recuperar uma realização de Xna sua escala original. 2.1. VARIÁVEIS COM SUPORTE FINITO 19 Esse procedimento pode parecer um pouco abstrato neste momento, já que a noção de inversa de uma função acumulada fica mais clara quando lidamos com variáveis contínuas. Por isso, retornaremos a esse método mais adiante, ao estudarmos a simulação de variáveis contínuas via inversão. Antes, porém, vale formalizar essa ideia de maneira geral. Exercício 10. Seja X uma variável aleatória com função de distribuição acumulada FX. Considere U ∼ Uniforme(0,1)e defina Y=F−1 X(U),onde F−1 X(u) = min{x:FX(x)≥u}. Prove que Y tem a mesma distribuição que X. Esse resultado mostra que, a partir de uma variável uniforme, podemos simular qualquer outra distribuição usando a CDF e sua inversa generalizada. Com essa ferramenta em mãos, passamos agora ao estudo de algumas distribuições discretas fundamentais, que servirão de exemplo concreto dessa ideia. 2.1 Variáveis com suporte finito Comecemos com o caso em que Xassume um número finito de valores x1,x2, . . . , xm, cada um com probabilidade pj=P(X=xj). Por exemplo, suponha que p1=0.20, p2=0.15, p3=0.25, p4=0.40. Uma maneira direta de simular Xé gerar U∼Uniforme(0,1)e aplicar: • Se U<0.20, definir X=1 e pare; • Se U<0.35, definir X=2 e pare; • Se U<0.60, definir X=3 e pare; • Caso contrário, definir X=4. Embora possamos reordenar os testes para tornar a verificação mais eficiente, a ideia central permanece a mesma: dividir o intervalo (0, 1)em partes de comprimentos pje identificar onde Ucaiu. 20 CAPÍTULO 2. VARIÁVEIS DISCRETAS E COMO SIMULÁ-LAS De forma geral, se Xé uma variável com suporte finito S={x1,x2, . . . , xm}, sua distribuição é completamente determinada pela função de probabilidade pX(xk) = P(X=xk),xk∈S, a qual satisfaz pX(k)≥0 para todo k∈S,∑ k∈S pX(k) = 1. Exemplo 2. Seja S ={x1,x2, . . . , xK}um conjunto de K valores distintos. Dizemos que X tem distribuição uniforme discreta em S quando pX(xi) = 1 K,i=1, 2, . . . , K. Nesse caso, cada valor é igualmente provável e temos K ∑ i=1 pX(xi) = K ∑ i=1 1 K=1. Um caso especial é a uniforme discreta nos inteiros 1, 2, . . . , n, em que P(X=j) = 1 n,j=1, . . . , n. Neste cenário, o método se torna extremamente simples: basta gerar U∼Uniforme(0,1)e definir X=⌊nU⌋+1, onde ⌊x⌋indica a parte inteira de x(maior inteiro menor ou igual a x). De fato, X=jse e somente se j−1≤nU <j, o que ocorre com probabilidade 1 n. Variáveis uniformes discretas são particularmente importantes em simulação, pois permitem gerar inteiros equiprováveis de forma extremamente eficiente. 2.2 Bernoulli A distribuição de Bernoulli modela experimentos com dois resultados possíveis, tipicamente denominados “sucesso” (valor 1) e “fracasso” (valor 0). Dizemos que X∼Bernoulli(p)se P(X=1) = peP(X=0) = 1−p, 2.3. DISTRIBUIÇÃO BINOMIAL 21 onde 0 ≤p≤1 representa a probabilidade de sucesso. A função de probabilidade (pmf) pode ser escrita de forma compacta como pX(k) = pk(1−p)1−k,k∈ {0,1}. As principais características dessa distribuição são: E[X] = p, Var(X) = p(1−p). Exercício 11. Prove as propriedades acima, isto é, calcule a esperança e a variância de uma variável Bernoulli. No contexto de simulação, a Bernoulli é um caso particular da uniforme discreta em {0,1} com probabilidades 1−pep, respectivamente. O algoritmo é simples: sorteamos U∼Uniforme(0,1) e definimos X=   1, se U≤p, 0, caso contrário. Exercício 12. Mostre que o procedimento acima gera corretamente uma variável Bernoulli, isto é, verifique que P(X=1) = p e P(X=0) = 1−p. 2.3 Distribuição binomial A distribuição binomial modela o número de sucessos em nrepetições independentes de um experimento de Bernoulli com probabilidade de sucesso p. Sejam X1,X2, . . . , Xnvariáveis aleatórias independentes, todas com distribuição Bernoulli(p). Definimos X= n ∑ i=1 Xi. Nesse caso, dizemos que X∼Binomial(n,p), cuja função de probabilidade é P(X=k) = n kpk(1−p)n−k,k=0, 1, . . . , n. As principais propriedades são: E[X] = np, Var(X) = np(1−p). Exercício 13. Prove as propriedades acima. 2.3.1 Simulando via Bernoullis Uma forma simples e direta de simular uma variável aleatória binomial é a partir de variáveis de Bernoulli independentes. Recorde que se X∼Binomial(n,p), então Xpode ser escrito como X= n ∑ i=1 Bi, 22 CAPÍTULO 2. VARIÁVEIS DISCRETAS E COMO SIMULÁ-LAS onde B1,B2, . . . , Bnsão variáveis independentes e identicamente distribuídas, cada uma com Bi∼Bernoulli(p). Assim, o algoritmo de simulação da binomial segue naturalmente: 1. Para i=1, . . . , n, gerar Bi∼Bernoulli(p); 2. Retornar X=∑n i=1Bi. Em outras palavras, uma variável binomial conta o número de sucessos em ntentativas independentes, cada uma com probabilidade de sucesso p. Portanto, simular uma binomial se reduz a repetir nvezes o procedimento de simulação da Bernoulli e somar os resultados. 2.3.2 Simulando via identidade recursiva Uma alternativa mais eficiente utiliza o método da inversão, aproveitando a identidade recursiva da função massa de probabilidade da Binomial. Se X∼Binomial(n,p), então P(X=i) = n ipi(1−p)n−i,i=0, 1, . . . , n. Essas probabilidades satisfazem uma relação de recorrência simples. De fato, começando em P(X=i+1) = n i+1pi+1(1−p)n−i−1, observamos que 1 n i+1=n! (i+1)!(n−i−1)!=n−i i+1 n! i!(n−i)!=n−i i+1n i. Substituindo essa relação, P(X=i+1) = n−i i+1n ipi+1(1−p)n−i−1. Reorganizando, P(X=i+1) = n−i i+1·p 1−pP(X=i). Assim, conhecendo P(X=0) = (1−p)n, podemos calcular P(X=1),P(X=2), . . . de forma recursiva, sem reavaliar coeficientes binomiais nem potências. Isso leva ao seguinte algoritmo de simulação via inversão: 1. Gerar U∼Uniforme(0,1); 1Por exemplo, se n=10, i=6 e i+1=7, então 10! 7! 3! =10! ·4 7·6! ·4·3! =10! 6! ·4! 4 7 2.3. DISTRIBUIÇÃO BINOMIAL 23 2. Inicializar o índice i=0, a probabilidade atual pi= (1−p)ne a soma acumulada F=pi; 3. Enquanto U>F, atualizar pi+1=n−i i+1·p 1−ppi,i←i+1, F←F+pi; 4. Retornar X=i. Esse procedimento verifica primeiro se X=0, depois se X=1, e assim por diante, até encontrar o valor de Xsorteado. Em média, o número de passos necessários é aproximadamente 1+np, o que pode ser bem mais eficiente do que gerar nvariáveis de Bernoulli quando né grande. Exemplo 3. Considere n =5e p =0.3. Temos P(X=0)=(1−0.3)5=0.16807. Suponha que geramos U =0.4. Como U >0.16807, passamos ao próximo valor: p1=5−0 1·0.3 0.7 ·0.16807 ≈0.36015, F=0.16807 +0.36015 =0.52822. Agora U =0.4 <F, logo o algoritmo retorna X =1. Portanto, neste caso específico, o sorteio resultou em exatamente um sucesso entre as cinco tentativas. 2.3.3 Aspectos computacionais A escolha do método para simular variáveis binomiais tem implicações diretas em termos de eficiência. Dois fatores fundamentais influenciam o desempenho: o número de tentativas ne a probabilidade de sucesso p. No método da soma de Bernoullis, o custo de cada amostra é proporcional a n, já que é necessário realizar nsorteios independentes. Esse custo não depende do valor de p: tanto para valores pequenos quanto grandes de p, o algoritmo precisa sempre gerar todas as nBernoullis. Já no método da inversão recursiva, o número médio de passos é da ordem de 1 +np, pois o procedimento acumula probabilidades até ultrapassar o valor sorteado U. Quando pé pequeno, 24 CAPÍTULO 2. VARIÁVEIS DISCRETAS E COMO SIMULÁ-LAS o valor típico da variável Xtambém é pequeno, e o algoritmo tende a parar cedo, podendo ser competitivo em relação à soma de Bernoullis. Por outro lado, quando pé moderado ou grande, o valor esperado np cresce e, com ele, o número de passos, tornando a inversão significativamente mais lenta. Figura 2.3: Comparação de tempo de execução (em segundos) entre o método da inversão recursiva e a soma de Bernoullis para n=100 e N=2000 amostras, variando p. A Figura 2.3 ilustra essa comparação em implementações com loops explícitos, para n=100 e diferentes valores de p. Enquanto o tempo da soma de Bernoullis cresce linearmente apenas com ne não é afetado por p, o tempo do método da inversão cresce proporcionalmente a np, aumentando de forma acentuada à medida que pse aproxima de 1. Na prática, bibliotecas como NumPy utilizam algoritmos especializados para a binomial, ainda mais rápidos do que ambos os métodos discutidos aqui, de modo que a utilidade principal desses algoritmos é didática e comparativa, permitindo compreender os diferentes custos computacionais associados a cada abordagem. 2.3.4 Número médio de passos em algoritmos de inversão recursiva Nos algoritmos recursivos de inversão, a lógica é sempre a mesma: dado um número aleatório U∼Uniforme(0,1), acumulamos as probabilidades da distribuição até que a soma ultrapasse U. O valor de Xsorteado é exatamente o índice kem que essa condição se verifica pela primeira vez. Assim, se o valor sorteado é X=k, o algoritmo precisou verificar todos os valores 0,1,2,. . . , k− 1 e só então aceitou k. Isso significa que o número total de passos é S=k+1. Como Xé a variável aleatória que estamos simulando, temos E[S] = E[X+1] = E[X] + 1. 2.4. DISTRIBUIÇÃO GEOMÉTRICA 25 Esse resultado é geral para qualquer algoritmo de inversão recursiva que inicie a busca no valor mínimo do suporte e avance de forma sequencial. No caso da binomial X∼Bin(n,p), por exemplo, o número esperado de passos é E[S] = 1+np, uma vez que E[X] = np. Portanto, o custo médio do algoritmo está diretamente ligado ao valor esperado da distribuição sorteada: distribuições concentradas em valores pequenos produzem simulações muito rápidas, enquanto distribuições centradas em valores grandes exigem proporcionalmente mais passos. 2.4 Distribuição geométrica A distribuição geométrica modela o número de ensaios de Bernoulli até a ocorrência do primeiro sucesso. Seja pa probabilidade de sucesso em cada tentativa, com 0 <p≤1. Definimos Xcomo o número de ensaios necessários até o primeiro sucesso. Dizemos que X∼Geom(p)se P(X=k) = (1−p)k−1p,k=1, 2,3,. . . Nesse caso: E[X] = 1 p, Var(X) = 1−p p2. Exercício 14. Prove que a função de probabilidade acima satisfaz ∑∞ k=1P(X=k) = 1. Exercício 15. Prove as propriedades acima. 2.4.1 Simulando via Bernoullis A distribuição geométrica modela o número de tentativas até a ocorrência do primeiro sucesso, em uma sequência de experimentos de Bernoulli independentes com probabilidade p∈(0,1)de sucesso. Essa definição leva naturalmente a um algoritmo de simulação baseado em Bernoullis. A ideia é repetir experimentos de Bernoulli até obter sucesso pela primeira vez: 1. Inicializar o contador X←1; 2. Gerar B∼Bernoulli(p); 3. Enquanto B=0, repetir: •X←X+1; • Gerar novo B∼Bernoulli(p); 4. Retornar X. Note que esse procedimento reflete exatamente a definição da variável: contar quantas tentativas são necessárias até que ocorra o primeiro sucesso. 32 CAPÍTULO 2. VARIÁVEIS DISCRETAS E COMO SIMULÁ-LAS Exercício 18. Prove que a função de probabilidade acima é válida, isto é, que ∑∞ n=rP(X=n) = 1. Exercício 19. Prove as fórmulas da média e variância usando o fato de que X é a soma de r variáveis independentes com distribuição Geom(p). 2.6.1 Simulando via Bernoullis Uma forma direta de gerar uma variável Binomial Negativa é simular sucessivos ensaios de Bernoulli(p) até obter o r-ésimo sucesso. De fato, por definição, Xrepresenta o número total de ensaios necessários até a ocorrência de rsucessos. Assim, o algoritmo pode ser descrito da seguinte forma: 1. Inicializar n=0 (contador de ensaios) e s=0 (contador de sucessos); 2. Enquanto s<r: (a) Gerar B∼Bernoulli(p); (b) Atualizar n←n+1; (c) Se B=1, atualizar s←s+1; 3. Retornar X=n. Esse método é conceitualmente simples e corresponde exatamente à definição da distribuição Binomial Negativa. No entanto, quando pé pequeno e ré grande, o número esperado de ensaios E[X] = r/ppode ser elevado, tornando o algoritmo computacionalmente mais custoso. 2.6.2 Simulando via soma de geométricas Recorde que se X∼NegBin(r,p), então Xpode ser decomposto como X=X1+X2+···+Xr, onde X1, . . . , Xrsão variáveis independentes com Xi∼Geom(p),i=1, . . . , r, no suporte {1, 2, . . .}(número de ensaios até o primeiro sucesso). Assim, podemos simular uma Binomial Negativa somando rgeométricas independentes, cada uma gerada via o método da inversão: Xi=log(1−Ui) log(1−p)+1, Ui∼Uniforme(0, 1). O algoritmo de simulação segue: 1. Para i=1, . . . , r, gerar Xi∼Geom(p)via inversão; 2. Retornar X=∑r i=1Xi. 2.6. DISTRIBUIÇÃO BINOMIAL NEGATIVA 33 2.6.3 Simulando via inversão recursiva Outra forma de simular a Binomial Negativa é aplicar diretamente o método da inversão, aproveitando a relação de recorrência da sua função de probabilidade. Se X∼NegBin(r,p), então P(X=n) = n−1 r−1pr(1−p)n−r,n=r,r+1,r+2, . . . Essas probabilidades satisfazem a seguinte relação recursiva: P(X=n+1) P(X=n)=n n−r+1(1−p). Exercício 20. Prove a identidade recursiva acima. Portanto, conhecendo P(X=r) = pr, podemos calcular recursivamente as demais probabilidades. Isso leva ao seguinte algoritmo: 1. Gerar U∼Uniforme(0,1); 2. Inicializar n=r,pn=pr,F=pn; 3. Enquanto U>F, atualizar pn+1=pn·n n−r+1(1−p),n←n+1, F←F+pn+1; 4. Retornar X=n. O número esperado de passos Tno método recursivo não coincide diretamente com E[X], pois o algoritmo já inicia em n=r, que é o menor valor possível para a variável X∼NegBin(r,p). De fato, se o sorteio resultar em X=n, o número de passos dados pelo algoritmo é T= (n−r) + 1, pois começamos verificando o valor n=r(primeiro passo) e avançamos até alcançar n. Assim, em termos de valor esperado, E[T] = E[X−r] + 1=r p−r+1. Esse termo −raparece porque, embora E[X] = r/p, o procedimento de inversão não percorre todos os valores desde 0, mas já parte de r. Quando pé pequeno, E[T]pode ainda ser bastante grande, tornando o método recursivo lento. Nessas situações, a versão ingênua baseada na soma de geométricas pode ser mais eficiente na prática. 34 CAPÍTULO 2. VARIÁVEIS DISCRETAS E COMO SIMULÁ-LAS 2.6.4 Por que o nome “Binomial Negativa”? O nome Binomial Negativa tem origem na conexão com a expansão binomial para expoentes negativos. Para um inteiro n≥0 e p+q=1, o teorema binomial fornece 1= (p+q)n= n ∑ k=0n kpkqn−k. Essa identidade estende-se a expoentes reais (expansão binomial generalizada): (1−p)−α= ∞ ∑ k=0α+k−1 kpk,|p|<1. Tomando α=r∈ {1, 2, . . .}, (1−p)−r= ∞ ∑ k=0r+k−1 kpk. Os coeficientes (r+k−1 k)são precisamente os que aparecem na parametrização da Binomial Negativa em termos do número de falhas kantes do r-ésimo sucesso: P(Y=k) = r+k−1 kpr(1−p)k,k=0, 1,2, . . . Para ver a equivalência com a forma escrita em função do número total de ensaios n, detalhamos a reparametrização. Defina n=r+k(isto é, k=n−r). Então r+k−1 k=(r+k−1)! k!(r−1)!=(n−1)! (n−r)!(r−1)!=n−1 n−r. Pela simetria binomial, (a b)=(a a−b); aplicando com a=n−1 e b=n−robtemos n−1 n−r=n−1 (n−1)−(n−r)=n−1 r−1. Substituindo k=n−rem P(Y=k)e usando as igualdades acima, P(X=n) = P(Y=n−r) = r+ (n−r)−1 n−rpr(1−p)n−r =n−1 r−1pr(1−p)n−r,n=r,r+1, . . . Mostramos, assim, passo a passo, que as duas formas da PMF — em função de k(falhas) ou de n(ensaios) — são exatamente equivalentes; trata-se apenas de uma reparametrização. 2.7 Distribuição hipergeométrica A distribuição hipergeométrica modela experimentos de seleção sem reposição a partir de uma população finita contendo dois tipos de elementos. Por exemplo, suponha uma urna com N+M bolas, das quais Nsão claras e Msão escuras. Retiramos, de forma aleatória e sem reposição, uma amostra de tamanho n. Seja Xo número de bolas claras na amostra. 2.7. DISTRIBUIÇÃO HIPERGEOMÉTRICA 35 Nesse caso, cada subconjunto de tamanho né igualmente provável, e a probabilidade de observar exatamente kbolas claras é P(X=k) = (N k)( M n−k) (N+M n), max(0, n−M)≤k≤min(n,N). Dizemos então que X∼Hipergeom(N,M,n). As principais propriedades dessa distribuição são: E[X] = n·N N+M, Var(X) = n·N N+M·M N+M·N+M−n N+M−1. Exercício 21. Prove que a função de probabilidade acima é válida, isto é, que min(n,N) ∑ k=max(0, n−M) P(X=k) = 1. Exercício 22. Prove as fórmulas da média e variância acima. Dica: considere o sorteio sequencial das n bolas e defina Xicomo a variável indicadora do evento “a i-ésima bola é clara”. Para a variância, use a decomposição Var(X) = n ∑ i=1 Var(Xi) + 2∑ 1≤i<j≤n Cov(Xi,Xj). 2.7.1 Simulando a Hipergeométrica Uma maneira natural de simular X∼Hipergeom(N,K,n)é reproduzir o sorteio sem reposição. Basta imaginar uma população com Ksucessos e N−Kfracassos, e retirar nelementos dela. O valor de Xserá o número de sucessos observados. Isso leva ao seguinte algoritmo: 1. Construir a população formada por Kuns (sucessos) e N−Kzeros (fracassos); 2. Sortear nelementos dessa população sem reposição; 3. Definir Xcomo a soma dos elementos sorteados; 4. Retornar X. Para realizar o sorteio sem reposição, podemos usar um procedimento eficiente baseado no embaralhamento parcial de Fisher–Yates. A ideia é que não precisamos embaralhar toda a população, apenas selecionar nelementos distintos de forma aleatória. O algoritmo funciona assim: 1. Coloque os elementos da população em um vetor de tamanho N; 2. Para cada posição i=1, 2, . . . , n: (a) Sorteie um índice juniformemente entre ieN; (b) Troque os elementos das posições iej. 3. Os nprimeiros elementos do vetor agora constituem a amostra sem reposição. Esse método garante que cada subconjunto de tamanho ntem a mesma probabilidade de ser escolhido, e é mais eficiente do que embaralhar toda a população. 36 CAPÍTULO 2. VARIÁVEIS DISCRETAS E COMO SIMULÁ-LAS Figura 2.4: Distribuição Hipergeométrica com parâmetros N=50 (população total), K=20 (número de sucessos) e n=10 (tamanho da amostra). Tabela de referência Distribuição Técnica utilizada Dica/Obs Suporte finito Inversão simples Separar em intervalos Bernoulli Inversão simples Caso particular do suporte finito com m=2 Binomial Soma de Bernoullis ou inversão recursiva Relação de recorrência evita coeficientes binomiais Poisson Inversão recursiva Relação pi+1=λ i+1pievita fatoriais Geométrica Inversão direta (CDF) Retornar X=⌊log(1−U) log(1−p)⌋+1 Binomial negativa Soma de geométricas ou inversão recursiva Soma de rgeométricas independentes Hipergeométrica Sorteio sem reposição (Fisher–Yates) Selecionar nelementos distintos e contar sucessos Capítulo 3 Variáveis contínuas e como simulá-las Nosso objetivo agora é estudar algoritmos para simular variáveis aleatórias contínuas, isto é, variáveis cuja distribuição é descrita por uma função densidade de probabilidade. Como no caso discreto, o ponto de partida será sempre o mesmo: assumimos que temos acesso a uma variável U∼Uniforme(0,1), e construiremos a partir dela procedimentos para gerar amostras de outras distribuições. A principal diferença em relação ao caso discreto é que, para variáveis contínuas, muitas vezes não é possível escrever a função de distribuição acumulada (CDF) de forma explícita, ou mesmo obter sua inversa em forma fechada. Com isso, diversos métodos alternativos são necessários. Neste capítulo, organizamos os métodos de simulação em três grandes grupos: •Métodos por inversão: funcionam diretamente a partir da CDF da distribuição; •Métodos por rejeição ou aceitação: baseiam-se em gerar propostas e aceitar com certa probabilidade; •Métodos por transformação: aplicam funções determinísticas a variáveis já conhecidas. 3.1 Método da Inversão O método da inversão é uma das formas mais diretas de simular variáveis aleatórias contínuas. A ideia central é simples: se conhecemos a função de distribuição acumulada (CDF) Fde uma variável contínua X, e se essa função é estritamente crescente, então podemos inverter Fe definir X=F−1(U), com U∼Uniforme(0,1). Exemplo 7. Prove o fato acima. Esse procedimento garante que Xterá exatamente a distribuição desejada, pois a probabilidade de Xcair em qualquer intervalo será proporcional ao comprimento correspondente no domínio de U. Esse método é particularmente útil quando a inversa de Fpode ser escrita de forma explícita, como ocorre com as distribuições Exponencial, Uniforme e Pareto, por exemplo. 37 38 CAPÍTULO 3. VARIÁVEIS CONTÍNUAS E COMO SIMULÁ-LAS 3.1.1 Distribuição exponencial A distribuição exponencial modela o tempo de espera até a ocorrência de um evento em um processo de Poisson, isto é, um processo no qual eventos ocorrem de forma contínua e independente, a uma taxa constante. Seja λ>0 a taxa de ocorrência dos eventos. Dizemos que X∼Exponencial(λ)se fX(x) = λe−λx,x≥0. A função de distribuição acumulada (CDF) é dada por: FX(x) = P(X≤x) = 1−e−λx,x≥0. As principais características da distribuição são: E[X] = 1 λ, Var(X) = 1 λ2. Exercício 23. Verifique que fX(x)é uma densidade de probabilidade, isto é, R∞ 0fX(x)dx =1. Exercício 24. Prove as expressões da média e variância acima. A distribuição exponencial é um exemplo clássico onde o método da inversão pode ser aplicado diretamente. Sabemos que a CDF é F(x) = 1−e−λx, e queremos encontrar sua inversa, portanto, basta resolver a equação F(x) = U, com U∼ Uniforme(0,1). Isso nos leva a: 1−e−λx=U⇒x=−1 λlog(1−U). Como 1 −U∼Uniforme(0,1), podemos reescrever de forma equivalente: X=−1 λlog(U), com U∼Uniforme(0,1). Assim, o algoritmo para simular uma variável X∼Exponencial(λ)via inversão é: 1. Gerar U∼Uniforme(0,1); 2. Calcular X=−1 λlog(U); 3. Retornar X. Relação com a distribuição geométrica A distribuição exponencial pode ser vista como o análogo contínuo da distribuição geométrica. Na distribuição geométrica, X∼Geom(p), interpretamos Xcomo o número de tentativas independentes até a ocorrência do primeiro sucesso, em uma sequência de ensaios de Bernoulli com probabilidade pde sucesso. 3.1. MÉTODO DA INVERSÃO 39 A distribuição exponencial, por sua vez, modela o tempo contínuo até a ocorrência de um evento, sob uma taxa constante λ>0. Embora uma seja discreta e a outra contínua, existe uma relação direta entre essas duas distribuições, que pode ser formalizada por um limite. Seja Xn∼Geom(pn), com pn=λ/n, e defina a variável reescalada Tn=Xn n. A variável Tnrepresenta o tempo até o primeiro sucesso quando fazemos ntentativas por unidade de tempo, cada uma com probabilidade de sucesso pn=λ/n. À medida que n→∞, as tentativas se tornam mais frequentes e individualmente menos prováveis, mas o número esperado de sucessos por unidade de tempo permanece constante: n·pn=λ. Vamos mostrar que Tnconverge em distribuição para uma variável exponencial de parâmetro λ. De fato, temos: P(Tn>t) = PXn n>t=P(Xn>⌊nt⌋). Como Xné geométrica com parâmetro pn=λ/n, segue que: P(Xn>k) = (1−pn)k, logo P(Tn>t) = 1−λ n⌊nt⌋. Quando n→∞, vale que ⌊nt⌋ ∼ nt, e obtemos: 1−λ nnt −→ e−λt. Portanto, P(Tn≤t)→1−e−λt, que é a função de distribuição acumulada da exponencial Exp(λ). Isso conclui a demonstração da convergência. 40 CAPÍTULO 3. VARIÁVEIS CONTÍNUAS E COMO SIMULÁ-LAS Essa convergência tem uma interpretação intuitiva. Inicialmente, a variável Xnconta o número de tentativas até o sucesso. Se cada tentativa leva um tempo fixo de 1/nsegundos, então o tempo total até o sucesso é Tn=Xn/n. A divisão por nserve justamente para transformar o número de tentativas em tempo contínuo. Por exemplo, se cada tentativa leva 0.01 segundo e o sucesso ocorre na 17ª tentativa, então o tempo até o sucesso foi 17 ×0.01 =0.17 segundos. À medida que ncresce, as tentativas são feitas cada vez mais rapidamente (a cada 1/nunidades de tempo), e a chance de sucesso em cada uma cai proporcionalmente (pn=λ/n). O resultado final é que o tempo total até o sucesso — Tn— se aproxima de uma variável contínua exponencial com taxa λ. Essa relação também pode ser observada diretamente nas fórmulas de inversão utilizadas para simulação. Seja U∼Uniforme(0, 1). A inversão da CDF da exponencial dá: T=−1 λln(U). Já no caso da geométrica Xn∼Geom(pn), a fórmula de inversão baseada na CDF discreta é: Xn=ln(U) ln(1−pn), com pn=λ n. Dividindo por n, temos: Tn=Xn n≈1 n·ln(U) ln(1−λ/n). Sabemos que para ngrande, ln(1−λ/n)≈ −λ n, então: Tn≈ −1 λln(U), o que mostra que, no limite, a fórmula de simulação da geométrica reescalada tende para a fórmula da exponencial. Relação com a Poisson A distribuição exponencial pode ser entendida como o análogo contínuo da distribuição geométrica, e sua relação com a distribuição de Poisson surge naturalmente ao considerarmos divisões finas de um intervalo fixo em pequenos subintervalos com experimentos de Bernoulli raros. Considere o intervalo de tempo [0, 1]dividido em nsubintervalos de comprimento 1/n. Em cada subintervalo, ocorre um evento (ou sucesso) com probabilidade pn=λ/n, de forma independente. Este é exatamente o modelo da variável binomial Xn∼Binomial(n,λ/n), que conta o número total de eventos no intervalo. Sabemos que, quando n→∞, Xnd −→ Poisson(λ). 3.2. MÉTODO DA REJEIÇÃO-ACEITAÇÃO 41 Por outro lado, podemos perguntar: quanto tempo leva até o primeiro evento acontecer? A resposta a essa pergunta leva à distribuição exponencial. Seja Xn∼Geom(pn)com pn=λ/n, modelando o número de subintervalos até o primeiro sucesso. O tempo contínuo correspondente é então Tn=Xn n. Como visto anteriormente, temos Tnd −→ Exponencial(λ). Portanto, podemos pensar nessas distribuições da seguinte forma: • A distribuição Poisson modela o número total de eventos no intervalo. • A distribuição geométrica modela a posição discreta do primeiro sucesso. • A distribuição exponencial modela o tempo contínuo até o primeiro evento. 3.2 Método da rejeição-aceitação Embora o método da inversão funcione muito bem para distribuições cuja função de distribuição acumulada (CDF) possa ser invertida de forma analítica ou computacionalmente eficiente, ele se torna inviável em casos como o da distribuição normal padrão. A função de distribuição acumulada da normal, denotada por Φ(x), não possui inversa em forma fechada, o que impede a aplicação direta da fórmula X=Φ−1(U). Embora existam aproximações numéricas para Φ−1, elas podem ser computacionalmente custosas ou introduzir erros de arredondamento. Nesses casos, recorre-se a métodos alternativos que não exigem a inversão da CDF — como o método da rejeição-aceitação. O método de aceitação-rejeição é uma técnica geral para gerar variáveis aleatórias com uma dada densidade f(x), partindo de uma densidade auxiliar g(x)mais simples, da qual é fácil simular. A ideia central é gerar candidatos a partir de ge aceitá-los com uma certa probabilidade que depende da razão f(x)/g(x). 48 CAPÍTULO 3. VARIÁVEIS CONTÍNUAS E COMO SIMULÁ-LAS • Para α=1, a densidade é simplesmente a Exponencial decrescente. • Para α>1, a densidade é unimodal, com máximo em (α−1)θ. Já o parâmetro de escala θatua como fator multiplicativo, alongando ou comprimindo a distribuição. A média e a variância crescem proporcionalmente a ele, conforme: E[X] = αθ, Var(X) = αθ2. 3.3.1 Simulando quando αé inteiro Quando o parâmetro de forma é inteiro, α=k∈N, a distribuição Gamma(k,θ)(escala θ> 0) recebe o nome de Erlang. Ela pode ser obtida como a soma de kvariáveis exponenciais independentes. Na parametrização por taxa λ=1/θ: X∼Gamma(k,λ)⇐⇒ Xd = k ∑ i=1 Ei,Eii.i.d. ∼Exp(λ). O algoritmo de simulação é o seguinte: 1. Fixe k∈Neλ>0 (ou θ=1/λ). 2. Gere Ei∼Exp(λ)de forma independente, para i=1, . . . , k. 3. Calcule X=∑k i=1Ei. O resultado segue X∼Gamma(k,λ). 3.3.2 Simulando quando α>1via aceitação–rejeição com Exponencial Considere X∼Gamma(α,λ)com α>1 (parametrização por taxa λ). A densidade alvo é f(x) = λα Γ(α)xα−1e−λx,x>0. 3.3. DISTRIBUIÇÃO GAMMA 49 Usaremos como proposta Y∼Exp(µ), com densidade g(x) = µe−µx,x>0. Para que o método de aceitação–rejeição seja válido, precisamos de uma constante ctal que f(x)≤c g(x)para todo x>0. O quociente f(x) g(x)=λα Γ(α)µxα−1e−(λ−µ)x mostra que é necessário ter µ<λ, pois caso contrário o termo exponencial não decai e o quociente não tem máximo finito. Quando µ<λ, o máximo ocorre em x∗=α−1 λ−µ, com valor c(µ) = λα Γ(α)µα−1 λ−µα−1 e−(α−1). A probabilidade de aceitação por tentativa é 1/c(µ)e, portanto, o número médio de tentativas até aceitar uma amostra é c(µ). Interpretando em tempo contínuo como um processo de Poisson de taxa 1, o thinning com probabilidade 1/c(µ)gera um processo aceito com taxa 1/c(µ), de modo que o tempo médio entre aceitações é c(µ). O algoritmo é o seguinte: 1. Escolha µ∈(0, λ), por exemplo µ=µ∗=λ/αque minimiza c(µ). 2. Gere Y∼Exp(µ)eU∼Unif(0,1). 3. Aceite X=Yse U≤f(Y)/(c(µ)g(Y)); caso contrário, volte ao passo 2. 50 CAPÍTULO 3. VARIÁVEIS CONTÍNUAS E COMO SIMULÁ-LAS O valor aceito Xtem distribuição Gamma(α,λ). A tabela a seguir mostra a constante c∗e a taxa de aceitação 1/c∗para alguns valores de αna escolha ótima µ=λ/α(os valores independem de λ): α µ∗=λ/αc∗=αα Γ(α)e−(α−1)aceitação (1/c∗) 1.5 2 3λ1.2573 0.7953 21 2λ1.4715 0.6796 31 3λ1.8270 0.5473 51 5λ2.3848 0.4193 81 8λ3.0355 0.3294 12 1 12λ3.7306 0.2681 3.3.3 Simulando quando α<1 Para 0 <α<1, uma forma simples de simular Γ(α,θ)é usar a identidade se G∼Γ(α+1, θ)eU∼Unif(0, 1)(indep.), X=G U1/α∼Γ(α,θ). Para entender essa relação, considere as variáveis independentes (U,G)com U∼Unif(0,1) eG∼Γ(α+1, θ), 0 <α<1. Defina a transformação (x,g) = T(u,g) = g u1/α,g, cuja inversa é (u,g) = T−1(x,g) = (x/g)α,g. O suporte transformado é x>0 e g>x(pois u∈(0,1)implica x/g∈(0,1)). A densidade conjunta de (U,G)é fU,G(u,g) = fU(u)fG(g) = 1(0,1)(u)gαe−g/θ Γ(α+1)θα+1,g>0. 3.4. DISTRIBUIÇÃO BETA 51 Pela fórmula de mudança de variável, fX,G(x,g) = fU,G(x/g)α,gdet ∂(u,g) ∂(x,g). Calculamos o jacobiano usando a inversa u= (x/g)α: ∂u ∂x=αxα−1g−α,∂u ∂g=−αxαg−(α+1),∂g ∂x=0, ∂g ∂g=1, logo det ∂(u,g) ∂(x,g)=∂u ∂x·1−0=αxα−1g−α. Portanto, fX,G(x,g) = 1(x,∞)(g)gαe−g/θ Γ(α+1)θα+1αxα−1g−α=1(x,∞)(g)αxα−1 Γ(α+1)θα+1e−g/θ. Integrando em gpara obter a marginal de X: fX(x) = Z∞ xfX,G(x,g)dg =αxα−1 Γ(α+1)θα+1Z∞ xe−g/θdg =αxα−1 Γ(α+1)θα+1θe−x/θ. Usando Γ(α+1) = αΓ(α), obtemos fX(x) = xα−1e−x/θ Γ(α)θα,x>0, que é exatamente a densidade Γ(α,θ). Logo, X=G U1/α∼Γ(α,θ). O algoritmo de simulação é: 1. Dado 0 <α<1 e θ>0, defina α′=α+1. 2. Gere G∼Γ(α′,θ)(por exemplo, via Marsaglia–Tsang, pois α′>1). 3. Gere U∼Unif(0,1), independente de G. 4. Retorne X=G U1/α. Então X∼Γ(α,θ). 3.4 Distribuição Beta A distribuição Beta é uma das mais importantes distribuições contínuas em estatística, definida no intervalo unitário [0, 1]e parametrizada por dois parâmetros de forma α>0 e β>0. Sua densidade é dada por f(x) = 1 B(α,β)xα−1(1−x)β−1, 0 <x<1, onde B(α,β) = Z1 0uα−1(1−u)β−1du =Γ(α)Γ(β) Γ(α+β) é a função Beta de Euler. 52 CAPÍTULO 3. VARIÁVEIS CONTÍNUAS E COMO SIMULÁ-LAS A interpretação intuitiva da distribuição Beta é como um modelo de incerteza sobre probabilidades. Se pensamos em xcomo a probabilidade de sucesso em uma sequência de ensaios de Bernoulli, a Beta aparece naturalmente como distribuição a posteriori em modelos Bayesianos conjugados: começando com uma priori Beta(α,β), após observar ssucessos e ffracassos, a posteriori é Beta(α+s,β+f). A forma da densidade é bastante flexível: • Para α,β<1, a densidade concentra-se nos extremos 0 e 1. • Para α=β=1, temos a uniforme no intervalo (0,1). • Para α>1 e β>1, a densidade é unimodal, com máximo em (α−1)/(α+β−2). Os momentos principais são: E[X] = α α+β, Var(X) = αβ (α+β)2(α+β+1). Essas fórmulas mostram como αeβpodem ser interpretados como “pseudocontagens” de sucessos e fracassos, de forma que α+βcontrola a concentração da distribuição em torno da média. 3.4.1 Simulando a Beta via aceitação–rejeição com proposta uniforme A distribuição Beta(α,β)tem suporte em (0,1), de modo que uma escolha natural de proposta éY∼Unif(0, 1). A densidade da uniforme é g(y) = 1 para 0 <y<1, e precisamos de uma constante ctal que f(y)≤c g(y) = c, 0 <y<1. O algoritmo de aceitação–rejeição é: 1. Gere Y∼Unif(0,1)eU∼Unif(0,1)independentes. 3.5. TRANSFORMAÇÕES DE VARIÁVEIS ALEATÓRIAS 53 2. Aceite X=Yse U≤f(Y)/c, caso contrário repita o passo 1. O valor aceito Xterá distribuição Beta(α,β). A escolha ótima de cpode ser caracterizado em três casos: • Se α>1 e β>1, a densidade é unimodal, com modo em y∗=α−1 α+β−2, e portanto c=f(y∗) = 1 B(α,β)(y∗)α−1(1−y∗)β−1. • Se α=β=1, a distribuição é uniforme, logo c=1. • Se α<1 ou β<1, a densidade diverge em uma das extremidades (0 ou 1), e assim sup f(y) = ∞. Nesse caso não existe constante finita ce o método de aceitação–rejeição com uniforme como proposta não pode ser aplicado. 3.4.2 Simulando a Beta via Gammas independentes O método mais utilizado e geral para simular variáveis Beta(α,β)explora a relação entre as distribuições Beta e Gama. Seja G1∼Γ(α,1),G2∼Γ(β,1), independentes. Então vale a identidade X=G1 G1+G2∼Beta(α,β). A prova segue do fato de que o vetor normalizado (G1,G2)/(G1+G2)tem distribuição Dirichlet(α,β), e portanto sua primeira coordenada é Beta(α,β). Outra forma é calcular a densidade conjunta de (X,T), com X=G1 G1+G2eT=G1+G2, e verificar que a marginal de Xcoincide com a densidade da Beta. O algoritmo é simples e eficiente: 1. Gere G1∼Γ(α,1)eG2∼Γ(β,1)de forma independente. 2. Retorne X=G1 G1+G2. Esse procedimento funciona para qualquer α,β>0, inclusive quando são menores que 1, ao contrário do método de aceitação–rejeição com proposta uniforme. 3.5 Transformações de Variáveis Aleatórias Neste capítulo estudaremos transformações de variáveis e vetores aleatórios. A ideia central é a seguinte: dado um modelo probabilístico inicial, frequentemente precisamos aplicar funções a variáveis ou vetores aleatórios para obter novas quantidades de interesse. O objetivo, então, é caracterizar a distribuição resultante após a transformação, seja ela univariada ou multivariada. 54 CAPÍTULO 3. VARIÁVEIS CONTÍNUAS E COMO SIMULÁ-LAS 3.5.1 Geração de Normais via Método de Box–Muller Seja U∼Unif(0,2π)eT∼Expo(1)independentes, com densidade conjunta fU,T(u,t) = 1 2πe−t,u∈(0,2π),t>0. Definimos a transformação X=√2TcosU,Y=√2TsinU. O objetivo é determinar a densidade conjunta de (X,Y). Para isso usamos a fórmula de mudança de variáveis fX,Y(x,y) = fU,T(u,t)det ∂(u,t) ∂(x,y), onde (u,t)é obtido a partir de (x,y). Primeiro observamos que x2+y2=2t(cos2u+sin2u) = 2t, de modo que t=1 2(x2+y2),u=arctany x(ajustado para o quadrante correto). Assim, a transformação é invertível. O próximo passo é calcular o jacobiano da transformação direta (u,t)7→ (x,y). Temos ∂x ∂u=−√2tsin u,∂x ∂t=1 √2tcos u, ∂y ∂u=√2tcos u,∂y ∂t=1 √2tsin u. Logo, J= −√2tsin u1 √2tcos u √2tcos u1 √2tsin u . O determinante é det(J) = −√2tsin u 1 √2tsin u−1 √2tcos u√2tcos u. 3.5. TRANSFORMAÇÕES DE VARIÁVEIS ALEATÓRIAS 55 Simplificando, det(J) = −sin2u−cos2u=−1. Portanto, |det(J)|=1. Aplicando a fórmula de mudança de variáveis, fX,Y(x,y) = fU,T(u,t)det ∂(u,t) ∂(x,y)=1 2πe−t·1. Substituindo t=1 2(x2+y2), fX,Y(x,y) = 1 2πexp−1 2(x2+y2). Finalmente, notamos que fX,Y(x,y) = 1 √2πe−x2/2 1 √2πe−y2/2, o que mostra que XeYsão independentes e ambos têm distribuição Normal padrão. Com essa dedução, concluímos que (X,Y)definidos acima são variáveis independentes com distribuição Normal padrão. Assim, o método de Box–Muller pode ser usado diretamente para gerar Normais a partir de variáveis Uniformes e Exponenciais. Na prática, o algoritmo segue os seguintes passos: 1. Gere U∼Unif(0,2π). 56 CAPÍTULO 3. VARIÁVEIS CONTÍNUAS E COMO SIMULÁ-LAS 2. Gere T∼Expo(1). 3. Calcule X=√2TcosUeY=√2TsinU. 4. Então XeYsão independentes e possuem distribuição N(0,1). Intuitivamente, o que estamos fazendo é gerar um par de variáveis Normais independentes (X,Y)e representá-las em coordenadas polares: o ângulo é sorteado uniformemente e o raio vem de uma distribuição que garante a forma circular da densidade Normal. 3.5.2 Geração da normal bivariada Nosso objetivo agora é mostrar como simular uma Normal bivariada (Z,W)com marginais N(0,1)e correlação ρ, onde −1<ρ<1. A ideia é construir (Z,W)a partir de variáveis independentes mais simples. Sejam X,Y∼ N(0,1)independentes. Definimos Z=X,W=ρX+τY,τ=q1−ρ2. Então ZeWtêm marginais N(0,1)e Corr(Z,W) = ρ. Este é um procedimento construtivo que permite gerar diretamente a Normal bivariada a partir de duas Normais independentes. Para verificar a validade da construção, vamos obter a densidade conjunta de (Z,W). A densidade de (X,Y)é fX,Y(x,y) = 1 2πexp−1 2(x2+y2). Como a transformação é z=x,w=ρx+τy, a inversa é x=z,y=w−ρz τ. Aplicando a fórmula de mudança de variáveis, fZ,W(z,w) = fX,Y(x,y)det ∂(x,y) ∂(z,w), com (x,y)dados pela inversa acima. O jacobiano da transformação inversa é ∂(x,y) ∂(z,w)=  ∂x ∂z∂x ∂w ∂y ∂z∂y ∂w = 1 0 −ρ τ1 τ , portanto det ∂(x,y) ∂(z,w)=1 τ. Substituindo em fZ,W, fZ,W(z,w) = 1 2πτ exp−1 2z2+w−ρz τ2. 3.5. TRANSFORMAÇÕES DE VARIÁVEIS ALEATÓRIAS 57 Fazendo as contas e lembrando que ρ2+τ2=1, obtemos fZ,W(z,w) = 1 2πτ exp−1 2τ2(z2−2ρzw +w2). Essa é exatamente a forma conhecida da densidade Normal bivariada com matriz de covariância Σ= 1ρ ρ1!,|Σ|=1−ρ2=τ2. De fato, podemos escrever fZ,W(z,w) = 1 2π|Σ|1/2 exp−1 2(z,w)Σ−1(z,w)⊤. Do ponto de vista de simulação, esse resultado mostra que basta gerar X,Y∼ N(0, 1)independentes (e.g. via Box–Muller) e aplicar a transformação acima. O algoritmo é: 1. Gere X∼ N(0, 1). 2. Gere Y∼ N(0, 1)independentemente. 3. Calcule Z=X,W=ρX+q1−ρ2Y. 4. O par (Z,W)tem distribuição Normal bivariada com matriz de covariância Σ. Generalizando para normal multivariada No caso bivariado, construímos Z W!= 1 0 ρp1−ρ2  X Y!,X,Y∼ N(0,1)independentes. Chamando a matriz de transformação de A, temos Z W!=A X Y!,A= 1 0 ρp1−ρ2 . Como X,Ysão independentes com variância 1, temos Cov Z W!!=ACov X Y!!A⊤=AI2A⊤=AA⊤. Fazendo as contas, AA⊤= 1 0 ρp1−ρ2  1ρ 0p1−ρ2 = 1ρ ρ1 , que é exatamente a matriz de covariância da Normal bivariada desejada. 64 CAPÍTULO 4. SIMULAÇÃO VIA MONTE CARLO Geramos X1, . . . , Xn∼Uniforme(0,1)e usamos a identidade I=Ehe−X2i. Assim, o estimador de Monte Carlo é ˆ In=1 n n ∑ i=1 e−X2 i. O valor aproximado da integral é I ≈0.7468. Exemplo 10. Considere X ∼ N(0,1)e o evento A ={X>1}. Queremos estimar a probabilidade p=P(X>1). Geramos X1, . . . , Xn∼ N(0, 1)e usamos o estimador de Monte Carlo ˆ pn=1 n n ∑ i=1 1 {Xi>1}. Pela Lei dos Grandes Números, ˆ pn→p quando n →∞. O valor verdadeiro é p=1−Φ(1)≈0.1587, onde Φdenota a CDF da normal padrão. Exemplo 11. Podemos estimar o valor de πpor simulação de Monte Carlo usando uma interpretação geométrica. Considere o quadrado [0, 1]×[0, 1]e o quarto de círculo de raio 1centrado na origem, definido por x2+y2≤1. A área do quarto de círculo é π/4. Assim, se gerarmos pontos (Xi,Yi)uniformemente distribuídos no quadrado, a fração que cai dentro do círculo aproxima a razão entre as áreas, isto é, número de pontos no círculo número total de pontos ≈π 4. Logo, o estimador de Monte Carlo é ˆ πn=4×1 n n ∑ i=1 1 {X2 i+Y2 i≤1}. Pela Lei dos Grandes Números, ˆ πn→πquando n →∞. Exemplo 12. Considere X ∼Uniforme(0,1). Queremos estimar simultaneamente E[X]eVar [X]por simulação. Geramos X1, . . . , Xnindependentes e usamos ˆ µn=1 n n ∑ i=1 Xi,ˆ σ2 n=1 n−1 n ∑ i=1 (Xi−ˆ µn)2. 4.1. ESTIMANDO MÉDIAS 65 Exemplo 13. Considere a integral em duas dimensões I=Z1 0Z1 0e−(x2+y2)dx dy. Geramos (Xi,Yi)∼Uniforme([0,1]2)e usamos o estimador ˆ In=1 n n ∑ i=1 e−(X2 i+Y2 i). Exemplo 14. Considere X ∼Exponencial(1)e queremos estimar Ee−X. Geramos X1, . . . , Xnindependentes e calculamos ˆ µn=1 n n ∑ i=1 e−Xi. O valor verdadeiro é Ee−X=1 2. Exemplo 15. Considere a integral I=Z∞ 0 e−x 1+xdx. Essa integral não tem forma fechada simples, mas pode ser expressa como uma esperança sob uma distribuição conveniente. Observe que f(x) = e−x1{x>0}é a densidade de uma distribuição Exponencial(1). Assim, podemos escrever I=Z∞ 0 e−x 1+xdx =Z∞ 0 1 1+xf(x)dx =E1 1+X,X∼Exponencial(1). Logo, podemos estimar I por Monte Carlo gerando X1, . . . , Xn∼Exponencial(1)e calculando ˆ In=1 n n ∑ i=1 1 1+Xi. O valor verdadeiro da integral é aproximadamente I ≈0.5963. Exemplo 16. Considere a integral I=Z∞ 0 sin(x) xdx. Não há uma densidade de probabilidade aparecendo explicitamente, mas podemos introduzir uma para reescrever a integral como uma esperança. Escolha, por conveniência, a densidade exponencial f(x) = e−x1{x>0}. Então, I=Z∞ 0 sin(x) xdx =Z∞ 0 sin(x) xe−xf(x)dx =Esin(X) Xe−X,X∼Exponencial(1). Assim, a integral pode ser estimada por Monte Carlo sem precisar integrar diretamente uma função oscilatória: ˆ In=1 n n ∑ i=1 sin(Xi) Xie−Xi,Xi∼Exponencial(1). O valor exato da integral é I =π 2≈1.5708. 66 CAPÍTULO 4. SIMULAÇÃO VIA MONTE CARLO 4.2 Intervalos de Confiança Uma estimativa obtida por simulação de Monte Carlo é aleatória. Mesmo quando o estimador é não tendencioso, seu valor varia a cada execução devido à variabilidade amostral. Por isso, é importante quantificar essa incerteza por meio de um intervalo de confiança. Seja ˆ µn=1 n∑n i=1h(Xi)um estimador de Monte Carlo para µ=E[h(X)]. Pelo Teorema Central do Limite, √nˆ µn−µ σ d −→ N(0,1), onde σ2=Var [h(X)]. Assim, para ngrande, P|ˆ µn−µ| ≤ z1−α/2 σ √n≈1−α, onde z1−α/2 é o quantil da normal padrão. Na prática, a variância σ2é desconhecida. Usamos a estimativa amostral s2 n=1 n−1 n ∑ i=1 (h(Xi)−ˆ µn)2. Substituindo σpor sn, obtemos o intervalo de confiança assintótico ˆ µn±z1−α/2 sn √n. O intervalo representa a faixa de valores plausíveis para µ, dada a variabilidade da amostra. Para um nível de confiança de 95%, usamos z0.975 ≈1.96, obtendo ˆ µn±1.96 sn √n. Exemplo 17. Considere novamente a estimativa I=Ehe−X2i,X∼Uniforme(0,1). Queremos construir um intervalo de confiança para I com base em n =105simulações. (1) O estimador de Monte Carlo é ˆ In=1 n n ∑ i=1 e−X2 i. A média amostral obtida foi ˆ In=0.7472. (2) O desvio padrão amostral das observações e−X2 ié sn=0.289. (3) O erro padrão do estimador é EP =sn √n=0.289 √100000 =0.000914. 4.2. INTERVALOS DE CONFIANÇA 67 (4) Pela regra empírica 68–95–99.7, sabemos que: –cerca de 68% das observações estão dentro de 1desvio padrão da média; –cerca de 95% estão dentro de 2desvios padrão; –e cerca de 99.7% estão dentro de 3desvios padrão. Assim, podemos construir um intervalo aproximado de 95% de confiança usando dois desvios padrão em vez de 1.96. (5) O termo de margem é então 2×EP =2×0.000914 =0.00183. (6) O intervalo de confiança é [0.7472 −0.00183, 0.7472 +0.00183 ]=[0.7454, 0.7490]. Em outras palavras, esperamos que cerca de 95% das repetições do experimento de Monte Carlo produzam valores de ˆ Indentro de dois erros padrão da média verdadeira. O valor teórico I =0.7468 está de fato dentro desse intervalo. Exercício 25. Ache intervalos de confianças para todos os exemplos anteriores. 68 CAPÍTULO 4. SIMULAÇÃO VIA MONTE CARLO Capítulo 5 Redução de variância Em um estudo de simulação, é comum que se deseje estimar um parâmetro θassociado a um modelo estocástico. Para isso, o modelo é executado a fim de gerar uma variável de saída X, cuja esperança é θ=E[X]. Realizam-se então nrepetições independentes da simulação, sendo que a i-ésima repetição fornece o valor Xi. A partir dessas observações, a estimativa natural de θé a média amostral X=1 n n ∑ i=1 Xi. Note que Xé um estimador não viesado de θ, de modo que E[X] = θ. Assim, o erro quadrático médio do estimador coincide com sua variância: MSE(X) = E(X−θ)2=Var(X) = Var(X) n. Portanto, se for possível construir um outro estimador não viesado de θcom variância menor do que a de X, obteremos uma estimativa mais eficiente. Este é o ponto de partida para as técnicas de redução de variância que discutiremos a seguir. 5.1 Uso de variáveis antitéticas Considere o problema de estimar θ=E[X]por simulação. Se gerarmos duas observações X1e X2, identicamente distribuídas com esperança θ, uma estimativa natural é a média ˆ θ=X1+X2 2. A variância desse estimador pode ser escrita como Var(ˆ θ) = 1 4Var(X1+X2) = 1 4Var(X1) + Var(X2) + 2Cov(X1,X2). Como X1eX2têm a mesma distribuição, Var(X1) = Var(X2), segue que Var(ˆ θ) = 1 2Var(X1) + 1 2Cov(X1,X2). Portanto: 69 70 CAPÍTULO 5. REDUÇÃO DE VARIÂNCIA • se X1eX2forem independentes, Cov(X1,X2) = 0 e Var(ˆ θ) = 1 2Var(X1); • se conseguirmos construir X1eX2de modo que a covariância seja negativa, então a variância de ˆ θserá ainda menor. A questão, então, é: como gerar dois valores X1eX2com a mesma distribuição, mas negativamente correlacionados? Suponha que X1seja função de mnúmeros aleatórios independentes, isto é, X1=h(U1, . . . , Um), onde U1, . . . , Umsão independentes e uniformemente distribuídos em (0,1). Observe que, se U∼U(0,1), então também 1 −U∼U(0, 1). Assim, se definirmos X2=h(1−U1, . . . , 1 −Um), teremos que X2possui a mesma distribuição que X1. Além disso, como 1 −Ué negativamente correlacionado com U, é razoável esperar que X2 seja negativamente correlacionado com X1. Para tornar a ideia mais clara, considere o caso em que X1depende apenas de uma variável uniforme. Seja U∼U(0, 1)e uma função monótona crescente h:[0,1]→R. Definimos X1=h(U),X2=h(1−U). Note que X1eX2têm a mesma distribuição, pois Ue 1 −Usão identicamente distribuídos. Além disso, se Uassume um valor grande, então X1=h(U)também será grande, mas nesse caso 1 −Userá pequeno, de modo que X2=h(1−U)será pequeno. Assim, valores altos de X1 tendem a estar associados a valores baixos de X2, e vice-versa, o que implica correlação negativa. No caso particular em que h(u) = u, temos X1=U,X2=1−U. Claramente, E[X1] = E[X2] = 1 2, de modo que E[X1]E[X2] = 1 4. Por outro lado, E[X1X2] = Z1 0u(1−u)du =Z1 0(u−u2)du =1 2−1 3=1 6. Assim, Cov(X1,X2) = 1 6−1 4=−1 12 <0, mostrando explicitamente a correlação negativa entre X1eX2. Esse raciocínio se estende naturalmente para funções hmonótonas em várias variáveis. Sejam U1, . . . , Umvariáveis independentes uniformes em (0,1)e definamos X1=h(U1, . . . , Um),X2=h(1−U1, . . . , 1 −Um), 5.1. USO DE VARIÁVEIS ANTITÉTICAS 71 com hcrescente em cada coordenada. Nesse caso, X1é uma função crescente do vetor (U1, . . . ,Um), enquanto X2é decrescente. Considerando g(U1, . . . , Um) = −X2, vemos que gtambém é crescente em cada coordenada. Ora, quando duas funções de um mesmo conjunto de variáveis independentes são monótonas no mesmo sentido (ambas crescentes ou ambas decrescentes), seus valores tendem a variar em conjunto, de modo que a covariância é não-negativa. Aplicando esse raciocínio a X1eg, concluímos que Cov(X1,−X2)≥0, o que implica Cov(X1,X2)≤0. Portanto, para qualquer função hcrescente em cada coordenada, o par (X1,X2)é negativamente correlacionado, e o uso de variáveis antitéticas reduz (ou, no pior caso, não aumenta) a variância do estimador. Em resumo, o método das variáveis antitéticas consiste em explorar a correlação negativa entre pares de simulações para reduzir a variância do estimador. Em vez de gerar duas réplicas independentes X1eX2, construímos o par de forma que ambas tenham a mesma distribuição, mas sejam negativamente correlacionadas. A variável antitética é dado por X′=X1+X2 2, o qual satisfaz E[X′] = θ, mas possui variância menor ou igual à do estimador baseado em amostras independentes. Um algoritmo simples para aplicar o método pode ser descrito da seguinte forma: 1. Gere U1, . . . ,Um∼Uniforme(0, 1)independentes. 2. Calcule X1=h(U1, . . . , Um). 3. Calcule X2=h(1−U1, . . . , 1 −Um). 4. Defina a variável antitética como X′=X1+X2 2. Sempre que utilizamos o método da inversão para gerar variáveis aleatórias, podemos aplicar diretamente a técnica das variáveis antitéticas. De fato, se U∼U(0,1)gera a variável desejada via a transformação X=F−1(U), então 1 −Utambém é uniforme em (0, 1), e portanto X′= F−1(1−U)tem a mesma distribuição de X. A grande vantagem é que, em vez de gerar duas variáveis independentes U1eU2para obter duas amostras de X, basta gerar uma única variável uniforme U. Com ela, obtemos simultaneamente o par antitético (X,X′), o que não apenas economiza custo computacional como também pode reduzir a variância do resultado final. 72 CAPÍTULO 5. REDUÇÃO DE VARIÂNCIA Exemplo 18. Considere a geração de uma variável aleatória exponencial com parâmetro λ>0. Pelo método da inversão, se U ∼U(0,1), então X=−1 λlog(U) segue a distribuição Exp(λ). Para aplicar o método das variáveis antitéticas, em vez de gerar duas variáveis independentes U1,U2∼ U(0,1), usamos o par (U, 1 −U). Assim, obtemos X1=−1 λlog(U),X2=−1 λlog(1−U). Definimos, então, a variável final como a média Z=X1+X2 2. O algoritmo é: 1. Gere U ∼Uniforme(0,1). 2. Calcule X1=−1 λlog(U). 3. Calcule X2=−1 λlog(1−U). 4. Defina Z = (X1+X2)/2. No método independente, como cada variável exponencial pode assumir valores próximos de zero (quando U →1), a média também pode se aproximar de zero. Já no método antitético temos, supondo λ=1, Z=−1 2log U(1−U). 5.1. USO DE VARIÁVEIS ANTITÉTICAS 73 Como U(1−U)≤1/4, segue que Z≥1 2log4 =log2 ≈0.693. Ou seja, a variável construída por antitéticos nunca assume valores menores que log2. Esse resultado explica por que, ao comparar os histogramas, a média independente pode assumir valores próximos de zero, enquanto a antitética tem suporte a partir de log 2. Além disso, no experimento com n=105, o erro quadrático médio foi aproximadamente 0.505 no caso independente e apenas 0.174 no caso antitético, mostrando a expressiva redução de variância obtida pelo método. Exemplo 19. Considere a integral I=Z∞ 0log(1+x2)e−xdx. Observe que o termo e−xcorresponde à densidade de uma variável X ∼Exp(1). Assim, podemos reescrever a integral como I=Elog(1+X2),X∼Exp(1). Portanto, a solução via Monte Carlo é imediata: basta gerar amostras Xi∼Exp(1), calcular log(1+ X2 i)e tirar a média. O algoritmo segue os passos: 1. Gerar Ui∼U(0,1). 2. Transformar em Xi=−log(Ui). 3. Calcular log(1+X2 i)e tirar a média. Para reduzir a variância, podemos usar variáveis antitéticas. Nesse caso, ao invés de gerar apenas Ui, usamos também 1−Ui. Isso produz Xi=−log(Ui),X′ i=−log(1−Ui), e então o estimador final é ˆ Iant =1 n n ∑ i=1 1 2log(1+X2 i) + log(1+ (X′ i)2). Note que nem sempre variáveis antitéticas reduzem a variância: essa técnica é mais eficaz quando a função aplicada às amostras (aqui, log(1+x2)) é monotônica, pois nesse caso os pares (U,1 −U)tendem a gerar correlação negativa entre os valores simulados. Exemplo 20. Considere a integral J=Z∞ −∞ex1 √2πe−x2/2 dx. O integrando envolve a densidade da normal padrão N(0,1), logo podemos escrever J=E[eZ],Z∼N(0,1). O valor exato é conhecido: J=e1/2. Para estimar J via Monte Carlo, seguimos os passos: 80 CAPÍTULO 5. REDUÇÃO DE VARIÂNCIA Como V2∼U(−1,1)e é independente de V1, temos PV2 2≤1−v2=P−p1−v2≤V2≤p1−v2. Portanto, E[I|V1=v] = Z√1−v2 −√1−v2 1 2dx =p1−v2. Assim, obtemos o estimador condicionado Z=q1−V2 1, que satisfaz E[Z] = π/4, mas tem variância menor do que I. Finalmente, se U ∼U(0,1), temos V1=2U−1, e portanto Z=q1−(2U−1)2. Logo, podemos simular πa partir do estimador ˆ π=4 n n ∑ i=1q1−(2Ui−1)2, que é mais eficiente do que usar diretamente o indicador I. Exemplo 23. Considere a seguinte modelagem para a altura em uma população. Seja Y ∼Bernoulli(p) a variável que indica o sexo do indivíduo, em que Y =0representa mulher e Y =1representa homem. Condicionalmente a Y, a altura X tem distribuição normal X|Y=0∼N(µf,σ2 f),X|Y=1∼N(µm,σ2 m). Nosso objetivo é estimar a média da população, θ=E[X]=(1−p)µf+pµm. Se gerarmos indivíduos completos, isto é, sorteando Y e depois X |Y, o estimador de Monte Carlo é ˆ θsimples =1 n n ∑ i=1 Xi, que é não viesado para θ, mas apresenta variância Var(X) = (1−p)σ2 f+pσ2 m+p(1−p)(µm−µf)2. Em vez de usar diretamente X, podemos aplicar condicionamento. Nesse caso, registramos Z =E[X| Y], ou seja, µfse Y =0eµmse Y =1. O estimador correspondente é ˆ θcond =1 n n ∑ i=1 Zi, que também é não viesado para θ, mas com variância Var(Z) = p(1−p)(µm−µf)2, estritamente menor do que Var(X), pois elimina a variabilidade interna de cada grupo (σ2 feσ2 m). Assim, ao invés de considerar a altura ruidosa de cada indivíduo, utilizamos a média condicional do grupo, que é mais estável e resulta em um estimador mais eficiente. Capítulo 6 Amostragem por importância Considere uma variável aleatória Xcom densidade f(x). Nosso objetivo é calcular θ=E[h(X)] = Zh(x)f(x)dx. Em alguns casos, uma simulação direta de X∼fpode ser ineficiente: • pode ser difícil gerar amostras segundo f; • a variância de h(X)sob fpode ser grande; • ou ainda uma combinação desses fatores. Uma alternativa é escolher uma outra densidade g(x)tal que f(x) = 0 sempre que g(x) = 0. Nesse caso, podemos reescrever θ=Zh(x)f(x) g(x)g(x)dx =Egh(X)f(X) g(X), onde Egdenota esperança em relação à densidade g. Assim, se gerarmos X1, . . . , Xn∼g, um estimador natural é ˆ θ=1 n n ∑ i=1 h(Xi)f(Xi) g(Xi). Se a densidade instrumental gfor bem escolhida, a variância do peso h(X)f(X)/g(X)pode ser bem menor do que a variância de h(X)sob f, resultando em uma estimação mais eficiente. Note que f(x)eg(x)representam as probabilidades relativas de se observar xquando X∼f ou X∼g. Quando X∼g, em geral a razão f(x)/g(x)é menor que 1, mas como Egf(X) g(X)=1, ela ocasionalmente assume valores grandes, podendo gerar alta variância. A ideia central da amostragem por importância é escolher gde modo que nesses pontos onde f(x)/g(x)é grande, a função h(x)seja pequena (ou mesmo nula). Dessa forma, o produto W(X) = h(X)f(X) g(X) 81 82 CAPÍTULO 6. AMOSTRAGEM POR IMPORTÂNCIA permanece controlado, evitando explosões na variância. Esse raciocínio mostra por que a técnica é especialmente eficaz na estimação de probabilidades raras. Nesse caso, h(x)é uma função indicadora de um conjunto Apouco provável sob f. Se escolhermos gde modo que Aseja mais frequente, então: • para x∈A, temos h(x) = 1 e a razão f(x)/g(x)é moderada; • para x/∈A, temos h(x) = 0, logo não importa se f(x)/g(x)é grande. Assim, o estimador se torna muito mais estável e com variância reduzida, o que torna a amostragem por importância uma ferramenta poderosa para lidar com eventos raros. 6.1 Densidades Inclinadas (Tilted Densities) Uma questão central em amostragem por importância é a escolha da densidade instrumental g(x). Uma família bastante útil é a das densidades inclinadas, definidas a partir da função geradora de momentos. Seja X∼fuma variável aleatória com função geradora de momentos M(t) = Ef[etX] = Zetx f(x)dx. Definição 1. A densidade inclinada de f, associada ao parâmetro t ∈R, é definida por ft(x) = etx f(x) M(t). Intuitivamente, a densidade ftdá mais peso a valores grandes de Xquando t>0 e mais peso a valores pequenos quando t<0. Em muitos casos, ftpertence à mesma família paramétrica de f, mas com parâmetros modificados. Alguns exemplos: •Normal. Se X∼N(µ,σ2), então ftéN(µ+tσ2,σ2). Nesse caso, o tilt desloca a média, concentrando a massa de probabilidade à direita quando t>0 e à esquerda quando t<0. •Exponencial. Se X∼Exp(λ), então fté Exp(λ−t), válido para t<λ. Aqui, o tilt altera o comportamento da cauda: para t>0, a distribuição decai mais rápido (cauda mais leve), enquanto para t<0 a cauda se torna mais pesada. •Gama. Se X∼Gamma(α,β), então fté Gamma(α,β−t), válido para t<β. Assim como na exponencial (caso particular da gama), o tilt controla a espessura da cauda, deixando-a mais leve quando t>0 e mais pesada quando t<0. •Poisson. Se X∼Poisson(λ), então fté Poisson(λet). Nesse caso, o tilt modifica a média exponencialmente: para t>0, a distribuição se desloca para valores grandes, enquanto para t<0 se concentra em valores pequenos. 6.1. DENSIDADES INCLINADAS (TILTED DENSITIES) 83 •Binomial. Se X∼Binomial(n,p), então fté Binomial(n,pt)com pt=pet 1−p+pet. Aqui, o tilt altera diretamente a probabilidade de sucesso: quando t>0 temos pt>p, o que força mais sucessos, e quando t<0 temos pt<p, forçando mais fracassos. Exercício 26. Prove as afirmações acima. Exemplo 24 (Estimando probabilidades raras).Sejam X1, . . . , Xnvariáveis aleatórias independentes com densidades (funções de massa ou de probabilidade) fi, para i =1, . . . , n. Defina S= n ∑ i=1 Xi,µ= n ∑ i=1 E[Xi]. Nosso objetivo é estimar a probabilidade de que S seja maior do que um limiar a, onde a ≫µ, isto é, θ=P(S>a) = E1{S>a}. Quando a é muito maior que a média µ, esse evento é raro, e portanto uma simulação direta via Monte Carlo ingênuo é ineficiente, pois apenas uma fração ínfima das amostras contribui para o cálculo do estimador. Uma alternativa é utilizar a técnica de amostragem por importância com densidades inclinadas. Seja fi,t(x) = etx fi(x) Mi(t),Mi(t) = E[etXi], a densidade inclinada de Xi, onde t >0é um parâmetro comum a todas as variáveis. Ao simular cada Xi segundo fi,t, obtemos que θ=Et"1{S>a}exp(−tS) n ∏ i=1 Mi(t)#, onde Etdenota esperança sob as densidades inclinadas. A partir dessa representação, segue naturalmente um estimador de Monte Carlo: ˆ θ=1 N N ∑ j=1 1{S(j)>a}exp−tS(j)n ∏ i=1 Mi(t), onde S(j)=∑n i=1X(j) ie cada X(j) ié simulado de fi,t. A escolha de t é crucial: se for muito pequeno, a distribuição inclinada pouco difere da original e o evento {S>a}continua raro. Se for muito grande, os pesos podem se tornar instáveis, aumentando a variância. O critério usual é escolher t de forma que Et[S]≈a, ou seja, deslocar a média da soma sob a medida inclinada para próximo do limiar a. Dessa forma, amostras são concentradas justamente nas regiões que mais contribuem para o evento raro, aumentando a eficiência do método. O algoritmo pode ser resumido da seguinte forma: 84 CAPÍTULO 6. AMOSTRAGEM POR IMPORTÂNCIA 1. Escolha t >0de modo que Et[S]≈a. 2. Para j =1, . . . , N: (a) Gere X(j) 1, . . . , X(j) nindependentemente segundo as densidades inclinadas fi,t. (b) Calcule S(j)=∑n i=1X(j) i. (c) Associe o peso W(j)=1{S(j)>a}exp(−tS(j)) n ∏ i=1 Mi(t). 3. Estime θpor ˆ θ=1 N N ∑ j=1 W(j). No caso particular em que cada Xi∼N(0,1), temos que a soma S =∑n i=1Xié normal N(0, n). Quando n =1, por exemplo, estimar P(S>5)é um evento extremamente raro, pois o valor exato dessa probabilidade é θ=P(S>5)≈2.87 ×10−7. Um procedimento de Monte Carlo ingênuo, baseado apenas em amostrar de N(0,1), é ineficiente: em uma simulação com N =200,000 repetições, a estimativa obtida foi de aproximadamente 5.0 ×10−6, um valor que não coincide com o verdadeiro devido à raridade do evento. Aplicando a técnica de densidades inclinadas, obtemos que a tilted density é ft(x) = etx f(x) M(t),M(t) = exp1 2t2, o que implica que fté a densidade de uma normal N(t,1). Ao escolher t =5, a distribuição inclinada desloca a média exatamente para o limiar de interesse. O peso associado a cada amostra é dado por W=1{S>5}exp(−tS +1 2t2). Nesse caso, a estimativa via amostragem por importância com N =200,000 simulações foi ˆ θtilt ≈2.87 ×10−7, em perfeito acordo com o valor teórico. 6.2. DESIGUALDADE DE CHERNOFF 85 6.2 Desigualdade de Chernoff A inclinação exponencial (exponential tilting) introduzida anteriormente também pode ser entendida como uma mudança de medida aplicada diretamente às probabilidades. Essa é, essencialmente, a mesma ideia usada no exemplo normal, em que substituímos f=N(0,1)por sua versão inclinada fλ=N(λ,1), de modo que o evento raro {X>10}se torne típico sob a nova distribuição. Seja Xuma variável aleatória com densidade f, e defina a densidade inclinada: fλ(x) = eλxf(x) Z(λ),Z(λ) = Ef[eλX]. Podemos então expressar qualquer probabilidade como uma esperança sob essa nova medida: P(X≥a) = Zx≥af(x)dx =Zx≥a f(x) fλ(x)fλ(x)dx =Eλf(X) fλ(X)1{X≥a}. Usando a definição de fλ, a razão de verossimilhança é f(X) fλ(X)=e−λX+ψ(λ), e portanto P(X≥a) = Eλhe−λX+ψ(λ)1{X≥a}i=Z(λ)Eλhe−λX1{X≥a}i. Essa identidade expressa a probabilidade de um evento raro como uma esperança sob a medida inclinada fλ. Em princípio, essa igualdade poderia ser usada para estimação — poderíamos simular de fλe calcular a média dos pesos e−λX+ψ(λ)1{X≥a}, exatamente como em importance sampling. No entanto, se o objetivo não é estimar, mas limitar a probabilidade, podemos substituir o peso aleatório e−λXpor um limite superior determinístico que vale no evento de interesse. No evento {X≥a}, temos e−λX≤e−λa. Aplicando essa desigualdade dentro da esperança obtemos: P(X≥a)≤e−λaEλ[eψ(λ)1{X≥a}] = e−λa+ψ(λ)Pλ(X≥a). Como Pλ(X≥a)≤1, chegamos finalmente a P(X≥a)≤exp−λa+ψ(λ). Esse passo transforma a identidade exata do importance sampling em um limite superior determinístico — a Desigualdade de Chernoff. Mostra que a mesma inclinação exponencial usada para redução de variância em estimação Monte Carlo também fornece uma maneira analítica elegante de controlar probabilidades de eventos raros. O limite obtido acima depende do parâmetro λ. Diferentes valores de λcorrespondem a diferentes distribuições inclinadas fλe, portanto, a diferentes mudanças de medida. Para obter o limite mais apertado, minimizamos o expoente: P(X≥a)≤inf λ>0exp−λa+ψ(λ). O valor ótimo λ⋆satisfaz a condição de primeira ordem: ψ′(λ⋆) = a. 86 CAPÍTULO 6. AMOSTRAGEM POR IMPORTÂNCIA Para entender a condição para λ⋆ótimo, calculemos a derivada da função log-partição. A partir de ψ(λ) = log Zeλxf(x)dx, derivando em relação a λobtemos: ψ′(λ) = Rxeλxf(x)dx Reλxf(x)dx . Essa expressão pode ser reconhecida como a média de Xsob a densidade inclinada fλ(x)∝ eλxf(x): ψ′(λ) = Eλ[X]. Portanto, a derivada da função log-partição coincide com o valor esperado de Xsob a inclinação exponencial. No valor ótimo λ⋆, temos: Eλ[X] = ψ′(λ⋆) = a, o que significa que, sob a inclinação ótima, a média de Xé igual ao limiar a. Em termos probabilísticos, isso mostra que a distribuição fλ⋆torna o evento {X≥a}típico — seu valor médio já se encontra na fronteira da região rara que estamos tentando estudar. 6.3 Variância sob Inclinação Exponencial Já vimos que a medida inclinada fλtorna a região de interesse, como X>a, típica. A inclinação ótima é alcançada quando ψ′(λ) = a, pois ψ′(λ) = Eλ[X]é a média de Xsob a distribuição inclinada. Em outras palavras, o parâmetro λdesloca a distribuição de modo que sua esperança coincida com o ponto que queremos estimar. Para avaliar a qualidade dessa reponderação, uma quantidade natural a estudar é a variância do estimador sob a medida inclinada. Se a distribuição inclinada permanece altamente concentrada em torno de sua média, os pesos do importance sampling são estáveis e o estimador é eficiente. Por outro lado, se a lei inclinada é muito dispersa, os pesos flutuam fortemente e o estimador sofre com alta variância. Assim, a concentração da distribuição inclinada fornece uma medida direta da qualidade do estimador de importance sampling. Essa concentração é capturada pela segunda derivada da função log-partição. De fato, a partir de ψ(λ) = logEheλXi, temos ψ′(λ) = EXeλX E[eλX]=Eλ[X], que representa a média de Xsob a distribuição inclinada. Derivando novamente, ψ′′(λ) = EX2eλX E[eλX]− EXeλX E[eλX]!2 =EλX2−(Eλ[X])2=Varλ[X]. Portanto, a curvatura da log-partição quantifica quão concentrada é a medida inclinada em torno de sua média e, consequentemente, quão eficiente será o estimador de importance sampling. 6.3. VARIÂNCIA SOB INCLINAÇÃO EXPONENCIAL 87 Em alguns casos, essa variância pode ser uniformemente limitada para todos os valores de λ. Por exemplo, quando Xé uma variável aleatória limitada tal que X∈[a,b], a variância sob qualquer medida inclinada satisfaz Varλ[X]≤(b−a)2 4. De fato, a distância de Xao ponto médio do intervalo (a,b)é sempre menor que metade do comprimento do intervalo, isto é, X−a+b 2≤b−a 2. Seja m=a+b 2. Então, (X−m)2≤b−a 22 . Tomando expectativas sob a lei inclinada, obtemos Eλ(X−m)2≤b−a 22 . Além disso, para qualquer constante c, Varλ[X]=min c∈R Eλ(X−c)2≤Eλ(X−m)2, de modo que Varλ[X]≤b−a 22 . Isso mostra que, para qualquer variável aleatória limitada, a variância sob inclinação exponencial permanece uniformemente controlada. Em particular, a curvatura da função log-partição — que determina tanto a concentração da lei inclinada quanto a estabilidade do estimador de importance sampling — não pode crescer sem limite. Agora suponha que Xé centrada, isto é, E[X]=0. Então, por definição, ψ(0) = logEhe0·Xi=0, ψ′(0) = E[X]=0. Pela expansão de Taylor de segunda ordem de ψ, existe algum θ∈(0, λ)tal que ψ(λ) = ψ(0) + ψ′(0)λ+λ2 2ψ′′(θ) = λ2 2Varθ[X]≤λ2 2sup θ∈(0,λ) Varθ[X]. Em particular, conhecer o comportamento da variância sob inclinação permite controlar toda a forma da função log-partição. Se a variância inclinada permanece uniformemente limitada, a curvatura de ψtambém é limitada, e os momentos exponenciais de Xcrescem no máximo quadraticamente em λ. No caso especial de variáveis limitadas, combinando isso com o limite uniforme sobre a variância inclinada obtemos ψ(λ)≤λ2(b−a)2 8. Esse resultado mostra que, sempre que a variância sob inclinação exponencial é uniformemente limitada, a função log-partição cresce no máximo quadraticamente em λ. O crescimento quadrático da log-partição é precisamente a marca do comportamento sub-Gaussiano. Na próxima seção, formalizamos essa conexão e mostramos como ela permite controlar probabilidades de eventos raros mesmo quando a função log-partição exata é desconhecida. 88 CAPÍTULO 6. AMOSTRAGEM POR IMPORTÂNCIA 6.4 Variáveis Sub-Gaussianas e Desigualdade de Hoeffding Suponha que desejamos aplicar inclinação exponencial em importance sampling, mas a função log-partição ψ(λ) = logE[eλX] é desconhecida ou muito difícil de calcular exatamente. Nesse caso, muitas vezes basta conhecer um limite superior para ψ(λ). Se pudermos encontrar uma função simples que domina a logpartição verdadeira, ainda podemos controlar probabilidades de eventos raros e obter limites exponenciais. Por exemplo, suponha que sabemos que ψ(λ)≤σ2λ2 2,∀λ∈R. Isso significa que os momentos exponenciais de Xcrescem no máximo como os de uma variável normal com variância σ2. Dizemos então que Xésub-Gaussiana com parâmetro de variância σ2. Recorde que, sob inclinação exponencial, a probabilidade de um evento raro pode ser escrita como P(X≥a) = Eλhe−λX+ψ(λ)1{X≥a}i=eψ(λ)Eλhe−λX1{X≥a}i. Se a função log-partição exata é desconhecida, podemos substituí-la por qualquer limite superior válido. Usando a condição sub-Gaussiana acima, o argumento de Chernoff fornece P(X≥a)≤inf λ>0e−λa+ψ(λ)≤inf λ>0e−λa+σ2λ2 2. Minimizando o expoente em relação a λ, obtemos λ⋆=a/σ2, o que dá P(X≥a)≤exp−a2 2σ2. Portanto, qualquer variável aleatória cuja função log-partição é limitada por uma função quadrática apresenta cauda do tipo Gaussiana. Mesmo que não possamos realizar importance sampling exato sem conhecer a constante de normalização eψ(λ), a desigualdade acima fornece uma estimativa analítica precisa da probabilidade de evento raro. Como vimos na seção anterior, variáveis limitadas satisfazem um limite uniforme na variância inclinada, e, portanto, sua função log-partição cresce no máximo quadraticamente. Isso significa que qualquer variável limitada é automaticamente sub-Gaussiana, com parâmetro σ2=(b−a)2 4. Essa observação leva diretamente a um dos resultados mais fundamentais na teoria das desigualdades de concentração, conhecido como Lema de Hoeffding, que formaliza esse fato e fornece limites exponenciais explícitos para variáveis limitadas. Teorema 6 (Lema de Hoeffding).Seja X uma variável aleatória tal que X ∈[a,b]eE[X]=0. Então, para todo λ∈R, ψ(λ) = logEheλXi≤λ2(b−a)2 8. Em particular, X é sub-Gaussiana com parâmetro de variância σ2= (b−a)2/4. 6.5. POR QUE A INCLINAÇÃO EXPONENCIAL? 89 O lema de Hoeffding implica imediatamente um limite exponencial para a soma de variáveis limitadas independentes. Teorema 7 (Desigualdade de Hoeffding).Sejam X1, . . . , Xnvariáveis aleatórias independentes tais que Xi∈[ai,bi]eE[Xi]=0para todo i. Então, para qualquer t >0, P n ∑ i=1 Xi≥t!≤exp−2t2 ∑n i=1(bi−ai)2. Demonstração. Pelo lema de Hoeffding, cada Xisatisfaz EheλXii≤expλ2(bi−ai)2 8. Como os Xisão independentes, Eheλ∑n i=1Xii= n ∏ i=1 EheλXii≤exp λ2 8 n ∑ i=1 (bi−ai)2!. Aplicando o limite de Chernoff, P n ∑ i=1 Xi≥t!≤inf λ>0exp −λt+λ2 8 n ∑ i=1 (bi−ai)2!. Minimizando o expoente em relação a λ, obtemos λ⋆=4t ∑n i=1(bi−ai)2, e substituindo esse valor, P n ∑ i=1 Xi≥t!≤exp−2t2 ∑n i=1(bi−ai)2, o que conclui a prova. Essa desigualdade mostra que somas de variáveis aleatórias limitadas independentes exibem concentração do tipo Gaussiana: suas caudas decaem tão rapidamente quanto e−ct2, com uma constante determinada apenas pela largura dos intervalos [ai,bi]. No contexto de importance sampling, isso significa que, quando cada componente do estimador é limitado, o estimador como um todo permanece estável — a variância efetiva da medida inclinada não pode explodir. 6.5 Por que a Inclinação Exponencial? Ao realizar importance sampling, pode surgir a pergunta: por que usar a inclinação exponencial em vez de qualquer outra distribuição com a mesma média? Afinal, muitas reponderações podem satisfazer Eg[X]=a. O que torna a inclinação exponencial especial? Para entender isso, fazemos um breve desvio pela dualidade convexa. Dada uma função convexa ψ:Rd→R∪{+∞}, seu conjugado convexo (ou dual de Fenchel) é definido por ψ∗(y) = sup x∈Rd{⟨y,x⟩−ψ(x)}. 96 CAPÍTULO 7. CADEIAS DE MARKOV E MCMC Assim, q(5) 13 =52 243. A matriz de transição Qcodifica a distribuição condicional de X1dado o estado inicial da cadeia. Especificamente, a i-ésima linha de Qé a PMF condicional de X1dado X0=i, escrita como um vetor linha. De forma análoga, a i-ésima linha de Qncorresponde à PMF condicional de Xndado X0=i. Para obter as distribuições marginais de X0,X1, . . ., precisamos especificar não apenas a matriz de transição, mas também as condições iniciais da cadeia. O estado inicial X0pode ser especificado de forma determinística, ou de forma aleatória segundo alguma distribuição. Seja t= (t1,t2, . . . , tM)a PMF de X0, vista como vetor linha, em que ti=P(X0=i). Proposição 1 (Distribuição marginal de Xn).Seja t = (t1,t2, . . . , tM)o vetor de probabilidades iniciais, com ti=P(X0=i). Então a distribuição marginal de Xné dada por tQn. Em particular, a j-ésima componente do vetor tQné P(Xn=j). Demonstração. Pela lei da probabilidade total, condicionando em X0, a probabilidade de a cadeia estar no estado japós npassos é P(Xn=j) = M ∑ i=1 P(X0=i)P(Xn=j|X0=i) = M ∑ i=1 tiq(n) ij . Mas essa soma corresponde exatamente à j-ésima componente do vetor tQn, pela definição de multiplicação de matrizes. Exemplo 27 (Distribuições marginais de uma cadeia de Markov com 4 estados).Considere novamente a cadeia de Markov com 4 estados representada na Figura anterior. Suponha que as condições iniciais sejam dadas por t=1 4,1 4,1 4,1 4, isto é, a cadeia começa com igual probabilidade em cada um dos quatro estados. Seja Xna posição da cadeia no tempo n. A distribuição marginal de X1é tQ =1 41 41 41 4      1/3 1/3 1/3 0 0 0 1/2 1/2 0 1 0 0 0 1/2 0 0 1/2      =5 24,1 3,5 24,1 4. A distribuição marginal de X5é tQ5=1 41 41 41 4      853/3888 509/1944 52/243 395/1296 173/864 85/432 31/108 91/288 37/144 29/72 1/9 11/48 499/2592 395/1296 71/324 245/864      =3379 15552,2267 7776,101 486,1469 5184. Neste caso utilizamos o computador para realizar as multiplicações de matrizes. 7.1. CADEIAS DE MARKOV (RESUMO) 97 7.1.1 Classificação dos estados Nesta parte introduziremos a terminologia usada para descrever as várias características de uma cadeia de Markov. Os estados de uma cadeia podem ser classificados como recorrentes ou transientes, dependendo de o processo retornar ou não a eles ao longo do tempo. Além disso, cada estado possui um período, que é um número inteiro positivo que resume a quantidade de passos que pode decorrer entre visitas sucessivas a esse estado. Essas características são importantes porque determinam o comportamento de longo prazo da cadeia de Markov, que será estudado mais adiante. Os conceitos de recorrência e transiência são melhor ilustrados com um exemplo. Na cadeia de Markov mostrada à esquerda da Figura, uma partícula se movendo entre os estados continuará visitando todos os quatro estados indefinidamente, pois é possível transitar de qualquer estado para qualquer outro. Em contraste, considere a cadeia à direita da Figura, e suponha que a partícula comece no estado 1. Durante algum tempo, a cadeia pode permanecer no triângulo formado pelos estados 1,2,3, mas eventualmente atingirá o estado 4. A partir do momento em que entra no estado 4, a cadeia nunca mais retorna a 1, 2,3, e passa a se mover apenas entre os estados 4,5,6 para sempre. Assim, os estados 1,2,3 são transientes, enquanto os estados 4,5,6 são recorrentes. Definição 5 (Estados recorrentes e transientes).Um estado i de uma cadeia de Markov é dito recorrente se, partindo de i, a probabilidade de que a cadeia eventualmente retorne a i é igual a 1. Caso contrário, o estado é dito transiente, o que significa que, se a cadeia começar em i, existe probabilidade positiva de nunca mais retornar a i. Embora a definição de estado transiente apenas exija que haja probabilidade positiva de nunca retornar ao estado, podemos dizer algo mais forte: sempre que existir probabilidade positiva de abandonar ipara sempre, a cadeia inevitavelmente deixará o estado iem algum momento. Além disso, é possível caracterizar a distribuição do número de retornos ao estado. Proposição 2 (Número de retornos a um estado transiente é Geométrico).Seja i um estado transiente de uma cadeia de Markov. Suponha que a probabilidade de nunca retornar a i, partindo de i, seja 98 CAPÍTULO 7. CADEIAS DE MARKOV E MCMC p>0. Então, partindo de i, o número de vezes que a cadeia retorna a i antes de sair para sempre é uma variável aleatória com distribuição Geom(p). Demonstração. A demonstração segue pela interpretação da distribuição Geométrica. Cada vez que a cadeia está em i, temos um ensaio de Bernoulli: ocorre “sucesso” se a cadeia sair de ipara sempre, e ocorre “falha” se a cadeia eventualmente retornar a i. Esses ensaios são independentes pela propriedade de Markov. O número de retornos ao estado icorresponde ao número de falhas antes do primeiro sucesso, exatamente a história da distribuição Geométrica. E como uma variável Geométrica assume valores finitos com probabilidade 1, concluímos que, após um número finito de visitas, a cadeia deixará o estado ipara sempre. Se o número de estados não for muito grande, uma maneira de classificar estados como recorrentes ou transientes é desenhar o diagrama da cadeia de Markov e aplicar o mesmo tipo de raciocínio feito na análise dos exemplos anteriores. Um caso especial em que podemos concluir imediatamente que todos os estados são recorrentes ocorre quando a cadeia é irredutível, isto é, quando é possível ir de qualquer estado a qualquer outro. Definição 6 (Cadeia irredutível e redutível).Uma cadeia de Markov com matriz de transição Q é dita irredutível se, para quaisquer dois estados i e j, for possível ir de i até j em um número finito de passos, com probabilidade positiva. Isto é, para quaisquer estados i,j, existe um número inteiro n >0tal que a entrada (i,j)de Qné positiva. Uma cadeia que não é irredutível é chamada de redutível. Proposição 3 (Irredutibilidade implica recorrência de todos os estados).Em uma cadeia de Markov irredutível com espaço de estados finito, todos os estados são recorrentes. Demonstração. É claro que pelo menos um estado deve ser recorrente; se todos fossem transientes, a cadeia eventualmente abandonaria todos os estados para sempre, o que é impossível. Sem perda de generalidade, suponha que o estado 1 seja recorrente. Considere outro estado i. Pela definição de irredutibilidade, existe algum ntal que q(n) 1i>0. Assim, toda vez que a cadeia visita o estado 1, há uma probabilidade positiva de que, após npassos, ela esteja no estado i. Como a cadeia visita o estado 1 infinitas vezes (por recorrência), concluirá inevitavelmente no estado i. Além disso, partindo de i, a cadeia retorna ao estado 1, já que este é recorrente. Aplicando o mesmo argumento recursivamente, a cadeia visitará o estado iinfinitas vezes. Como ifoi arbitrário, segue que todos os estados são recorrentes. A recíproca da proposição anterior é falsa: é possível ter uma cadeia de Markov redutível em que todos os estados sejam recorrentes. Um exemplo é a cadeia ilustrada na Figura abaixo, que consiste em duas “ilhas” de estados. 7.1. CADEIAS DE MARKOV (RESUMO) 99 Outra forma de classificar estados é de acordo com seus períodos. O período de um estado resume quanto tempo pode se passar entre visitas sucessivas a esse estado. Definição 7 (Período de um estado, cadeia periódica e aperiódica).Operíodo de um estado i em uma cadeia de Markov é o máximo divisor comum (mdc) dos números de passos em que é possível retornar a i, partindo de i. Mais precisamente, o período de i é d(i) = gcd{n≥1 : (Qn)ii >0}. Se nunca for possível retornar a i, definimos d(i) = ∞. Um estado é dito aperiódico se d(i) = 1, e periódico caso contrário. A cadeia como um todo é chamada aperiódica se todos os seus estados forem aperiódicos, e periódica caso contrário. Exemplo 28 (Periodicidade em cadeias de Markov).Considere novamente as duas cadeias de Markov da Figura abaixo. Na cadeia com 6 estados (à direita), partindo do estado 1, é possível retornar a ele após 3passos, 6 passos, 9passos e assim por diante. No entanto, não é possível retornar ao estado 1em um número de passos que não seja múltiplo de 3. Portanto, o estado 1tem período 3. De forma análoga, os estados 2e3 também têm período 3. Por outro lado, os estados 4,5, 6 possuem período 1. Como nem todos os estados têm período 1, a cadeia é considerada periódica. Em contraste, na cadeia da Figura 2 (à esquerda), todos os estados são aperiódicos, de modo que a cadeia como um todo é aperiódica. 100 CAPÍTULO 7. CADEIAS DE MARKOV E MCMC IMPORTANTE!!! Vale destacar uma diferença importante: algumas propriedades pertencem aestados individuais da cadeia de Markov, enquanto outras pertencem à cadeia como um todo. A tabela abaixo resume essa distinção, destacando os conceitos principais lado a lado para facilitar a comparação. Propriedades de estados Propriedades da cadeia Estado recorrente: Partindo de i, a probabilidade de eventualmente retornar a ié 1. Cadeia irredutível: É possível ir de qualquer estado ipara qualquer estado jem um número finito de passos, com probabilidade positiva. Estado transiente: Partindo de i, existe probabilidade positiva de nunca mais retornar ai. Cadeia redutível: Não é possível ir de alguns estados ipara outros jem um número finito de passos, com probabilidade positiva. Estado aperiódico: O período do estado ié 1, ou seja, é possível retornar a iem tempos arbitrários suficientemente grandes. Cadeia aperiódica: Todos os estados da cadeia são aperiódicos. Estado periódico: O período do estado ié maior que 1, ou seja, o retorno a isó pode ocorrer em múltiplos de algum inteiro d> 1. Cadeia periódica: Pelo menos um estado da cadeia é periódico. 7.1.2 Distribuição estacionária Os conceitos de recorrência e transiência são fundamentais para compreender o comportamento de longo prazo de uma cadeia de Markov. Inicialmente, a cadeia pode passar algum tempo em estados transitórios; porém, com o tempo, ela tende a permanecer apenas nos estados recorrentes. Surge então uma pergunta natural: qual fração do tempo a cadeia passará em cada um desses estados recorrentes? A resposta é dada pela distribuição estacionária da cadeia, também chamada de distribuição em regime permanente. Veremos nesta seção que, para cadeias de Markov irredutíveis e aperiódicas, a distribuição estacionária descreve o comportamento assintótico da cadeia, independentemente das condições iniciais. Ela fornece tanto a probabilidade de longo prazo de estar em um determinado estado quanto a proporção de tempo que a cadeia passará nesse estado. Definição 8 (Distribuição estacionária).Um vetor linha s = (s1, . . . , sM), com si≥0e∑isi=1, é dito uma distribuição estacionária para uma cadeia de Markov com matriz de transição Q se ∑ i siqij =sj,para todo j. Esse sistema de equações lineares pode ser escrito de forma compacta como sQ =s. 7.1. CADEIAS DE MARKOV (RESUMO) 101 Recorde que, se sé a distribuição de X0, então sQ é a distribuição marginal de X1. Assim, a igualdade sQ =ssignifica que, se X0tem distribuição s, então X1também terá distribuição s. Pelo mesmo raciocínio, X2,X3, . . . também seguirão a mesma distribuição. Em outras palavras, uma cadeia de Markov cuja distribuição inicial é a distribuição estacionária spermanecerá nessa distribuição para sempre. Observação 2. Podemos ter uma interpretação intuitiva da distribuição estacionária a partir de uma simulação mental. Imagine que temos um número muito grande de partículas (por exemplo, um bilhão), e que a distribuição inicial dessas partículas entre os estados é proporcional à distribuição inicial da cadeia. Em seguida, fazemos todas as partículas evoluírem segundo a matriz de transição Q. Após um certo número de passos n, contamos quantas partículas estão em cada estado. Quando a cadeia atinge o regime estacionário, essas proporções se estabilizam: se contarmos novamente após aplicar Q mais uma vez, as frações relativas de partículas em cada estado permanecerão essencialmente as mesmas. Assim, a distribuição estacionária s representa justamente essa configuração de equilíbrio em que a aplicação de Q não altera mais as proporções — isto é, sQ =s. Exemplo 29 (Distribuição estacionária para uma cadeia com dois estados).Considere a matriz de transição Q= 1 32 3 1 21 2!. A distribuição estacionária é da forma s = (s,1 −s). Devemos então resolver (s,1 −s) 1 32 3 1 21 2!= (s,1 −s), o que é equivalente ao sistema      1 3s+1 2(1−s) = s, 2 3s+1 2(1−s) = 1−s. A única solução é s =3 7. Portanto, a distribuição estacionária única dessa cadeia de Markov é s=3 7,4 7. Mais geralmente, suponha que q12 =aeq21 =b, com 0<a<1e0<b<1. A matriz de transição é então Q= 1−a a b1−b!. Escrevendo s = (s1,s2), a equação sQ =s fornece o sistema linear      (1−a)s1+bs2=s1, as1+ (1−b)s2=s2. Ambas as equações se reduzem a as1=bs2. Como s2=1−s1, obtemos a solução única s=b a+b,a a+b. 102 CAPÍTULO 7. CADEIAS DE MARKOV E MCMC Existência, unicidade e convergência Surge naturalmente a questão: uma distribuição estacionária sempre existe? E, caso exista, ela é única? Para cadeias de Markov com espaço de estados finito, a resposta é afirmativa: sempre existe uma distribuição estacionária. Além disso, quando a cadeia é irredutível, essa distribuição é única. Teorema 8 (Existência e unicidade da distribuição estacionária).Para qualquer cadeia de Markov irredutível, existe uma única distribuição estacionária. Nessa distribuição, todos os estados possuem probabilidade positiva. Esse resultado decorre de um teorema clássico da álgebra linear conhecido como teorema de Perron–Frobenius. Além da existência e unicidade, também é importante entender a convergência para a distribuição estacionária. Já afirmamos de forma informal que a distribuição estacionária descreve o comportamento de longo prazo da cadeia: se a cadeia for executada por tempo suficiente, a distribuição marginal de Xntende à distribuição estacionária s. O resultado a seguir formaliza essa ideia. Teorema 9 (Convergência para a distribuição estacionária).Se (X0,X1, . . .)é uma cadeia de Markov irredutível e aperiódica com distribuição estacionária s e matriz de transição Q, então P(Xn=i)−→ siquando n →∞. Em termos matriciais, Qnconverge para uma matriz cujas linhas são todas iguais a s. Portanto, após um número suficientemente grande de passos, a probabilidade de que a cadeia esteja em um estado ise aproxima de si, independentemente das condições iniciais. Isso mostra que cadeias irredutíveis e aperiódicas são particularmente agradáveis de trabalhar, pois o seu comportamento assintótico é estável e previsível. Observação 3. De modo intuitivo, a condição adicional de aperiodicidade serve para evitar cadeias que apenas “giram em ciclos”, alternando de forma determinística entre grupos de estados. Por exemplo, em cadeias onde certos estados só são acessíveis após um número par de passos, enquanto outros apenas após um número ímpar, a convergência não ocorre sem essa hipótese. A combinação de irredutibilidade e aperiodicidade garante que a cadeia possa se misturar completamente no espaço de estados e, portanto, converge para sua distribuição estacionária. Exemplo 30 (Cadeia periódica).Considere a cadeia de Markov ilustrada abaixo, em que cada estado possui período 5. Q=         01000 00100 00010 00001 10000         . 7.1. CADEIAS DE MARKOV (RESUMO) 103 Pode-se verificar facilmente que s=1 5,1 5,1 5,1 5,1 5 é uma distribuição estacionária dessa cadeia, e que ela é única. No entanto, suponha que a cadeia comece em X0=1. Nesse caso, a distribuição de Xnatribui probabilidade 1ao estado (nmod 5) + 1e probabilidade 0a todos os outros estados. Em particular, a distribuição de Xnnão converge para s quando n →∞. De forma equivalente, a matriz Qnnão converge para uma matriz em que todas as linhas são iguais a s. As transições dessa cadeia são determinísticas, e portanto cada Qncontinua sendo uma matriz composta apenas por zeros e uns. Esse exemplo mostra que, embora a distribuição estacionária exista e seja única, a convergência não ocorre quando a cadeia é periódica. Observação 4. A condição de irredutibilidade é essencial para que a distribuição estacionária seja única e represente o comportamento de longo prazo da cadeia. Se a cadeia não for irredutível, o espaço de estados pode se decompor em vários subconjuntos fechados — ou seja, conjuntos de estados dos quais não é possível sair. Cada um desses subconjuntos pode ter a sua própria distribuição estacionária, o que implica que não há uma única distribuição que descreva o comportamento de toda a cadeia. Por exemplo, suponha que existam dois conjuntos de estados A e B tais que, uma vez que a cadeia entra em A, nunca mais pode ir para B, e vice-versa. Então, a probabilidade de longo prazo de estar em A ou em B dependerá da condição inicial. Isso impede a convergência para uma distribuição estacionária única. A irredutibilidade elimina esse problema: ela garante que todos os estados se comunicam entre si, isto é, para quaisquer i e j, existe algum número de passos n tal que (Qn)ij >0. Com isso, a cadeia pode eventualmente alcançar qualquer estado a partir de qualquer outro, o que assegura tanto a existência quanto a unicidade da distribuição estacionária e a convergência para ela. Tempo médio de retorno e comportamento de longo prazo Além de descrever o comportamento assintótico da cadeia, a distribuição estacionária também está relacionada ao tempo médio entre visitas a um estado específico. 104 CAPÍTULO 7. CADEIAS DE MARKOV E MCMC Teorema 10 (Tempo esperado de retorno).Considere uma cadeia de Markov irredutível com distribuição estacionária s. Seja rio tempo esperado para a cadeia retornar ao estado i, dado que ela começa em i. Então, si=1 ri. Esse resultado mostra que estados com maior probabilidade estacionária são visitados com mais frequência — em média, o tempo até retornar a eles é menor. Exemplo 31 (Comportamento de longo prazo de uma cadeia com dois estados).Considere novamente a cadeia de dois estados discutida anteriormente, cuja matriz de transição é Q= 1 32 3 1 21 2!. A distribuição estacionária é s=3 7,4 7. No longo prazo, a cadeia passará aproximadamente 3/7 do tempo no estado 1 e 4/7 do tempo no estado 2. Começando no estado 1, o tempo médio para retornar a esse estado é r1=7/3, em conformidade com o teorema acima, pois s1=1/r1. Além disso, as potências da matriz de transição convergem para uma matriz em que cada linha coincide com a distribuição estacionária: Qn= 1 32 3 1 21 2!n −→ 3 74 7 3 74 7!,quando n →∞. 7.1.3 Reversibilidade Vimos que a distribuição estacionária de uma cadeia de Markov é extremamente útil para compreender seu comportamento de longo prazo. Entretanto, em muitos casos pode ser computacionalmente difícil determinar essa distribuição, especialmente quando o espaço de estados é grande. Nesta seção, estudamos um caso especial importante em que é possível evitar o cálculo direto das equações de autovalor associadas à matriz de transição. Definição 9 (Reversibilidade).Seja Q = (qij)a matriz de transição de uma cadeia de Markov. Dizemos que a cadeia é reversível em relação a um vetor s = (s1, . . . , sM), com si≥0e∑isi=1, se siqij =sjqji,para todos os estados i,j. Essa equação é chamada de condição de equilíbrio detalhado (ou detailed balance condition). A intuição por trás da reversibilidade é a seguinte: uma cadeia reversível, iniciada segundo sua distribuição estacionária, se comporta da mesma forma independentemente de o tempo estar sendo observado para frente ou para trás. Mais precisamente, quando a cadeia está em equilíbrio, a probabilidade de sair do estado ie ir para o estado jem um passo é siqij, e essa probabilidade é exatamente igual à de sair de je voltar para i, que é sjqji. Em outras palavras, o fluxo de probabilidade de ipara jé o mesmo que o de jpara i: siqij =sjqji. 7.1. CADEIAS DE MARKOV (RESUMO) 105 Isso significa que, no regime estacionário, as transições “para frente” e “para trás” ocorrem com a mesma frequência média, de modo que, se observarmos a cadeia no tempo inverso, ela parecerá estatisticamente idêntica à original. Outra maneira de entender a reversibilidade é pensar na cadeia como um sistema com um grande número de partículas que se movem de forma independente de acordo com as probabilidades de transição. No longo prazo, a proporção de partículas em cada estado jé dada pela probabilidade estacionária sj. O equilíbrio estacionário garante que, em média, o fluxo de partículas que sai de cada estado é igual ao fluxo de partículas que entra nele. Mais precisamente, seja no número total de partículas e so vetor de proporções atuais de partículas em cada estado. Temos que sé estacionário se, e somente se, sj=∑ i siqij =sjqjj +∑ i=j siqij, para todo j. Multiplicando por n, obtemos nsj(1−qjj) = ∑ i=j nsiqij. O lado esquerdo representa o número médio de partículas que sairão do estado jno próximo passo, enquanto o lado direito representa o número médio de partículas que entrarão em j. Portanto, há um equilíbrio entre entrada e saída de partículas em cada estado. A condição de reversibilidade impõe uma forma ainda mais forte de equilíbrio: para cada par de estados distintos iej, nsiqij =nsjqji. O lado esquerdo corresponde ao número médio de partículas que vão de ipara j, e o lado direito ao número médio que vai de jpara i. Assim, a reversibilidade garante que, par a par, o fluxo entre dois estados é perfeitamente equilibrado. Proposição 4 (Reversibilidade implica estacionariedade).Se Q = (qij)é a matriz de transição de uma cadeia de Markov reversível em relação a um vetor s = (s1, . . . , sM)não negativo, com soma dos componentes igual a 1, então s é uma distribuição estacionária da cadeia. Demonstração. Temos ∑ i siqij =∑ i sjqji =sj∑ i qji =sj, onde a última igualdade decorre do fato de que a soma das probabilidades em cada linha de Q é igual a 1. Logo, sQ =s, e portanto sé estacionária. Esse é um resultado poderoso, pois frequentemente é mais simples verificar a condição de reversibilidade do que resolver o sistema completo de equações sQ =s. Um caso importante e simples de cadeia reversível ocorre quando a matriz de transição Qé simétrica. Se Qé simétrica, isto é, qij =qji para todos os i,j, então a distribuição estacionária é uniforme sobre o espaço de estados: s=1 M,1 M, . . . , 1 M. 112 CAPÍTULO 7. CADEIAS DE MARKOV E MCMC Após um número suficientemente grande de iterações, as variáveis Wn,Wn+1, . . . seguem aproximadamente a distribuição Beta(a,b). Essas amostras, entretanto, não são independentes: a cadeia gera uma sequência de variáveis correlacionadas, que oscilam em torno da distribuição alvo conforme o processo evolui. Perceba que no algoritmo geral de Metropolis–Hastings, a probabilidade de aceitação é dada por a(x,u) = min1, s(u)p(x|u) s(x)p(u|x), onde p(u|x)é a densidade da proposta, isto é, a probabilidade de propor o ponto u dado o estado atual x. O termo p(x|u) p(u|x)está presente para corrigir possíveis assimetrias da distribuição de proposta. Por exemplo, se é mais fácil propor de x para u do que o contrário, o fator p(x|u) p(u|x)compensa essa diferença, garantindo que o fluxo médio de partículas entre x e u permaneça equilibrado: s(x)p(u|x)a(x,u) = s(u)p(x|u)a(u,x). No caso da simulação da distribuição Beta, utilizamos um independence sampler, em que as propostas são independentes do estado atual: p(u|x) = p(u) = Unif(0,1). Consequentemente, p(x|u) = p(x) = Unif(0,1),e portanto p(x|u) p(u|x)=1. O fator de correção se cancela, e a fórmula de aceitação se reduz para a(x,u) = min1, s(u) s(x), que é exatamente a expressão usada no caso da Beta. Em resumo, o termo da proposta desaparece porque a distribuição de proposta é simétrica e independente do estado atual, de modo que não há assimetria a corrigir. Observação 6 (Período de aquecimento (burn-in)).Mesmo quando uma cadeia de Markov possui s como distribuição estacionária, isso significa apenas que s é um estado de equilíbrio: se a cadeia for iniciada com distribuição s, ela permanecerá em s para sempre. 7.4. AMOSTRAGEM DE GIBBS 113 Na prática, porém, a cadeia é iniciada em um ponto fixo X0=i ou segundo alguma distribuição inicial diferente de s. As primeiras iterações servem para que a cadeia se desloque gradualmente em direção ao equilíbrio, aproximando-se da distribuição estacionária. Durante esse período inicial, as distribuições de Xnainda refletem o estado inicial e não representam bem o comportamento estacionário. Esse intervalo é chamado de período de aquecimento, ou burn-in. Se denotarmos por p(n)o vetor de probabilidades no tempo n, temos p(n+1)=p(n)Q. Conforme n cresce, ocorre a convergência p(n)−→ s, isto é, a distribuição da cadeia tende à distribuição estacionária s. Somente após essa fase de convergência é que as amostras podem ser consideradas representativas do regime estacionário. Intuitivamente, podemos imaginar muitas cópias da cadeia evoluindo em paralelo. Inicialmente, há um acúmulo de partículas em certos estados e um déficit em outros. À medida que o tempo passa, os fluxos de transição entre estados se equilibram, até que a proporção de partículas em cada estado se estabilize segundo s. Descartar as primeiras amostras equivale a ignorar essa fase de ajuste até o equilíbrio. 7.4 Amostragem de Gibbs Aamostragem de Gibbs é um algoritmo de Monte Carlo empregado para gerar amostras aproximadas de uma distribuição conjunta. A ideia central é simples: atualizar sucessivamente uma variável de cada vez, amostrando-a de sua distribuição condicional dado o valor atual das demais. Esse método é particularmente útil quando as distribuições condicionais são fáceis de manipular e de amostrar diretamente. Considere o caso de duas variáveis aleatórias discretas XeY, com função de probabilidade conjunta pX,Y(x,y) = P(X=x,Y=y). Desejamos construir uma cadeia de Markov (Xn,Yn)cuja distribuição estacionária seja pX,Y. Existem duas versões principais do algoritmo de Gibbs, dependendo de como as variáveis são atualizadas: (1) o Gibbs sistemático, no qual as variáveis são atualizadas em ordem fixa e alternada; e (2) o Gibbs aleatório, no qual a variável a ser atualizada é escolhida aleatoriamente a cada iteração. O procedimento pode ser descrito da seguinte forma: 1. No passo atual, suponha que a cadeia esteja em (Xn,Yn) = (xn,yn). 2. Gere um novo valor xn+1a partir da distribuição condicional de Xdado Y=yn, isto é, xn+1∼p(x|Y=yn). 3. Em seguida, gere um novo valor yn+1a partir da distribuição condicional de Ydado X= xn+1: yn+1∼p(y|X=xn+1). 114 CAPÍTULO 7. CADEIAS DE MARKOV E MCMC 4. Atualize o estado da cadeia para (Xn+1,Yn+1) = (xn+1,yn+1). Repetindo esses passos indefinidamente, a cadeia (Xn,Yn)converge para a distribuição estacionária pX,Y. Na versão aleatória, escolhe-se a cada iteração qual variável será atualizada, com probabilidades iguais. O procedimento segue: 1. No passo atual, suponha que a cadeia esteja em (Xn,Yn) = (xn,yn). 2. Escolha aleatoriamente qual componente será atualizado: • com probabilidade 1/2, atualize X; • com probabilidade 1/2, atualize Y. 3. Se Xfor escolhido: (a) Gere xn+1∼p(x|Y=yn); (b) Defina (Xn+1,Yn+1) = (xn+1,yn). 4. Caso Yseja escolhido: (a) Gere yn+1∼p(y|X=xn); (b) Defina (Xn+1,Yn+1) = (xn,yn+1). Repetindo o processo, obtemos novamente uma cadeia cuja distribuição estacionária é pX,Y. O algoritmo de Gibbs se estende naturalmente para o caso de dvariáveis aleatórias. Nesse caso, o estado da cadeia é um vetor Wn= (W(1) n, . . . ,W(d) n). Em cada iteração: 1. Escolhe-se (de forma determinística ou aleatória) um índice j∈ {1, . . . , d}; 2. Amostra-se o componente W(j) nda distribuição condicional W(j) n+1∼pw(j)|w(−j), onde w(−j)denota todos os outros componentes fixos; 3. Mantêm-se os demais componentes inalterados. Observação 7. O amostrador de Gibbs pode ser interpretado como um caso especial do algoritmo de Metropolis–Hastings. Enquanto o Metropolis–Hastings requer uma distribuição proposta e uma etapa de aceitação ou rejeição, o Gibbs utiliza propostas que sempre são aceitas, pois cada amostra é retirada exatamente da distribuição condicional correta. Teorema 11 (Gibbs aleatório como caso particular de Metropolis–Hastings).O amostrador de Gibbs aleatório é um caso particular do algoritmo de Metropolis–Hastings, no qual toda proposta é sempre aceita. Em particular, isso implica que a distribuição estacionária do amostrador de Gibbs aleatório é exatamente a distribuição conjunta desejada. 7.4. AMOSTRAGEM DE GIBBS 115 Demonstração. Apresentaremos a demonstração no caso bidimensional, embora o argumento se estenda naturalmente a qualquer número de dimensões. Sejam XeYvariáveis aleatórias discretas cuja distribuição conjunta p(x,y) = P(X=x,Y=y) é a distribuição estacionária que queremos obter. O algoritmo de Metropolis–Hastings, no estado atual (x,y), procede da seguinte forma: propõe um novo estado (x′,y′)segundo uma distribuição proposta q((x′,y′)|(x,y)) e aceita essa proposta com probabilidade a((x,y),(x′,y′)) = min1, p(x′,y′)q((x,y)|(x′,y′)) p(x,y)q((x′,y′)|(x,y)) . No caso do Gibbs aleatório, a proposta consiste em escolher aleatoriamente uma das coordenadas e atualizá-la de acordo com sua distribuição condicional verdadeira. Mais precisamente: 1. Com probabilidade 1/2, atualiza-se Xa partir de p(x′|Y=y), mantendo Y′=y. 2. Com probabilidade 1/2, atualiza-se Ya partir de p(y′|X=x), mantendo X′=x. Vamos considerar o segundo caso, em que apenas Yé atualizado (o caso simétrico em que X é atualizado é análogo). Assim, q((x,y′)|(x,y)) = 1 2p(y′|x)eq((x,y)|(x,y′)) = 1 2p(y|x). Substituindo esses termos na fórmula de aceitação, obtemos: a((x,y),(x,y′)) = min1, p(x,y′)p(y|x) p(x,y)p(y′|x). Como p(x,y) = p(x)p(y|x), temos: p(x,y′)p(y|x) p(x,y)p(y′|x)=p(x)p(y′|x)p(y|x) p(x)p(y|x)p(y′|x)=1. Logo, a((x,y),(x,y′)) = 1. Isto é, toda proposta é sempre aceita. Consequentemente, o algoritmo de Metropolis–Hastings, com essa escolha específica de distribuição proposta, coincide exatamente com o amostrador de Gibbs aleatório, e ambos têm a mesma distribuição estacionária p(x,y). Exemplo 34 (O problema da galinha e dos ovos).Uma galinha põe um número N de ovos, onde N∼Poisson(λ). Cada ovo choca com probabilidade p, onde p é desconhecido e tem distribuição p∼Beta(a,b). Os parâmetros λ,a,b são conhecidos. O problema é que não observamos o número total de ovos N, apenas o número de ovos chocados, denotado por X. Nosso objetivo é estimar a esperança posterior E[p|X=x], isto é, a média de p após observar x ovos chocados. 116 CAPÍTULO 7. CADEIAS DE MARKOV E MCMC 1. A distribuição de X dado p é (pense no porquê) X|p∼Poisson(λp). Logo, a densidade posterior de p é proporcional a f(p|X=x)∝P(X=x|p)f(p)∝e−λp(λp)xpa−1(1−p)b−1. Essa distribuição não tem forma fechada conhecida, o que dificulta a amostragem direta. Para contornar isso, introduzimos a variável latente N, correspondente ao número total de ovos postos. 2. Condicionalmente a N =n e p, o número de ovos chocados segue X|N=n,p∼Binomial(n,p). Assim, ao condicionar em N, recuperamos a conjugação Beta–Binomial: p|X=x,N=n∼Beta(x+a,n−x+b). O fato de que essa forma condicional é simples motiva o uso da amostragem de Gibbs: alternaremos entre amostrar p dado N e N dado p. 3. Desejamos gerar amostras da distribuição conjunta de (p,N)condicionada a X =x. O algoritmo segue os seguintes passos: (a) Faça suposições iniciais para p e N. (b) Repita até a convergência: i. Atualização de p: p|X=x,N=n∼Beta(x+a,n−x+b). ii. Atualização de N:Seja Y =N−X o número de ovos que não chocaram. Condicionalmente a p, o número de ovos não chocados segue Y|p∼Poisson(λ(1−p)). Assim, sorteamos Y ∼Poisson(λ(1−p)) e definimos N=X+Y. 4. Após muitas iterações, obtemos amostras (p(1),N(1)),(p(2),N(2)), . . . extraídas aproximadamente de f(p,N|X=x). A esperança posterior é então estimada pela média amostral: E(p|X=x)≈1 T T ∑ t=1 p(t). Os resultados obtidos por simulação do amostrador de Gibbs são resumidos na Tabela 4. Nas simulações, utilizou-se λ=10 como valor esperado do número total de ovos postos, e um prior Beta(a,b) com a =b=1, correspondendo a uma distribuição uniforme sobre [0,1]para o parâmetro p. Foram 7.4. AMOSTRAGEM DE GIBBS 117 considerados diferentes valores observados de ovos chocados X ∈ {3, 5,7, 9}, enquanto os demais parâmetros permaneceram fixos. XE[p|X=x]E[N|X=x] 3 0.40 9.0 5 0.56 9.3 7 0.69 10.1 9 0.77 11.3 Observa-se que, à medida que o número de ovos chocados X aumenta, a média posterior de p cresce de forma aproximadamente monotônica, variando de cerca de 0.4 para 0.77. Esse comportamento é exatamente o esperado, pois X |p∼Poisson(λp): quanto maior o número de ovos chocados, maior deve ser a probabilidade de sucesso p. Além disso, a média posterior de N também aumenta levemente com X, o que é coerente com o modelo N =X+Y, em que Y |p∼Poisson(λ(1−p)). Ou seja, observar mais ovos chocados implica, em média, um número total ligeiramente maior de ovos postos. Esses resultados confirmam que o amostrador de Gibbs captura corretamente a relação entre X, p e N, produzindo inferências consistentes com o modelo teórico. 118 CAPÍTULO 7. CADEIAS DE MARKOV E MCMC Capítulo 8 Processos de Difusão 8.1 Movimento Browniano e SDE A solução de uma equação diferencial estocástica (SDE, do inglês stochastic differential equation) é um processo estocástico. Um processo estocástico é uma coleção de variáveis aleatórias Xt∈Rd que evoluem ao longo do tempo t≥0. Embora cada Xtseja aleatória para um instante fixo, o interesse está em compreender como valores em tempos diferentes se relacionam — isto é, como Xt+sdepende de Xt. Um exemplo fundamental de processo estocástico é o movimento browniano, denotado por Wt. Ele é caracterizado por duas propriedades essenciais: •Incrementos normais: os incrementos têm distribuição normal com variância proporcional ao intervalo de tempo: Wt+s−Wt∼ N(0, sId), para todo t,s≥0. •Incrementos independentes: para quaisquer t1>t2>t3, os incrementos Wt1−Wt2e Wt2−Wt3são independentes. Essas propriedades fazem do movimento browniano o modelo canônico de ruído contínuo, servindo como base para a definição das equações diferenciais estocásticas. Podemos aproximar numericamente uma trajetória de movimento browniano discretizando o tempo. Se tomamos passos igualmente espaçados de tamanho s>0, o incremento em cada passo é dado por Wt+s=Wt+√sε, onde ε∼ N(0, Id)é uma variável aleatória independente a cada passo. O parâmetro srepresenta otamanho do passo temporal — isto é, o intervalo entre duas amostragens consecutivas do processo. O fator √sgarante que a variância dos incrementos cresça linearmente com o tempo, como exige a definição de movimento browniano: Var(Wt+s−Wt) = sId. Essa relação mostra que quanto maior o intervalo s, maior a variabilidade esperada do incremento, refletindo a natureza difusiva do processo. 119 120 CAPÍTULO 8. PROCESSOS DE DIFUSÃO Para ilustrar, vejamos como podemos simular numericamente uma trajetória de movimento browniano. A ideia é construir uma sequência de valores (W0,Ws,W2s, . . . , WT)que satisfaça as propriedades do processo. 1. Escolha dos parâmetros: Defina o tempo total de simulação T>0 e o número de passos n. O tamanho do passo será s=T/n. 2. Inicialização: Comece com W0=0, que é a condição inicial típica do movimento browniano. 3. Geração dos incrementos: Para cada passo k=1, 2, . . . , n, gere um ruído gaussiano independente εk∼ N(0,1), e compute o incremento ∆Wk=√sεk. 4. Construção da trajetória: Atualize o valor do processo de forma recursiva: Wk=Wk−1+∆Wk. O vetor (W0,W1, . . . , Wn)representa uma amostra discreta da trajetória de Wtno intervalo [0, T]. O resultado é mostrado na Figura 8.1, que exibe duas trajetórias com diferentes tamanhos de passo s. Quanto menor o passo, mais suave e precisa é a aproximação da trajetória contínua de um movimento browniano. Figura 8.1: Simulação de duas trajetórias de movimento browniano com diferentes tamanhos de passo s. Exercício 27. Simule movimentos brownianos com diferente passos. Com base nesse ruído contínuo Wt, uma equação diferencial estocástica é uma equação que descreve a evolução de um processo Xtsegundo dXt=f(Xt,t)dt +g(Xt,t)dWt. 8.1. MOVIMENTO BROWNIANO E SDE 121 O primeiro termo, f(Xt,t)dt, representa a tendência média do movimento e é chamado de drift; o segundo termo, g(Xt,t)dWt, modela a difusão, responsável pelas flutuações aleatórias. Para obter uma intuição mais concreta sobre a dinâmica dessa equação, é útil pensar em sua forma discretizada no tempo. Consideremos pequenos intervalos de tempo ∆>0. O movimento browniano Btpossui incrementos Bt+∆−Bt∼ N(0, ∆Id), e, além disso, esses incrementos são independentes para intervalos disjuntos. Assim, podemos representar o incremento como Bt+∆−Bt=√∆εt, onde εt∼ N(0, Id)é uma variável aleatória independente a cada passo. Substituindo esse termo na SDE dXt=f(Xt,t)dt +g(Xt,t)dBt, obtemos a versão discreta aproximada: Xt+∆≈Xt+f(Xt,t)∆+g(Xt,t) (Bt+∆−Bt). Usando a expressão acima para o incremento browniano, isso se torna Xt+∆≈Xt+f(Xt,t)∆+g(Xt,t)√∆εt. Essa equação descreve como o processo Xtevolui passo a passo: a cada intervalo ∆, há um deslocamento determinístico dado por f(Xt,t)∆(o drift) e um deslocamento aleatório g(Xt,t)√∆εt (a difusão). Essa formulação é conhecida como o esquema de Euler–Maruyama, uma generalização estocástica do método de Euler para equações diferenciais ordinárias. À medida que ∆→0, a sequência Xt+∆converge, sob condições adequadas, para a solução contínua da SDE. Exercício 28. Considere o esquema de Euler–Maruyama para simular uma equação diferencial estocástica (SDE) da forma dXt=f(Xt)dt +σdBt, onde Bté um movimento browniano padrão e σ>0controla a intensidade do ruído. (a) Implemente a simulação de trajetórias para dois campos de drift diferentes: f1(x) = −x e f2(x) = x−x3. (b) Gere várias trajetórias para cada caso, mantendo o mesmo valor de σe do passo temporal ∆. (c) Compare os comportamentos obtidos: –Como o termo de drift influencia a dispersão das trajetórias? –Em que sentido o caso f2(x) = x−x3pode ser interpretado como um sistema com dois estados de equilíbrio? (d) Plote em um mesmo gráfico uma trajetória com drift e outra sem drift (isto é, f(x) = 0) para visualizar a diferença qualitativa entre difusão pura e dinâmica com força restauradora. 128 CAPÍTULO 8. PROCESSOS DE DIFUSÃO de modo que os níveis de ruído cubram uniformemente várias ordens de magnitude, indo de ruído forte (σmax) a ruído fraco (σmin). Para cada amostra xie para cada σk, podemos gerar várias versões corrompidas ˜ xi,j,k=xi+εi,j,k,εi,j,k∼ N(0, σ2 kI), e calcular os respectivos alvos de regressão ti,j,k=−εi,j,k σ2 k . Dessa forma, podemos expandir o conjunto de dados criando várias amostras supervisionadas (˜ xi,j,k,σk,ti,j,k), que descrevem, para diferentes níveis de ruído, a direção de denoising a ser aprendida. O próximo passo é ajustar um modelo de regressão sθque receba como entrada o ponto corrompido ˜ xe o valor de σ, e aprenda a prever o vetor t. Esse modelo pode ser uma rede neural, mas também algo mais simples, como uma árvore de decisão ou um modelo de regressão não linear. O objetivo é que, após o treinamento, o campo aprendido satisfaça aproximadamente sθ(˜ x,σ)≈ ∇˜ xlog qσ(˜ x), fornecendo uma boa estimativa do score da densidade suavizada. Com isso, temos um modelo capaz de indicar, para cada ponto corrompido, em que direção ele deve se mover para recuperar regiões de alta densidade de p(x). Em resumo, o treinamento do DSM pode ser realizado seguindo os seguintes passos: 1. Definir uma escala geométrica de valores de ruído σ1, . . . , σK, que vai do ruído mais forte (σmax) ao mais fraco (σmin). Essa escala define o quanto cada amostra será corrompida. 2. Para cada amostra xido conjunto de dados e para cada nível de ruído σk, gerar algumas versões corrompidas ˜ xi,j,k=xi+εi,j,k, com εi,j,k∼ N(0, σ2 kI). Isso aumenta o tamanho do conjunto de treinamento e ajuda o modelo a aprender a remover diferentes intensidades de ruído. 8.4. DENOISING SCORE MATCHING 129 3. Calcular o alvo de regressão para cada par (˜ xi,j,k,σk)como ti,j,k=−εi,j,k/σ2 k, que representa a direção na qual a amostra corrompida deve ser movida para retornar à distribuição original. 4. Montar um conjunto de dados supervisionado formado por pares de entrada e saída ([ ˜ xi,j,k,σk],ti,j,k). O valor de σké incluído como uma feature adicional para indicar o nível de ruído daquela amostra. 5. Ajustar um modelo de regressão — que pode ser simples, como uma Floresta Aleatória — usando essas amostras expandidas. O modelo deve aprender a prever ta partir de [˜ x,σ]. 6. Após o treinamento, o modelo resultante sθ(˜ x,σ)fornece uma aproximação do campo de score ∇˜ xlog qσ(˜ x), indicando para cada ponto corrompido em qual direção ele deve se mover para se aproximar de regiões de alta densidade de p(x). Exercício 32. Neste exercício, vamos implementar o treinamento de um modelo de Denoising Score Matching (DSM) em um conjunto de dados sintético. O objetivo é aprender o campo de score sθ(˜ x,σ)a partir de amostras corrompidas, conforme discutido em aula. 1. Gere um conjunto de dados bidimensional usando a função make_moons abaixo. 2. Construa uma escala geométrica de valores de ruído usando a função geometric_sigmas. 3. Para cada valor de σ, corrompa as amostras adicionando ruído gaussiano ε∼ N(0, σ2I), e calcule o alvo t =−ε/σ2. 4. Monte um conjunto de dados supervisionado contendo como entrada o par [˜ x,σ]e como saída o vetor t. 5. Treine um modelo de regressão à sua escolha (por exemplo, regressão linear, rede neural, ou floresta aleatória) para aprender a mapear [˜ x,σ]7→ t. 6. Fixe um valor de σe visualize o campo aprendido sobre uma grade bidimensional de pontos, comparando visualmente os resultados de diferentes modelos. 8.4.2 Etapa de Inferência via Langevin Uma vez treinado o modelo de score sθ(˜ x,σ), podemos utilizá-lo para gerar novas amostras de uma distribuição aproximando p(x). A ideia é usar o campo aprendido como uma estimativa do gradiente do logaritmo da densidade, e então realizar uma simulação do processo de Langevin Anelado (Annealed Langevin Dynamics — ALD). O método segue a dinâmica estocástica xt=xt−1+αi 2sθ(xt−1,σi) + √αizt,zt∼ N(0, I), onde cada nível de ruído σicontrola a escala das atualizações e αié o passo de integração proporcional a σ2 i. Em termos práticos, seguimos a sequência de sigmas do maior (σmax) ao menor (σmin), de modo que as primeiras iterações façam o ponto explorar amplamente o espaço, e as últimas permitam um refinamento local. O algoritmo pode ser descrito assim: 130 CAPÍTULO 8. PROCESSOS DE DIFUSÃO 1. Inicialização: Inicie x0como uma amostra de uma distribuição de ruído, por exemplo x0∼ N(0, σ2 maxI). 2. Iteração sobre os níveis de ruído: Para cada σida sequência geométrica: • Calcule o passo de integração αi=ε·(σ2 i/σ2 min), onde ε>0 é um parâmetro fixo de escala. • Repita Tvezes (por exemplo, T=10): x←x+αi 2sθ(x,σi) + √αiz,z∼ N(0, I). 3. Saída: Após percorrer todos os níveis de ruído, o vetor final xTé uma amostra aproximada da distribuição de interesse p(x). Figura 8.2: Exemplo de amostragem via Annealed Langevin Dynamics. As trajetórias começam em ruído grande e gradualmente convergem para as regiões de alta densidade de p(x). A intuição é que, nas primeiras escalas de ruído, o modelo aprende apenas a estrutura global de p(x)— as regiões de alta densidade —, e conforme σdiminui, o processo de Langevin refina as amostras nessas regiões, capturando detalhes finos da distribuição. Exercício 33. Neste exercício, você deverá implementar o processo de Annealed Langevin Dynamics (ALD) para gerar novas amostras de uma distribuição aproximando o conjunto de dados moons. 8.4. DENOISING SCORE MATCHING 131 1. Gere o conjunto de dados bidimensional usando a função make_moons da biblioteca sklearn.datasets. 2. Treine um modelo de Denoising Score Matching (DSM) utilizando uma escala geométrica de valores de ruído σ1>σ2>··· >σK, conforme descrito na seção anterior. 3. Implemente o algoritmo de Annealed Langevin Dynamics, usando o campo aprendido sθ(x,σ) para atualizar as amostras segundo xt=xt−1+αi 2sθ(xt−1,σi) + √αizt,zt∼ N(0, I), percorrendo os níveis de ruído do maior (σmax) ao menor (σmin). 4. Após a simulação, visualize lado a lado o conjunto de dados original e as amostras geradas pelo processo de Langevin. Compare visualmente se as amostras geradas reproduzem a estrutura característica do moons. Exercício 34. Neste exercício, você deverá pensar em como adaptar todo o processo de Denoising Score Matching (DSM) e a etapa de inferência via Annealed Langevin Dynamics (ALD) para o caso condicional, em que desejamos modelar a distribuição p(y|x). 1. Relembre que, no caso não condicional, o modelo sθ(˜ x,σ)é treinado para aproximar o score da densidade suavizada ∇˜ xlog qσ(˜ x), onde qσé obtida pela convolução de p(x)com ruído gaussiano. Pense em como essa ideia pode ser estendida para o caso condicional, em que queremos o score ∇ylog qσ(y|x). 2. Escreva como ficaria o conjunto de treinamento supervisionado para o modelo condicional. Dica: ao corromper as variáveis de saída yicom ruído gaussiano εi,j,k∼ N(0, σ2 kI), o alvo passa a ser ti,j,k=−εi,j,k/σ2 k, e o modelo deve receber [˜ yi,j,k,xi,σk]como entrada. 3. Implemente o treinamento de um modelo de score sθ(y,x,σ)que aprenda o campo condicional de denoising. Você pode usar um modelo simples, como uma Floresta Aleatória ou uma rede neural. 4. Adapte o algoritmo de Annealed Langevin Dynamics para o caso condicional, mantendo x fixo e atualizando apenas y: yt=yt−1+αi 2sθ(yt−1,x,σi) + √αizt,zt∼ N(0, I). 132 CAPÍTULO 8. PROCESSOS DE DIFUSÃO 5. Escolha um conjunto de dados simples para testar o modelo, por exemplo: •Gere pares (x,y)com y =sin(2x) + |x|ε,ε∼ N(0,1). •Treine o modelo de score condicional. •Use o ALD condicional para gerar novas amostras de y para valores fixos de x. 6. Visualize os resultados mostrando, para alguns valores fixos de x, as distribuições das amostras geradas de y |x, comparando-as com os valores verdadeiros observados. Capítulo 9 Bootstrap 9.1 Uma visão pragmática de Bootstrap O nome bootstrap tem origem indireta nas histórias fantásticas do Barão de Münchhausen, personagem do século XVIII conhecido por narrar feitos impossíveis. Em uma de suas aventuras, o Barão conta ter conseguido sair de um pântano puxando a si mesmo pelo cabelo (junto com o cavalo), uma façanha evidentemente absurda. A expressão inglesa posterior “to pull oneself up by one’s bootstraps” — erguer-se puxando as próprias botas — tornou-se uma metáfora para realizar algo sem ajuda externa, e foi essa a imagem que inspirou Efron (1979) ao nomear seu método: um procedimento que, metaforicamente, se ergue sozinho. A ideia central do método de bootstrap é substituir a incerteza sobre a distribuição populacional Fpela incerteza induzida pela distribuição empírica ˆ Fn. Formalmente, seja X1, . . . , Xn∼F uma amostra i.i.d. e seja o parâmetro de interesse θ=t(F), para algum funcional tdefinido em um espaço apropriado de distribuições de probabilidade. O estimador empírico é ˆ θ=t(ˆ Fn), com ˆ Fn(x) = 1 n n ∑ i=1 1{Xi≤x}. Obootstrap consiste em gerar amostras X∗ 1, . . . , X∗ ni.i.d. de ˆ Fn(isto é, reamostrar com reposição dos dados observados), e então computar ˆ θ∗=t(ˆ F∗ n). A distribuição condicional de ˆ θ∗dado os dados é usada como aproximação para a distribuição amostral de ˆ θ. Exemplo 35 (Média amostral).Considere uma amostra X1, . . . , Xn∼F e o estimador usual da média ˆ θ=¯ X=1 n n ∑ i=1 Xi. O objetivo é quantificar a incerteza de ˆ θ, isto é, como ela variaria se repetíssemos o experimento várias vezes. Em termos formais, queremos aproximar a distribuição amostral de ˆ θ, P√n(ˆ θ−θ)≤a, 133 134 CAPÍTULO 9. BOOTSTRAP onde θ=E[X]. Em situações simples, podemos obter essa distribuição de forma analítica: se F for normal com variância σ2, então √n(¯ X−θ)∼ N(0, σ2). Entretanto, o bootstrap permite estimar a variabilidade de ¯ Xsem supor nada sobre F. A ideia é construir uma amostra artificial que imite o que aconteceria se o experimento fosse repetido. O procedimento é o seguinte: (1) A partir da amostra observada X1, . . . , Xn, sorteie com reposição n observações X∗ 1, . . . , X∗ n. Cada amostra reamostrada define uma distribuição empírica ˆ F∗ n. (2) Calcule a média de cada amostra reamostrada: ¯ X∗=1 n n ∑ i=1 X∗ i. (3) Repita o processo B vezes (por exemplo, B =1000), obtendo ¯ X∗(1), . . . , ¯ X∗(B). O conjunto dessas médias forma uma aproximação empírica da distribuição de ¯ X. O desvio padrão das médias reamostradas, ˆ σboot =v u u t1 B−1 B ∑ b=1¯ X∗(b)−¯ ¯ X∗2,¯ ¯ X∗=1 B B ∑ b=1 ¯ X∗(b), é uma estimativa do erro padrão de ¯ X. Como consequência, é possível construir intervalos de confiança para θusando os quantis das médias bootstrap: Iboot 1−α=¯ X∗ (α/2),¯ X∗ (1−α/2), onde os termos entre parênteses denotam os quantis empíricos da distribuição das médias reamostradas. A Figura 35 compara a distribuição verdadeira da média amostral (obtida por simulação Monte Carlo) com a distribuição condicional gerada pelo bootstrap, para F =N(5, 22)e n =30. 9.1. UMA VISÃO PRAGMÁTICA DE BOOTSTRAP 135 Exemplo 36 (Mediana).Enquanto a média é um estimador linear e de fácil análise, a mediana apresenta um comportamento mais sutil. Se X1, . . . , Xn∼F e ˆ θ=ˆ F−1 n(0.5), sua variabilidade depende da densidade de F no ponto da mediana θ, pois pequenas flutuações em F se traduzem em variações maiores ou menores na posição onde F(x) = 1/2. Uma forma breve de derivar a variância assintótica é a seguinte. Como F(θ) = 1/2 eˆ Fn(ˆ θ) = 1/2, podemos relacionar ˆ θeθpor uma expansão local de F em torno de θ: F(ˆ θ)≈F(θ) + f(θ)( ˆ θ−θ), onde f(θ) = F′(θ). Como ˆ Fnconverge uniformemente para F (Teorema de Glivenko–Cantelli), é lícito substituir F(·)por ˆ Fn(·)nessa aproximação sem alterar o termo assintótico dominante — a diferença é da ordem op(1/√n). Substituindo F(ˆ θ)por ˆ Fn(ˆ θ) = 1/2, obtemos: 0≈ˆ Fn(θ)−F(θ) + f(θ)( ˆ θ−θ). Multiplicando por √n, √n(ˆ θ−θ)≈ − 1 f(θ)√nˆ Fn(θ)−F(θ). Pelo Teorema Central do Limite empírico, √nˆ Fn(θ)−F(θ)⇒ N0, F(θ)(1−F(θ))=N(0,1/4). Portanto, Var ˆ θ≈1 4n f(θ)2. O bootstrap oferece uma alternativa direta. Partindo da amostra observada, reamostra-se com reposição B vezes e calcula-se a mediana em cada reamostra, ˆ θ∗(b)=median(X∗(b) 1, . . . , X∗(b) n),b=1, . . . , B. 136 CAPÍTULO 9. BOOTSTRAP A variabilidade entre as medianas reamostradas fornece uma estimativa do erro padrão de ˆ θ, ˆ σboot =v u u t1 B−1 B ∑ b=1ˆ θ∗(b)−¯ θ∗2,¯ θ∗=1 B B ∑ b=1 ˆ θ∗(b). Dessa forma, mesmo sem conhecer f(θ)nem a forma de F, é possível avaliar empiricamente a incerteza da mediana. Exercício 35. Reproduza os experimentos dos exercícios anteriores. 9.2 Uma visão teórica de Bootstrap À primeira vista, o funcionamento do bootstrap pode parecer misterioso. Afinal, estamos tentando aproximar a variabilidade de um estimador — algo que depende da distribuição populacional desconhecida F— usando apenas a distribuição empírica ˆ Fn, construída a partir de uma única amostra. Em outras palavras, substituímos o próprio objeto que queremos inferir por uma aproximação baseada nos dados: o método “ergue-se” sobre si mesmo, exatamente como sugere seu nome. O fato de essa substituição produzir resultados válidos não é trivial. Não há, a princípio, nenhuma razão óbvia para que as flutuações de um estimador calculado a partir de ˆ Fnreflitam corretamente aquelas que surgiriam se repetíssemos o experimento sob F. Ainda assim, sob condições gerais, o bootstrap funciona — e de forma surpreendentemente robusta. 9.2.1 A desigualdade de Dvoretzky–Kiefer–Wolfowitz Antes de entender por que o bootstrap funciona, é preciso quantificar o quão bem a distribuição empírica ˆ Fnaproxima a verdadeira F. A desigualdade de Dvoretzky–Kiefer–Wolfowitz (DKW) é o ponto de partida: ela fornece uma garantia não assintótica, válida para qualquer n, de que as duas funções de distribuição estão próximas com alta probabilidade. Mais precisamente, para amostras i.i.d. X1, . . . , Xn∼F, vale que Psup xˆ Fn(x)−F(x)>ε≤Ce−nε2,∀ε>0. Essa desigualdade mostra que o erro uniforme entre Feˆ Fndecai exponencialmente com n. Em particular, ˆ Fnconverge quase certamente para F, o que é o conteúdo do teorema de Glivenko–Cantelli. Esse resultado é notável por duas razões. Primeiro, ele é completamente não assintótico: a probabilidade de desvio pode ser controlada explicitamente para qualquer n. Segundo, ele já sugere a ideia central do bootstrap — se ˆ Fnestá uniformemente próxima de F, então estimadores baseados em uma ou outra devem ter comportamentos muito semelhantes. O primeiro ingrediente para a prova da DKW é a desigualdade das diferenças finitas, também conhecida como desigualdade de McDiarmid. Ela fornece um limite de concentração para funções de variáveis independentes cujo valor não muda muito quando uma única observação é alterada. 9.2. UMA VISÃO TEÓRICA DE BOOTSTRAP 137 Teorema 12 (Desigualdade de McDiarmid).Sejam X1, . . . , Xnvariáveis independentes assumindo valores em um conjunto arbitrário X, e seja f :Xn→Ruma função tal que |f(x1, . . . , xi, . . . , xn)−f(x1, . . . , x′ i, . . . , xn)| ≤ ci para todo i e para todos os valores possíveis das variáveis. Então, para todo ε>0, P(f(X1, . . . , Xn)−E[f(X1, . . . , Xn)]≥ε)≤exp−2ε2 ∑n i=1c2 i. A ideia é simples: se cada variável individual tem influência limitada sobre o valor final de f, então f(X1, . . . , Xn)não pode se desviar muito de sua média. Essa desigualdade é uma generalização do lema de Hoeffding para funções simétricas e não lineares das observações. No caso da DKW, aplicamos esse resultado à função f(X1, . . . , Xn) = sup x|ˆ Fn(x)−F(x)|. Observe que alterar um único Ximuda no máximo um termo da soma que define ˆ Fn(x), e portanto o valor de fsó pode variar em 1/n. Assim, podemos tomar ci=1/npara todo i, o que dá n ∑ i=1 c2 i=n×1 n2=1 n. Substituindo isso na desigualdade de McDiarmid, obtemos Psup x|ˆ Fn(x)−F(x)|−Esup x|ˆ Fn(x)−F(x)|≥ε≤e−2nε2. Esse é o passo essencial da prova da DKW: ele mostra que a distância uniforme entre ˆ FneF está fortemente concentrada em torno de sua média. O passo seguinte consiste em controlar o valor esperado Esupx|ˆ Fn(x)−F(x)|. Para controlar o valor esperado de supx|ˆ Fn(x)−F(x)|, precisamos entender primeiro um caso mais simples: o comportamento do valor esperado do máximo de um número finito de variáveis aleatórias. 144 CAPÍTULO 9. BOOTSTRAP Capítulo 10 Estratégias para acelerar códigos em Python 10.1 Profiling com cProfile Antes de otimizar, é essencial medir onde o tempo realmente está sendo gasto. O módulo cProfile, da biblioteca padrão do Python, permite gerar um perfil de execução mostrando quantas vezes cada função foi chamada e quanto tempo ela consumiu. O uso mais simples é direto pelo terminal, aplicando o profiler a um script: 1# Executa o script inteiro e mostra estatisticas 2python -m cProfile meu_script . py 3 4# Salva os resultados em arquivo para analise posterior 5python -m cProfile -o saida . prof meu_script . py Rodando pelo terminal O arquivo gerado pode ser inspecionado com o módulo pstats, que permite ordenar e filtrar resultados: 1python -m pstats saida . prof 2# Comandos uteis no prompt do pstats : 3# sort time ( ordena pelo tempo interno da funcao ) 4# sort cumtime ( ordena pelo tempo acumulado ) 5# stats 20 ( mostra as 20 funcoes mais custosas ) 6# callers func ( quem chama ’func ’) 7# callees func ( quem ’ func ’ chama ) Explorando com pstats (terminal) Também é possível usar cProfile dentro do código, o que facilita em notebooks ou quando queremos medir apenas um trecho específico: 1import cProfile , pstats , io 2 3pr = cProfile . Profile () 145 146 CAPÍTULO 10. ESTRATÉGIAS PARA ACELERAR CÓDIGOS EM PYTHON 4pr. enable () 5 6# --- codigo a ser medido --- 7resultado = algoritmo_pesado () 8# ---------------------------- 9 10 pr. disable () 11 s = io. StringIO () 12 ps = pstats . Stats (pr , stream =s). sort_stats ("cumtime") 13 ps. print_stats (10) # mostra as 10 funcoes mais custosas 14 print (s. getvalue ()) Usando cProfile dentro do código As duas métricas principais são: •time: tempo gasto apenas dentro da função, sem contar chamadas internas. •cumtime: tempo acumulado, incluindo todas as funções chamadas. Em geral, começa-se ordenando por cumtime para encontrar o caminho mais caro da execução. Depois, olhar o time ajuda a identificar funções individuais que valem otimização. A seguir montamos um experimento simples para evidenciar como o cProfile ajuda a localizar gargalos: comparamos uma multiplicação de matrizes feita de forma ingênua em Python (três laços) com a versão vetorizada do NumPy (delegada à BLAS). O código abaixo implementa as duas versões e usa uma função auxiliar para rodar o profiler em cada uma delas, exibindo as funções mais custosas. O script pode ser salvo como profile_matmul.py. 1import numpy as np 2import math 3import cProfile , pstats , io 4import time 5 6# Versao ingenua : 3 loops em Python 7def matmul_naive(A, B): 8n, m = A. shape 9m2 , p = B. shape 10 assert m == m2 11 C = np.zeros ((n, p)) 12 for iin range (n): 13 for jin range (p): 14 s = 0.0 15 for kin range (m): 16 s += A[i, k] * B[k, j] 17 C[i, j] = s 18 return C 19 10.1. PROFILING COM CPROFILE 147 20 # Versao NumPy ( vetorizada / BLAS ) 21 def matmul_numpy(A, B): 22 return A@B 23 24 def profile_func ( func , *args , top =15) : 25 pr = cProfile . Profile () 26 pr. enable () 27 t0 = time . perf_counter () 28 result = func (* args ) 29 t1 = time . perf_counter () 30 pr. disable () 31 s = io. StringIO () 32 ps = pstats . Stats (pr , stream =s). sort_stats ("cumtime") 33 ps. print_stats ( top ) 34 print (f"\n>>>> Tempo total ( parede ): {t1 - t0 :.3f} s\n") 35 return s. getvalue () 36 37 A = np. random . rand (n, n) 38 B = np. random . rand (n, n) 39 40 print (" ==== Profiling matmul_naive ==== ") 41 out_naive = profile_func ( matmul_naive , A , B) 42 print ( out_naive ) 43 44 print (" ==== Profiling matmul_numpy ==== ") 45 out_np = profile_func ( matmul_numpy , A , B) 46 print (out_np) Esse script pode ser executado normalmente com python profile_matmul.py. Outra forma é rodar o profiler diretamente no terminal, usando python -m cProfile -o saida.prof profile_matmul.py. Nesse caso o resultado fica salvo em saida.prof, e podemos explorá-lo depois com o módulo pstats de forma interativa, usando comandos como sort cumtime,stats 20 ou callers matmul_naive. Rodando a versão ingênua, a saída típica mostra que praticamente todo o tempo foi consumido dentro de matmul_naive: 17 function calls in 8.532 seconds 2 3Ordered by: cumulative time 4ncalls tottime percall cumtime percall filename : lineno ( function ) 51 8.523 8.523 8.523 8.523 profile_matmul .py :11( matmul_naive) 61 0.009 0.009 0.009 0.009 {built -in method builtins . print } 7... Ao comparar com a versão vetorizada, vemos que a execução termina em milésimos de segundo, com o tempo todo acumulado em matmul_numpy: 148 CAPÍTULO 10. ESTRATÉGIAS PARA ACELERAR CÓDIGOS EM PYTHON 17 function calls in 0.020 seconds 2 3Ordered by: cumulative time 4ncalls tottime percall cumtime percall filename : lineno ( function ) 51 0.019 0.019 0.019 0.019 profile_matmul .py :28( matmul_numpy) 6... Os números exatos variam conforme o tamanho das matrizes e a biblioteca BLAS instalada, mas o padrão é claro: a implementação ingênua em Python puro consome segundos de CPU, enquanto a versão NumPy é milhares de vezes mais rápida. As colunas do profiler têm significados diferentes. O campo ncalls mostra o número de chamadas à função. O tottime corresponde ao tempo gasto apenas dentro da função, sem contar chamadas internas. Já o cumtime indica o tempo acumulado incluindo funções chamadas dentro dela. Em geral, ordenar por cumtime ajuda a encontrar o caminho mais custoso da execução, enquanto olhar para tottime revela funções “folha” particularmente lentas. Quando esse mesmo código é rodado em um notebook Jupyter, o output tende a ficar mais “poluído”, aparecendo referências a asyncio,zmq e outros componentes do kernel. Isso acontece porque o profiler mede tudo o que roda no processo, não apenas a nossa função. Para uma visão limpa e didática, vale a pena executar o script direto no terminal. 10.2 Paralelização com joblib.Parallel A biblioteca joblib fornece uma forma simples de paralelizar loops embaraçosamente paralelos em Python, isto é, situações em que várias tarefas independentes podem ser executadas ao mesmo tempo. A ideia básica é escrever um laço for como uma compreensão preguiçosa de chamadas a uma função via delayed, e despachar essas tarefas para Parallel, que se encarrega de distribuílas entre diferentes trabalhadores. 1from joblib import Parallel , delayed 2from math import sqrt 3 4# aplicar sqrt a 0^2 , 1^2 , ... , 9^2 em paralelo 5res = Parallel ( n_jobs =4) ( 6delayed ( sqrt )(i **2) 7for iin range (10) 8) Receita de bolo No exemplo acima, o parâmetro n_jobs define quantos trabalhadores serão usados (tipicamente o número de CPUs lógicas da máquina). A função delayed apenas empacota a chamada para que ela possa ser enviada a um worker, enquanto Parallel recolhe todas as tarefas e coordena sua execução. Uma forma intuitiva de entender esse mecanismo é pensar em uma cozinha: se temos apenas um cozinheiro (um for sequencial), cada prato é preparado do início ao fim antes do próximo 10.2. PARALELIZAÇÃO COM JOBLIB.PARALLEL 149 começar. Já com vários cozinheiros (workers), cada um recebe um prato e trabalha nele independentemente, de modo que vários ficam prontos ao mesmo tempo. Essa estratégia funciona muito bem, mas há alguns cuidados: se uma tarefa demora muito enquanto outras são rápidas, pode haver desequilíbrio entre os workers; por outro lado, se existem milhares de tarefas minúsculas, o custo de despachá-las pode ser maior que o ganho da paralelização. Para reduzir esse problema, ojoblib agrupa chamadas em lotes (batching), enviando várias de uma vez só. Outro detalhe importante está no backend usado. Em Python, o Global Interpreter Lock (GIL) impede que várias threads executem código Python puro ao mesmo tempo. Por isso, o backend padrão (loky) cria processos separados, que contornam o GIL e escalam bem em cálculos pesados. Já o backend threading mantém as tarefas no mesmo processo, sendo útil em funções que passam a maior parte do tempo esperando I/O ou que já liberam o GIL (como operações NumPy). Existe ainda o multiprocessing, mas o loky tende a ser mais robusto. 1# Uso de threads porque a funcao processa_io 2# passa a maior parte do tempo esperando rede. 3res = Parallel ( n_jobs =8, backend =" threading ")( 4delayed ( processa_io )(u) for uin urls 5) Exemplo com I/O Em resumo: use loky (padrão) para tarefas CPU-bound, threading para tarefas I/O-bound, e sempre ajuste o número de jobs de acordo com o hardware disponível. Paralelizar acelera muito, mas nem sempre compensa: quando as tarefas são pequenas demais, o overhead pode superar o benefício. Um exemplo clássico de tarefa CPU-bound é calcular números primos ou executar operações pesadas de álgebra linear. Nesses casos, vale usar o backend padrão: 1from joblib import Parallel , delayed 2import math 3 4def eh_primo (n): 5for iin range (2, int ( math. sqrt (n)) +1) : 6if n % i == 0: 7return False 8return True 9 10 nums = range (10**6 , 10**6+1000) 11 res = Parallel ( n_jobs =4) ( delayed ( eh_primo )(n) for nin nums) Exemplo CPU-bound Aqui, cada worker testa um conjunto de números independentemente. Quanto mais núcleos disponíveis, mais rápido o processamento. Já um exemplo I/O-bound seria baixar várias páginas da web. Cada tarefa fica a maior parte do tempo esperando a rede, e usar processos separados não traz vantagem; nesse caso o backend threading é mais leve: 150 CAPÍTULO 10. ESTRATÉGIAS PARA ACELERAR CÓDIGOS EM PYTHON 1import requests 2urls = [" https :// httpbin . org / delay /1 "] * 20 3 4def baixa ( url ): 5return requests . get ( url ). status_code 6 7res = Parallel ( n_jobs =8, backend =" threading ")( 8delayed ( baixa )(u) for uin urls 9) Exemplo I/O-bound Se cada requisição demora cerca de 1 segundo, com 8 threads as 20 requisições terminam em poucos segundos, em vez de mais de 20. Por fim, um caso em que a paralelização atrapalha é quando as tarefas são rápidas demais, por exemplo calcular o quadrado de números pequenos: 1def quadrado (n): 2return n*n 3 4nums = range (1000) 5res = Parallel ( n_jobs =4) ( delayed ( quadrado )(n) for nin nums) Exemplo de overhead Aqui o custo de organizar as tarefas, mandar para os workers e reunir os resultados é maior do que simplesmente rodar um for sequencial. Nesse cenário, a paralelização pode ser mais lenta. 10.3 Compilação Just-In-Time com Numba Numba é um compilador JIT (Just-In-Time) para Python focado em acelerar código numérico. Ele “traduz” funções Python (que operam sobre tipos e arrays compatíveis) para código nativo via LLVM, reduzindo drasticamente o overhead dos laços em Python puro. A ideia prática é simples: decorar funções críticas com @njit (ou @jit(nopython=True)), evitar objetos Python dentro dessas funções e, quando fizer sentido, ativar paralelização com parallel=True eprange. O primeiro cuidado ao medir é lembrar do custo de compilação: na primeira chamada de cada assinatura de tipos, Numba compila a função (demora mais). Depois disso, as chamadas seguintes usam o código nativo já gerado. O exemplo abaixo acelera uma multiplicação de matrizes ingênua (três laços) sem recorrer ao NumPy @. Primeiro mostramos a versão njit sequencial; em seguida, a variação paralela (parallel=True +prange). Usamos perf_counter para mostrar o tempo da primeira chamada (com compilação) e das chamadas seguintes (sem compilação). 1import numpy as np 2import time 3from numba import njit , prange 10.3. COMPILAÇÃO JUST-IN-TIME COM NUMBA 151 4 5# Versao Python pura ( referencia ) 6def matmul_naive(A, B): 7n, m = A. shape 8m2 , p = B. shape 9assert m == m2 10 C = np.zeros ((n, p)) 11 for iin range (n): 12 for jin range (p): 13 s = 0.0 14 for kin range (m): 15 s += A[i, k] * B[k, j] 16 C[i, j] = s 17 return C 18 19 # Versao Numba : nopython mode (sem objetos Python dentro ) 20 @njit 21 def matmul_numba(A, B): 22 n, m = A. shape 23 m2 , p = B. shape 24 C = np.zeros ((n, p)) 25 for iin range (n): 26 for jin range (p): 27 s = 0.0 28 for kin range (m): 29 s += A[i, k] * B[k, j] 30 C[i, j] = s 31 return C 32 33 # Versao Numba paralela : requer parallel = True e uso de prange 34 @njit ( parallel = True ) 35 def matmul_numba_parallel(A, B): 36 n, m = A. shape 37 m2 , p = B. shape 38 C = np.zeros ((n, p)) 39 for iin prange (n): # <-- prange permite paralelizar esse loop externo 40 for jin range (p): 41 s = 0.0 42 for kin range (m): 43 s += A[i, k] * B[k, j] 44 C[i, j] = s 45 return C 46 47 # Benchmark simples : separa " primeira chamada " e " repetidas " 48 def bench ( func , * args , repeat =3, label =""): 49 # primeira chamada ( inclui compilacao JIT quando aplicavel ) 50 t0 = time . perf_counter () 152 CAPÍTULO 10. ESTRATÉGIAS PARA ACELERAR CÓDIGOS EM PYTHON 51 out = func (* args ) 52 t1 = time . perf_counter () 53 print (f"{ label } [1a chamada ]: { t1 - t0 :.3 f} s") 54 55 # chamadas seguintes (ja compilado ) 56 best = float (" inf") 57 for _in range (repeat): 58 t0 = time . perf_counter () 59 func(*args) 60 t1 = time . perf_counter () 61 best = min(best , t1 - t0) 62 print (f"{ label } [ melhor chamada subsequente ]: { best :.3 f} s") 63 return out 64 65 if __name__ == " __main__ ": 66 n = 600 67 A = np. random . rand (n, n) 68 B = np. random . rand (n, n) 69 70 # Referencia Python puro ( lento ) 71 bench ( matmul_naive , A , B , label =" naive ( Python )") 72 73 # Numba sequencial 74 bench ( matmul_numba , A , B , label =" Numba ( njit )") 75 76 # Numba paralelo 77 bench ( matmul_numba_parallel , A , B , label =" Numba ( parallel )") Acelerando loops com Numba (@njit) Na prática, você deverá observar algo assim: a versão Python pura leva segundos; a versão @njit cai para frações (ou poucos segundos em matrizes grandes) após a compilação; a versão paralela tende a ganhar mais em máquinas com vários núcleos, desde que o tamanho do problema justifique o overhead de criar e sincronizar threads. Nem todo laço se beneficia de parallel=True; se o problema é pequeno, o custo extra pode superar o ganho. Outro modo útil de Numba é compilar funções elementwise com @vectorize, criando uma ufunc ao estilo NumPy; isso permite aplicar a função diretamente sobre arrays, com broadcast, sem escrever laços em Python. O exemplo a seguir define uma ufunc para uma transformação escalar simples e a aplica a um array grande. 1import numpy as np 2from numba import vectorize , float64 3 4@vectorize ([ float64 ( float64 ) ]) 5def transform (x): 6# alguma transformacao escalar ( exemplo ) 7return (x * x + 0.5) / (x + 1.0) 8 10.4. PARALELISMO SIMPLES EM BASH 153 9x = np. random . rand (1 _000_000 ) 10 y = transform (x) # aplica como ufunc , sem lacos explicitos em Python UFunc com @vectorize (estilo NumPy) Algumas recomendações práticas ao usar Numba: (i) mantenha dentro das funções JIT apenas operações suportadas (aritmética, indexação NumPy, algumas funções math/numpy); (ii) evite objetos Python (listas que crescem, dicionários, set) e chamadas que exijam o interpretador; (iii) prefira arrays com dtype numéricos (float64, int64, etc.) e formatos contíguos; (iv) tome cuidado com alocação excessiva dentro do laço; (v) ative parallel=True apenas após confirmar que o gargalo é CPU-bound e que o tamanho do problema compensa a paralelização; (vi) lembre-se do “aquecimento”: meça separando a primeira chamada (com compilação) das seguintes; (vii) quando a função estabilizar, @njit(cache=True) pode salvar o binário no disco e reduzir o tempo de compilação em execuções futuras (útil em scripts). Por fim, se você já tem uma versão vetorizada eficiente em NumPy (que usa BLAS), muitas vezes ela será tão rápida quanto (ou mais rápida que) reimplementar em Numba, a menos que o seu padrão de acesso/cálculo seja muito específico. O ponto forte do Numba é acelerar laços e lógicas numéricas que seriam lentas em Python puro, mantendo o código próximo ao original, sem partir direto para C/C++. 10.4 Paralelismo simples em Bash O Bash permite escrever pequenos scripts para automatizar tarefas repetitivas. Um dos recursos mais uteis é a possibilidade de rodar varios comandos em paralelo, sem esperar um terminar para comecar o proximo. Para isso usamos o operador &. No exemplo abaixo, usamos o comando sleep (que apenas dorme por alguns segundos) para simular tarefas demoradas. Cada chamada ao sleep é enviada ao plano de fundo com &, de modo que o laco continua imediatamente para a proxima iteracao. 1#!/ bin /bash 2 3for iin $(seq 1 5) 4do 5echo " Iniciando tarefa $i" 6sleep 3 & 7done 8 9echo " Todas as tarefas foram lancadas !" Rodando sleeps em paralelo Nesse script, as cinco tarefas comecam quase ao mesmo tempo e, apos cerca de tres segundos, todas terminam juntas. Se tirassemos o &, o script levaria cerca de 15 segundos, pois cada sleep 3seria executado em sequencia. Para visualizar essa diferenca, vejamos primeiro a execucao sequencial: 1#!/ bin /bash