Metodos Computacionais
Metodos Computacionais
10 de novembro de 2025
2
Sumário
1 Introdução 7
1.1 Um conselho: a importância de ser ruim antes de ser bom . . . . . . . . . . . . . . . 7
3
4 SUMÁRIO
3.7 Distribuição hipergeométrica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.7.1 Simulando a Hipergeométrica . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
6 Redução de variância 69
6.1 Uso de variáveis antitéticas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
6.2 O uso de variáveis de controle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
6.3 Redução de Variância por Condicionamento . . . . . . . . . . . . . . . . . . . . . . . 77
10 Bootstrap 133
10.1 Uma visão pragmática de Bootstrap . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
10.2 Uma visão teórica de Bootstrap . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
10.2.1 A desigualdade de Dvoretzky–Kiefer–Wolfowitz . . . . . . . . . . . . . . . . 136
10.2.2 Bootstrap e pontes brownianas . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
Introdução
7
8 CAPÍTULO 1. INTRODUÇÃO
Capítulo 2
n
P [ A1 ∪ · · · ∪ A n ] = ∑ P[ A i ].
i =1
P[ A ∩ B ]
P[ A | B ] = .
P[ B ]
Dois eventos A e B sã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
10 CAPÍTULO 2. ELEMENTOS BÁSICOS DE PROBABILIDADE
De forma equivalente, a independência entre A e B pode ser expressa afirmando que P[ A | B] =
P[ A], sempre que P[ B] ̸= 0.
Além disso, uma sequência de variáveis aleatórias é dita i.i.d. (independentes e identicamente
distribuídas) quando todas as variáveis da sequência são mutuamente independentes e seguem
a mesma distribuição de probabilidade.
Seguem algumas propriedades importantes:
P[ A ∪ B ] = P[ A ] + P[ B ] − P[ A ∩ B ] (regra da soma)
" #
n n
P ∑ P[ A i ]
[
Ai ≤ (desigualdade da união)
i =1 i =1
P[ B | A ]P[ A ]
P[ A | B ] = (fórmula de Bayes)
P[ B ]
" # " #
n −1
n\
P A i = P[ A1 ]P[ A2 | A1 ] · · · P A n |
\
Ai (regra da cadeia)
i =1 i =1
x 7 → P[ X = x ].
Em muitos contextos, o interesse recai não apenas sobre variáveis aleatórias individuais, mas
também sobre o relacionamento entre duas ou mais variáveis. Para descrever a dependência
entre X e Y, definimos a função de distribuição acumulada conjunta como
F ( x, y) = P{ X ≤ x, Y ≤ y},
p( x, y) = P{ X = x, Y = y}.
Se forem conjuntamente contínuas, existe uma função densidade de probabilidade conjunta f ( x, y) tal
que, para quaisquer conjuntos C, D ⊂ R,
ZZ
P{ X ∈ C, Y ∈ D } = f ( x, y) dx dy.
x ∈C, y∈ D
P{ X ∈ C, Y ∈ D } = P{ X ∈ C } P{Y ∈ D }.
De forma intuitiva, isso significa que conhecer o valor de uma das variáveis não altera a distri-
buição da outra.
No caso discreto, X e Y são independentes se, e somente se, para todo x, y,
P{ X = x, Y = y} = P{ X = x } P{Y = y}.
f ( x, y) = f X ( x ) f Y (y), ∀ x, y,
então
E[ I ] = 1 · P( A ) + 0 · P( A c ) = P( A ).
Portanto, a esperança de uma variável indicadora de um evento A é exatamente a probabilidade de que A
ocorra.
12 CAPÍTULO 2. ELEMENTOS BÁSICOS DE PROBABILIDADE
No caso contínuo, quando X possui uma função densidade de probabilidade f ( x ), a espe-
rança é dada por Z ∞
E[ X ] = x f ( x ) dx.
−∞
Além disso, dado uma função qualquer g, temos que:
Z ∞
E[ g( X )] = g( x ) f ( x ) dx.
−∞
Uma propriedade fundamental da esperança é sua linearidade. Isto é, para quaisquer variá-
veis aleatórias X e Y e constantes a, b ∈ R, temos:
E[ aX + bY ] = aE[ X ] + bE[Y ].
2.4 Variância
A variância de uma variável aleatória X é denotada por Var[ X ] e definida como
Var[ aX ] = a2 Var[ X ].
2.4.1 Covariância
A covariância entre duas variáveis aleatórias X e Y é denotada por Cov( X, Y ) e definida por
Exercício 4. Seja X uniforme no intervalo [−1, 1] e seja Y = X 2 . Mostre que Cov ( X, Y ) = 0 mas X, Y
não são independentes.
2.4. VARIÂNCIA 13
Observação 1. Considere uma variável aleatória contínua X centrada em zero, ou seja, E[ X ] = 0, com
densidade de probabilidade par e definida em um intervalo do tipo (− a, a), com a > 0. Seja Y = g( X )
para uma função g. A questão é: para quais funções g( X ) temos Cov( X, g( X )) = 0?
Sabemos que
Cov( X, g( X )) = E[ Xg( X )] − E[ X ]E[ g( X )].
Como E[ X ] = 0, segue que Cov( X, g( X )) = E[ Xg( X )]. Denotando a densidade de X por f ( x ), temos
Z a
Cov( X, g( X )) = xg( x ) f ( x )dx.
−a
Uma maneira de garantir que Cov( X, g( X )) = 0 é exigir que g( x ) seja uma função par. Assim,
xg( x ) f ( x ) será uma função ímpar e a integral em (− a, a) se anulará, ou seja,
Z a
xg( x ) f ( x )dx = 0.
−a
A covariância é uma forma bilinear simétrica e semi-definida positiva, com as seguintes pro-
priedades:
Além disso, vale a desigualdade de Cauchy-Schwarz, que afirma que para variáveis X e Y
com variância finita, q
|Cov( X, Y )| ≤ Var[ X ] Var[Y ].
Portanto, C(X) é a matriz cujos elementos são Cov( Xi , X j ). Além disso, é imediato mostrar
que
C(X) = E[XX⊤ ] − E[X] E[X]⊤ .
14 CAPÍTULO 2. ELEMENTOS BÁSICOS DE PROBABILIDADE
2.5 Desigualdades básicas de concentração
Nesta seção, apresentamos duas desigualdades fundamentais que estabelecem limites superiores
para a probabilidade de uma variável aleatória assumir valores distantes de sua média. Tais
resultados são amplamente utilizados em probabilidade, estatística e teoria da informação para
analisar o comportamento de caudas de distribuições.
A primeira delas é a Desigualdade de Markov, que fornece um limite simples para variáveis
aleatórias não-negativas em função apenas de sua esperança.
E[ X ]
P( X ≥ t ) ≤ .
t
Exercício 6. Prove a desigualdade de Markov. Dica: use o fato de que x
t ≥ I { x ≥ t }.
Teorema 2 (Desigualdade de Chebyshev). Seja X uma variável aleatória com valor esperado µ = E[ X ]
e variância finita Var( X ) = σ2 . Então, para todo ε > 0, vale:
σ2
P(| X − µ| ≥ ε) ≤ .
ε2
Teorema 3 (Lei Fraca dos Grandes Números). Seja ( Xn )n∈N uma sequência de variáveis aleatórias
independentes, todas com a mesma esperança µ e variância σ2 < ∞. Definindo a média amostral por
1 n
n i∑
Xn = Xi ,
=1
Exercício 8. Prove a Lei Fraca dos Grandes números utilizando a desigualdade de Chebyshev.
2.6. TEOREMAS ASSINTÓTICOS 15
Teorema 4 (Teorema Central do Limite). Seja X1 , . . . , Xn uma sequência de variáveis aleatórias i.i.d.
com esperança µ, variância σ2 e momento de ordem 3 finito. Definimos a média amostral como
1 n
n i∑
Xn = Xi .
=1
Então, √
n( X n − µ) d
−→ N (0, 1).
σ
Demonstração. Suponha, sem perda de generalidade, que µ = 0 e σ = 1. Defina
1 n 1 n
A n = √ ∑ Xi e Bn = √ ∑ Ni ,
n i =1 n i =1
i.i.d.
onde Ni ∼ N (0, 1) independentes de tudo. Note que Bn ∼ N (0, 1) para todo n.
d
Para provar que An → N (0, 1), é suficiente mostrar que, para qualquer função de teste f
suave e com crescimento controlado,
E[ f ( An )] − E[ f ( Bn )] −→ 0.
(0) (n)
Claramente, Cn = An e Cn = Bn .
Assim,
(0) (n)
E[ f ( An )] − E[ f ( Bn )] = E[ f (Cn )] − E[ f (Cn )]
n
= ∑ ∆k ,
k =1
onde
( k −1) (k)
∆k := E[ f (Cn ) − f (Cn )].
(k) ( k −1)
Passo 2: Isolando o termo que difere. Entre Cn e Cn , o único termo diferente é o k-ésimo.
Definamos
(k)
√1
Dn = n
N1 + · · · + Nk−1 + 0 + Xk+1 + · · · + Xn ,
(k) ( k −1)
isto é, a parte comum entre Cn e Cn , mas com o k-ésimo termo anulado.
Assim,
(k) (k) Nk ( k −1) (k) Xk
Cn = Dn + √
n
, Cn = Dn + √
n
.
16 CAPÍTULO 2. ELEMENTOS BÁSICOS DE PROBABILIDADE
Portanto, h i
(k) Xk (k) Nk
∆ k = E f Dn + √
n
− f Dn + √
n
.
(k)
Passo 3: Expansão de Taylor condicional. Fixe Dn = d. Aplicando Taylor em torno de d, temos:
Xk2 ′′ Xk3
Xk Xk ′
f d+ √ n
= f ( d ) + √
n
f ( d ) + 2n f ( d ) + 6n3/2
f (3) ( d + ξ X ) ,
N2 Nk3
Nk
f d+ √ n
= f (d) + √Nnk f ′ (d) + 2nk f ′′ (d) + 6n3/2 f (3) ( d + ξ N ) ,
√ √
para alguns ξ X , ξ N entre 0 e Xk / n ou Nk / n.
Subtraindo,
Xk Nk
f d+ √ n
− f d + √
n
= √1n f ′ (d)( Xk − Nk ) + 1 ′′ 2
2n f ( d )( Xk − Nk2 ) + Rk (d),
onde
Rk (d) = 1
6n3/2
Xk3 f (3) (d + ξ X ) − Nk3 f (3) (d + ξ N ) .
Passo 4: Tomando esperança condicional. Voltamos para
h i
(k) Xk (k) Nk
∆ k = E f Dn + √ n
− f Dn + √
n
.
(k)
Usando a decomposição anterior e condicionando em Dn , temos:
h i
(k)
∆k = E √1n f ′ ( Dn )( Xk − Nk )
h i
1 ′′ (k)
+ E 2n f ( Dn )( Xk − Nk )
2 2
+ E[ R k ].
(k)
Agora, como Xk e Nk são independentes de Dn , obtemos:
(k) (k)
E[ f ′ ( Dn )( Xk − Nk )] = E[ f ′ ( Dn )] · (E[ Xk ] − E[ Nk ]) = 0,
(k) (k)
E[ f ′′ ( Dn )( Xk2 − Nk2 )] = E[ f ′′ ( Dn )] · (E[ Xk2 ] − E[ Nk2 ]) = 0.
Portanto, só resta
∆ k = E[ R k ].
Passo 5: Controle do resto. Do termo Rk , temos
1
| Rk | ≤ 3/2 | Xk |3 sup | f (3) | + | Nk |3 sup | f (3) | .
6n
Tomando esperança,
C
|E[ Rk ]| ≤ 3/2 E[| X1 |3 ] + E[| N1 |3 ] ,
n
1 ( 3 )
onde C = 6 sup | f |.
Somando sobre k,
n
C
∑ ∆k E[| X1 |3 ] + E[| N1 |3 ] = O √1n → 0.
≤ n·
k =1
n3/2
Logo,
E[ f ( An )] − E[ f ( Bn )] → 0,
d
e como Bn ∼ N (0, 1) para todo n, segue que An → N (0, 1).
Capítulo 3
O ponto de partida do nosso curso será sempre o mesmo: só podemos utilizar variáveis unifor-
mes para gerar todas as demais distribuições. Ou seja, assumimos que temos disponível uma
variável aleatória
U ∼ Uniforme(0, 1),
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 X assuma
valores x1 , x2 , . . . , xm com probabilidades p1 , p2 , . . . , pm , onde
m
p k = P( X = x k ), pk ≥ 0, ∑ pk = 1.
k =1
k
Fk = ∑ pi , k = 1, . . . , m.
i =1
3. Retornar X = xk .
17
18 CAPÍTULO 3. VARIÁVEIS DISCRETAS E COMO SIMULÁ-LAS
A propriedade P( a < U < b) = b − a garante que
P( X = x k ) = p k .
Figura 3.1: Particionamento do intervalo (0, 1) para simular uma variável Bernoulli com p = 0.7.
Sorteia-se U ∼ Uniforme(0, 1); se U cair 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 3.2.
Figura 3.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 X uma 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
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
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.
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 p j e identificar onde
U caiu.
20 CAPÍTULO 3. VARIÁVEIS DISCRETAS E COMO SIMULÁ-LAS
p X ( x k ) = P( X = x k ), xk ∈ S,
a qual satisfaz
p X (k) ≥ 0 para todo k ∈ S, ∑ pX (k) = 1.
k∈S
X = ⌊nU ⌋ + 1,
3.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) = p e P( X = 0) = 1 − p,
3.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
Exercício 11. Prove as propriedades acima, isto é, calcule a esperança e a variância de uma variável
Bernoulli.
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.
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 X pode ser escrito como
n
X= ∑ Bi ,
i =1
22 CAPÍTULO 3. VARIÁVEIS DISCRETAS E COMO SIMULÁ-LAS
onde B1 , B2 , . . . , Bn são variáveis independentes e identicamente distribuídas, cada uma com
Bi ∼ Bernoulli( p).
2. Retornar X = ∑in=1 Bi .
Em outras palavras, uma variável binomial conta o número de sucessos em n tentativas inde-
pendentes, cada uma com probabilidade de sucesso p. Portanto, simular uma binomial se reduz
a repetir n vezes o procedimento de simulação da Bernoulli e somar os resultados.
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
n i
P( X = i ) = p (1 − p ) n −i , i = 0, 1, . . . , n.
i
observamos que 1
n−i n−i n
n n! n!
= = = .
i+1 ( i + 1) ! ( n − i − 1) ! i + 1 i! (n − i )! i+1 i
Substituindo essa relação,
n−i n
P( X = i + 1) = p i +1 (1 − p ) n − i −1 .
i+1 i
Reorganizando,
n−i p
P( X = i + 1) = · P( X = i ).
i+1 1− p
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:
n−i p
p i +1 = · pi , i ← i + 1, F ← F + pi ;
i+1 1− p
4. Retornar X = i.
5 − 0 0.3
p1 = · · 0.16807 ≈ 0.36015, F = 0.16807 + 0.36015 = 0.52822.
1 0.7
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.
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 n e 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 n sorteios 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 n Bernoullis.
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 3. VARIÁVEIS DISCRETAS E COMO SIMULÁ-LAS
o valor típico da variável X també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 3.3: Comparação de tempo de execução (em segundos) entre o método da inversão recur-
siva e a soma de Bernoullis para n = 100 e N = 2000 amostras, variando p.
A Figura 3.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 n e não é afetado por p, o tempo do método da inversão cresce proporcionalmente a np,
aumentando de forma acentuada à medida que p se 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.
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 X sorteado é exatamente o índice k em 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.
E[S] = E[ X + 1] = E[ X ] + 1.
3.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,
P( X = k ) = (1 − p)k−1 p, k = 1, 2, 3, . . .
Nesse caso:
1 1− p
E[ X ] = , Var( X ) = .
p p2
Exercício 14. Prove que a função de probabilidade acima satisfaz ∑∞
k =1 P( X = k ) = 1.
1. Inicializar o contador X ← 1;
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 tenta-
tivas são necessárias até que ocorra o primeiro sucesso.
26 CAPÍTULO 3. VARIÁVEIS DISCRETAS E COMO SIMULÁ-LAS
Exemplo 4. Se p = 0.3, então a probabilidade de obter um sucesso logo na primeira tentativa é 0.3. Se a
primeira tentativa falha, a segunda terá chance 0.3, e assim por diante.
Suponha que, ao simular, os primeiros valores de Bernoulli gerados foram 0, 0, 1. Isso indica duas falhas
seguidas e um sucesso na terceira tentativa. Portanto, o algoritmo retorna X = 3.
Assim, podemos usar o método da inversão para gerar X. Seja U ∼ Uniforme(0, 1). Definimos
X = j se
1 − (1 − p ) j −1 ≤ U < 1 − (1 − p ) j ,
ou seja,
(1 − p ) j < 1 − U ≤ (1 − p ) j −1 ,
o que equivale a
X = min{ j : (1 − p) j < 1 − U }.
Logo,
n log(1 − U ) o
X = min j : j > .
log(1 − p)
Portanto, obtemos a fórmula fechada
log(1 − U )
X = + 1.
log(1 − p)
Como 1 − U ∼ Uniforme(0, 1), podemos substituir 1 − U por U sem perda de generalidade,
resultando em
log(U )
X = + 1.
log(1 − p)
Exemplo 5. Considere uma variável geométrica X ∼ Geom( p) com p = 0.3 e seja U = 0.52 uma
realização de uma variável uniforme (0, 1). Usando a fórmula fechada da inversão, temos
log(1 − U )
X = + 1,
log(1 − p)
3.5. DISTRIBUIÇÃO DE POISSON 27
e, substituindo os valores, obtemos
log(0.48)
X = + 1 ≈ ⌊2.06⌋ + 1 = 3.
log(0.7)
F ( j ) = P( X ≤ j ) = 1 − (1 − p ) j .
Como F (2) = 0.51 < U = 0.52 mas F (3) = 0.657 ≥ 0.52, o menor j que satisfaz é j = 3.
e−λ λk
P( X = k ) = , k = 0, 1, 2, . . . ,
k!
onde λ > 0 representa a taxa média de ocorrências no intervalo considerado.
A média e a variância são dadas por
E[ X ] = λ, Var( X ) = λ.
Exercício 17. Prove as propriedades acima. Dica: para a variância, calcule primeiro E[ X ( X − 1)] e use o
fato de que
2
Var( X ) = E[ X ( X − 1)] + E[ X ] − E[ X ] .
28 CAPÍTULO 3. VARIÁVEIS DISCRETAS E COMO SIMULÁ-LAS
3.5.1 Simulação a Poisson via inversão e recursão
λi
P( X = i ) = e − λ , i = 0, 1, 2, . . .
i!
e essa expressão satisfaz a relação de recorrência
λ
P( X = i + 1) = P( X = i ).
i+1
i +1
Para ver isso, basta observar que P( X = i + 1) = e−λ (iλ+1)! . Separando um fator λ/(i + 1),
i
obtemos P( X = i + 1) = i+λ1 e−λ λi! , que nada mais é do que i+λ1 P( X = i ). Assim, conhecendo
p0 = P( X = 0) = e−λ , é possível calcular recursivamente p1 = λp0 , depois p2 = (λ/2) p1 , e
assim sucessivamente. Esse raciocínio evita a recomputação de fatoriais a cada passo e fornece
um procedimento numericamente mais estável.
O algoritmo clássico para gerar uma variável de Poisson com parâmetro λ funciona da se-
guinte maneira:
2. Inicializar i = 0, p0 = e−λ e F = p0 ;
λ
i ← i + 1, pi ← p i −1 , F ← F + pi ;
i
4. Retornar X = i.
Uma forma mais eficiente de implementar o método é iniciar a busca em torno do valor mais
provável da variável, que está próximo de λ. Seja m = ⌊λ⌋. Calcula-se a probabilidade acumulada
F ( m ) = P( X ≤ m ),
T = 1 + | X − m |.
30 CAPÍTULO 3. VARIÁVEIS DISCRETAS E COMO SIMULÁ-LAS
Como m ≈ λ, o custo médio pode ser aproximado por
E[ T ] = 1 + E[| X − λ|].
A distribuição de Poisson pode ser vista como um caso limite da distribuição binomial.
Seja X ∼ Binomial(n, p), que modela o número de sucessos em n tentativas independentes,
cada uma com probabilidade p de sucesso. Suponha agora que
Nesse regime, dizemos que a binomial entra no chamado limite de Poisson, e temos
d
Binomial(n, p) −
→ Poisson(λ).
n k
P( X = k ) = p (1 − p ) n − k .
k
λ n−k
k
n λ
P( X = k ) = 1− .
k n n
Para analisar o limite, consideramos cada fator separadamente. O coeficiente binomial pode
ser escrito como
n ( n − 1) · · · ( n − k + 1)
n
= .
k k!
Dividindo numerador e denominador por nk , temos
n k n ( n − 1) · · · ( n − k + 1) nk n n − 1 n−k+1
n
= · k
= · ··· .
k k! n k! n n n
3.6. DISTRIBUIÇÃO BINOMIAL NEGATIVA 31
Cada termo do produto no numerador pode ser escrito como
n−j j
= 1− , j = 0, 1, . . . , k − 1,
n n
e portanto
nk k−1
n 1 2
= · 1− 1− ··· 1− .
k k! n n n
Quando n → ∞, cada termo do produto tende a 1, de modo que
nk
n
∼ .
k k!
Além disso,
n −k
λ −λ λ
1− −→ e , 1− −→ 1
n n
Juntando os três resultados, obtemos
k
nk λ λk −λ
P( X = k ) ∼ e−λ · 1 = e .
k! n k!
Portanto,
e−λ λk
lim P( X = k ) =
,
n→∞ k!
que é exatamente a função de massa de probabilidade da distribuição Poisson(λ).
X = X1 + X2 + · · · + X r ,
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).
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, X representa o número total de ensaios necessários até a ocorrência de
r sucessos. Assim, o algoritmo pode ser descrito da seguinte forma:
2. Enquanto s < r:
3. Retornar X = n.
X = X1 + X2 + · · · + X r ,
Xi ∼ Geom( p), i = 1, . . . , r,
log(1 − Ui )
Xi = + 1, Ui ∼ Uniforme(0, 1).
log(1 − p)
2. Retornar X = ∑ri=1 Xi .
3.6. DISTRIBUIÇÃO BINOMIAL NEGATIVA 33
3.6.3 Simulando via inversão recursiva
Outra forma de simular a Binomial Negativa é aplicar diretamente o método da inversão, apro-
veitando a relação de recorrência da sua função de probabilidade.
Se X ∼ NegBin(r, p), então
n−1 r
P( X = n ) = p (1 − p ) n −r , n = r, r + 1, r + 2, . . .
r−1
P( X = n + 1) n
= (1 − p ).
P( X = n ) n−r+1
2. Inicializar n = r, pn = pr , F = pn ;
n
p n +1 = p n · (1 − p ), n ← n + 1, F ← F + p n +1 ;
n−r+1
4. Retornar X = n.
T = (n − r ) + 1,
r
E[ T ] = E[ X − r ] + 1 = − r + 1.
p
Esse termo −r aparece 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 3. VARIÁVEIS DISCRETAS E COMO SIMULÁ-LAS
3.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
n
n k n−k
1 = ( p + q) = ∑ n
p q .
k =0
k
Para ver a equivalência com a forma escrita em função do número total de ensaios n, deta-
lhamos a reparametrização. Defina n = r + k (isto é, k = n − r). Então
r+k−1 (r + k − 1) ! ( n − 1) ! n−1
= = = .
k k! (r − 1)! ( n − r ) ! (r − 1) ! n−r
r + (n − r ) − 1
P ( X = n ) = P (Y = n − r ) = p r (1 − p ) n −r
n−r
n−1
= p r (1 − p ) n −r , n = r, r + 1, . . .
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.
Exercício 22. Prove as fórmulas da média e variância acima. Dica: considere o sorteio sequencial das n
bolas e defina Xi como a variável indicadora do evento “a i-ésima bola é clara”. Para a variância, use a
decomposição
n
Var( X ) = ∑ Var(Xi ) + 2 ∑ Cov( Xi , X j ).
i =1 1≤ i < j ≤ n
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 n elementos distintos de forma aleatória. O algoritmo funciona assim:
Esse método garante que cada subconjunto de tamanho n tem a mesma probabilidade de ser
escolhido, e é mais eficiente do que embaralhar toda a população.
36 CAPÍTULO 3. VARIÁVEIS DISCRETAS E COMO SIMULÁ-LAS
Tabela de referência
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 rejeição ou aceitação: baseiam-se em gerar propostas e aceitar com certa
probabilidade;
Esse procedimento garante que X terá exatamente a distribuição desejada, pois a probabi-
lidade de X cair em qualquer intervalo será proporcional ao comprimento correspondente no
domínio de U.
Esse método é particularmente útil quando a inversa de F pode ser escrita de forma explícita,
como ocorre com as distribuições Exponencial, Uniforme e Pareto, por exemplo.
37
38 CAPÍTULO 4. VARIÁVEIS CONTÍNUAS E COMO SIMULÁ-LAS
4.1.1 Distribuição exponencial
FX ( x ) = P( X ≤ x ) = 1 − e−λx , x ≥ 0.
1 1
E[ X ] = , Var( X ) = .
λ λ2
R∞
Exercício 23. Verifique que f X ( x ) é uma densidade de probabilidade, isto é, 0
f X ( x ) dx = 1.
A distribuição exponencial é um exemplo clássico onde o método da inversão pode ser apli-
cado diretamente. Sabemos que a CDF é
F ( x ) = 1 − e−λx ,
1
1 − e−λx = U ⇒ x=− log(1 − U ).
λ
Como 1 − U ∼ Uniforme(0, 1), podemos reescrever de forma equivalente:
1
X=− log(U ), com U ∼ Uniforme(0, 1).
λ
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 X como o número de tentativas
independentes até a ocorrência do primeiro sucesso, em uma sequência de ensaios de Bernoulli
com probabilidade p de sucesso.
4.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
Xn
Tn = .
n
A variável Tn representa o tempo até o primeiro sucesso quando fazemos n tentativas por uni-
dade 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 espe-
rado de sucessos por unidade de tempo permanece constante: n · pn = λ.
Vamos mostrar que Tn converge em distribuição para uma variável exponencial de parâmetro
λ. De fato, temos:
Xn
P( Tn > t) = P > t = P( Xn > ⌊nt⌋).
n
Como Xn é geométrica com parâmetro pn = λ/n, segue que:
⌊nt⌋
λ
P( Xn > k ) = (1 − p n ) ,
k
logo P( Tn > t) = 1− .
n
Portanto,
P( Tn ≤ t) → 1 − e−λt ,
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á:
1
T=− ln(U ).
λ
Já no caso da geométrica Xn ∼ Geom( pn ), a fórmula de inversão baseada na CDF discreta é:
ln(U ) λ
Xn = , com pn = .
ln(1 − pn ) n
Dividindo por n, temos:
Xn 1 ln(U )
Tn = ≈ · .
n n ln(1 − λ/n)
Sabemos que para n grande,
λ
ln(1 − λ/n) ≈ − ,
n
então:
1
Tn ≈ −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.
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 n subintervalos de comprimento 1/n. Em
cada subintervalo, ocorre um evento (ou sucesso) com probabilidade pn = λ/n, de forma inde-
pendente. Este é exatamente o modelo da variável binomial
Xn ∼ Binomial(n, λ/n),
Xn
Tn = .
n
Como visto anteriormente, temos
d
Tn −
→ Exponencial(λ).
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 g e aceitá-los com uma certa probabilidade
que depende da razão f ( x )/g( x ).
42 CAPÍTULO 4. VARIÁVEIS CONTÍNUAS E COMO SIMULÁ-LAS
Suponha que desejamos gerar uma variável aleatória X com densidade alvo f ( x ), mas não
dispomos de um método direto para isso. Por outro lado, assumimos que sabemos simular uma
variável Y com densidade auxiliar g( x ), e que existe uma constante c > 0 tal que
f (x)
≤c para todo x.
g( x )
Essa condição garante que a função f está sempre abaixo da curva cg, ou seja, f ( x ) ≤ cg( x )
f (x)
para todo x. Além disso, isso assegura que a razão está sempre entre 0 e 1, podendo ser
cg( x )
interpretada como uma probabilidade de aceitação. Note que, ao integrar ambos os lados da
desigualdade f ( x ) ≤ cg( x ), obtemos f ( x ) dx ≤ cg( x ) dx, ou seja, 1 ≤ c, e portanto 1c ≤ 1.
R R
1. Gere um candidato Y ∼ g.
Para entender por que o método de rejeição-aceitação funciona, vamos construir uma intuição
passo a passo com um exemplo concreto. Suponha que queremos gerar uma variável aleatória
X ∼ f , com densidade definida por
Essa é uma densidade válida sobre o intervalo [0, 1], mas sua função de distribuição acumulada
F ( x ) não possui inversa em forma fechada já que envolve resolver uma equação polinomial de
grau 5, o que inviabiliza o uso direto do método da inversão. Por isso, recorremos ao método de
rejeição.
Nesse caso, utilizamos como densidade auxiliar a uniforme g( x ) = 1 sobre [0, 1], que é fácil
de simular. O procedimento funciona da seguinte forma:
1. Escolhemos uma constante c > 0 tal que f ( x ) ≤ cg( x ) para todo x ∈ [0, 1]. Como g( x ) = 1,
essa condição se torna f ( x ) ≤ c. Para garantir isso, basta determinar o valor máximo da
função f ( x ) no intervalo [0, 1], o que pode ser feito derivando:
3. Geramos um valor U ∼ Unif(0, 1), que usaremos para introduzir variabilidade vertical.
U · cg(Y ) < f (Y ),
f (x)
≤c para todo x ∈ X .
g( x )
Considere o algoritmo de geração:
Então a variável aleatória X, definida como o primeiro valor Y aceito, possui densidade f . Além disso,
o número total de iterações até a aceitação segue uma distribuição geométrica com parâmetro 1/c.
De forma análoga:
f (Y ) f (Y )
Z
P U< = P U< | Y = y g(y) dy
cg(Y ) cg(Y )
f (y)
Z
= P U< | Y = y g(y) dy
cg(y)
f (y)
Z
= P U< g(y) dy
cg(y)
f (y)
Z
= · g(y) dy
cg(y)
1 1
Z
= f (y) dy = .
c c
Substituindo numerador e denominador:
1 x
R
c −∞ f ( y ) dy
Z x
P( X ≤ x ) = 1
= f (y) dy = FX ( x ).
c −∞
ou seja, cada tentativa tem probabilidade 1/c de ser aceita. Portanto, o número de repetições até
obter um ponto aceito segue uma distribuição geométrica com parâmetro 1/c.
A distribuição normal padrão, denotada por N (0, 1), é uma das distribuições mais importantes
da estatística e da probabilidade. Sua densidade é dada por:
1 2
f ( x ) = √ e− x /2 , x ∈ R.
2π
2 2
f ( x ) = √ e− x /2 ,
2π
2 /2
e que, para x > 0, a função decresce com a cauda e− x . Por outro lado, a densidade da
distribuição exponencial com taxa λ = 1 é
g( x ) = e− x , x ≥ 0.
Note que:
r
f (x) 2 − x2 /2+ x
= e .
g( x ) π
Para aplicar o método da rejeição, precisamos encontrar o ponto de máximo da razão f ( x )/g( x ),
ou seja, maximizar a função x − x2 /2. Derivando:
x2
d
x− = 1 − x = 0 ⇒ x = 1.
dx 2
1. Gere Y ∼ Exp(1).
Finalmente, se X ∼ |N (0, 1)|, isto é, uma normal padrão truncada para valores positivos
(obtida via o algoritmo anterior).
Para gerar uma normal padrão simétrica Z ∼ N (0, 1), basta sortear um sinal S ∼ Bernoulli(1/2),
e definir:
X, se S = 1,
Z=
− X, se S = 0.
Dessa forma, Z tem distribuição simétrica em torno de zero, com densidade normal padrão,
como desejado.
Intuição geométrica
O método da rejeição pode ser visualizado como um processo de amostragem de pontos aleatórios
em uma região do plano, com o objetivo de “pintar” a curva da densidade alvo f ( x ).
Imagine que temos uma função auxiliar g( x ), da qual sabemos simular facilmente, e uma
constante de majoração c > 0 tal que f ( x ) ≤ cg( x ) para todo x. Isso nos permite usar cg( x )
como um envelope que cobre toda a curva de f ( x ).
A cada tentativa, sorteamos:
f (Y )
U< .
cg(Y )
1
f (x) = x α−1 e− x/θ , x > 0,
Γ(α) θ α
Uma maneira natural de interpretar essa distribuição é por analogia com modelos discretos.
No mundo discreto, a distribuição Geométrica mede o número de ensaios necessários até obser-
var o primeiro sucesso em uma sequência de Bernoullis. Se quisermos o número de ensaios até
o r-ésimo sucesso, obtemos a Binomial Negativa, que pode ser vista como a soma de variáveis
Geométricas independentes.
• Quando α é um inteiro maior que 1, a Gamma pode ser entendida como a soma de α
Exponenciais independentes.
O parâmetro de forma α regula quantos sucessos estão sendo acumulados e, portanto, afeta
diretamente a forma da distribuição:
E[ X ] = αθ, Var( X ) = αθ 2 .
Considere X ∼ Gamma(α, λ) com α > 1 (parametrização por taxa λ). A densidade alvo é
λα α−1 −λx
f (x) = x e , x > 0.
Γ(α)
4.3. DISTRIBUIÇÃO GAMMA 49
g( x ) = µe−µx , x > 0.
Para que o método de aceitação–rejeição seja válido, precisamos de uma constante c tal que
f ( x ) ≤ c g( x ) para todo x > 0. O quociente
f (x) λα
= x α−1 e−(λ−µ)x
g( x ) Γ(α) µ
mostra que é necessário ter µ < λ, pois caso contrário o termo exponencial não decai e o quoci-
ente não tem máximo finito. Quando µ < λ, o máximo ocorre em
α−1
x∗ = ,
λ−µ
com valor
α −1
α−1
λα
c(µ) = e−(α−1) .
Γ(α) µ λ−µ
Para 0 < α < 1, uma forma simples de simular Γ(α, θ ) é usar a identidade
Para entender essa relação, considere as variáveis independentes (U, G ) com U ∼ Unif(0, 1)
e G ∼ Γ(α + 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 ) é
gα e− g/θ
f U,G (u, g) = f U (u) f G ( g) = 1(0,1) (u) , g > 0.
Γ ( α + 1 ) θ α +1
4.4. DISTRIBUIÇÃO BETA 51
Pela fórmula de mudança de variável,
∂(u, g)
f X,G ( x, g) = f U,G ( x/g)α , g det .
∂( x, g)
∂u ∂u ∂g ∂g
= α x α −1 g − α , = −α x α g−(α+1) , = 0, = 1,
∂x ∂g ∂x ∂g
logo
∂(u, g) ∂u
det = · 1 − 0 = α x α −1 g − α .
∂( x, g) ∂x
Portanto,
gα e− g/θ α −1 − α α x α −1
f X,G ( x, g) = 1(x,∞) ( g) α x g = 1 ( x,∞ ) ( g ) e− g/θ .
Γ ( α + 1 ) θ α +1 Γ ( α + 1 ) θ α +1
x α−1 e− x/θ
f X (x) = , x > 0,
Γ(α) θ α
1
f (x) = x α −1 (1 − x ) β −1 , 0 < x < 1,
B(α, β)
onde
Γ(α)Γ( β)
Z 1
B(α, β) = uα−1 (1 − u) β−1 du =
0 Γ(α + β)
é a função Beta de Euler.
52 CAPÍTULO 4. VARIÁVEIS CONTÍNUAS E COMO SIMULÁ-LAS
A interpretação intuitiva da distribuição Beta é como um modelo de incerteza sobre proba-
bilidades. Se pensamos em x como 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 s sucessos e f fracassos, a
posteriori é Beta(α + s, β + f ).
α αβ
E[ X ] = , Var( X ) = .
α+β ( α + β )2 ( α + β + 1)
Essas fórmulas mostram como α e β podem ser interpretados como “pseudocontagens” de suces-
sos e fracassos, de forma que α + β controla a concentração da distribuição em torno da média.
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 c tal que
f (y) ≤ c g(y) = c, 0 < y < 1.
O algoritmo de aceitação–rejeição é:
O método mais utilizado e geral para simular variáveis Beta(α, β) explora a relação entre as
distribuições Beta e Gama. Seja
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.
1 −t
f U,T (u, t) = e , u ∈ (0, 2π ), t > 0.
2π
Definimos a transformação
√ √
X= 2T cos U, Y= 2T sin U.
de modo que
y
t = 12 ( x2 + y2 ), u = arctan x (ajustado para o quadrante correto).
∂x √ ∂x 1
= − 2t sin u, = √ cos u,
∂u ∂t 2t
∂y √ ∂y 1
= 2t cos u, = √ sin u.
∂u ∂t 2t
Logo, √
− 2t sin u √1 cos u
2t
J=√ .
2t cos u √1 sin u
2t
O determinante é
√ √
det( J ) = − 2t sin u √1 sin u − √1 cos u 2t cos u .
2t 2t
4.5. TRANSFORMAÇÕES DE VARIÁVEIS ALEATÓRIAS 55
Simplificando,
det( J ) = − sin2 u − cos2 u = −1.
Portanto,
|det( J )| = 1.
∂(u, t) 1 −t
f X,Y ( x, y) = f U,T (u, t) det = e · 1.
∂( x, y) 2π
Substituindo t = 21 ( x2 + y2 ),
1
exp − 21 ( x2 + y2 ) .
f X,Y ( x, y) =
2π
Finalmente, notamos que
1 2 1 −y2 /2
f X,Y ( x, y) = √ e−x /2 √ e ,
2π 2π
o que mostra que X e Y sã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 4. VARIÁVEIS CONTÍNUAS E COMO SIMULÁ-LAS
2. Gere T ∼ Expo(1).
√ √
3. Calcule X = 2T cos U e Y = 2T sin U.
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
q
Z = X, W = ρX + τY, τ = 1 − ρ2 .
a inversa é
w − ρz
x = z, y= .
τ
Aplicando a fórmula de mudança de variáveis,
∂( x, y)
f Z,W (z, w) = f X,Y ( x, y) det ,
∂(z, w)
portanto
∂( x, y) 1
det = .
∂(z, w) τ
Substituindo em f Z,W ,
1 1 w−ρz 2
f Z,W (z, w) = exp − z2 + .
2πτ 2 τ
4.5. TRANSFORMAÇÕES DE VARIÁVEIS ALEATÓRIAS 57
Fazendo as contas e lembrando que ρ2 + τ 2 = 1, obtemos
1 1 2 2
f Z,W (z, w) = exp − 2 (z − 2ρzw + w ) .
2πτ 2τ
Essa é exatamente a forma conhecida da densidade Normal bivariada com matriz de covari-
ância !
1 ρ
Σ= , | Σ | = 1 − ρ2 = τ 2 .
ρ 1
De fato, podemos escrever
1
−1 ⊤
f Z,W (z, w) = exp − 1
( z, w ) Σ ( z, w ) .
2π |Σ|1/2 2
Do ponto de vista de simulação, esse resultado mostra que basta gerar X, Y ∼ N (0, 1) inde-
pendentes (e.g. via Box–Muller) e aplicar a transformação acima. O algoritmo é:
3. Calcule q
Z = X, W = ρX + 1 − ρ2 Y.
Fazendo as contas,
1 0 1 ρ 1 ρ
AA⊤ = p p = ,
ρ 1 − ρ2 0 1 − ρ2 ρ 1
X = ( X1 , . . . , Xd )⊤ ∼ Nd (0, Id ),
Z = AX.
Então
Cov( Z ) = A Cov( X ) A⊤ = AId A⊤ = AA⊤ .
Portanto, dado Σ simétrica definida positiva, basta encontrar A tal que Σ = AA⊤ . Assim,
podemos gerar
Z ∼ Nd (0, Σ).
Z = µ + AX.
Para encontrar A a partir de uma matriz de covariância Σ simétrica definida positiva, existem
diferentes decomposições possíveis:
Σ = LL⊤ ,
Σ = QΛQ⊤ ,
A = QΛ1/2 ,
com
p p
Λ1/2 = diag( λ1 , . . . , λ d ).
0.3 0.5 1
4.5. TRANSFORMAÇÕES DE VARIÁVEIS ALEATÓRIAS 59
Para gerar amostras desse vetor, seguimos a ideia de representar uma Normal multivariada como
transformação linear de Normais independentes. Seja
X = ( X1 , X2 , X3 )⊤ ∼ N3 (0, I3 ),
Σ = LL⊤ ,
então
Z = LX
terá distribuição N3 (0, Σ).
Z3 X3
Dizemos que Q segue a distribuição qui-quadrado com k graus de liberdade, denotada por
Q ∼ χ2k .
60 CAPÍTULO 4. VARIÁVEIS CONTÍNUAS E COMO SIMULÁ-LAS
A densidade dessa distribuição pode ser derivada observando que Zi2 ∼ Γ 12 , 2 , e que a soma
Logo, a densidade é
1 k
f (q) = q 2 −1 e−q/2 , q > 0.
2k/2 Γ(k/2)
Os principais momentos são
E[ Q] = k, Var( Q) = 2k.
Essa distribuição aparece com frequência em estatística, por exemplo em testes de hipóteses e
intervalos de confiança, pois estatísticas do tipo “soma de quadrados de erros padronizados”
têm exatamente essa forma.
Então
k
Q= ∑ Zi2 ∼ χ2k .
i =1
O algoritmo é:
1. Fixe k ∈ N.
i.i.d.
2. Gere Z1 , . . . , Zk ∼ N (0, 1).
Outra forma de simular uma variável χ2k é usar sua equivalência com a distribuição Gama. Sabe-
mos que
Q ∼ χk ⇐⇒ Q ∼ Γ 2 , 2 ,
2 k
O algoritmo é:
4.5. TRANSFORMAÇÕES DE VARIÁVEIS ALEATÓRIAS 61
1. Fixe k ∈ N.
2. Gere G ∼ Γ k
2, 2 .
3. Retorne Q = G.
Esse método é geralmente mais eficiente em termos computacionais, já que muitas bibliote-
cas numéricas possuem rotinas otimizadas para a geração de variáveis Gama com parâmetros
arbitrários. Assim, para valores grandes de k, pode ser preferível usar diretamente a simulação
via Gama em vez da soma de muitos Normais.
A distribuição t de Student foi introduzida em 1908 por William Gosset, que trabalhava como
mestre cervejeiro na Guinness. Por restrições da empresa, ele publicou seus resultados sob o
pseudônimo Student, dando origem ao nome da distribuição. Essa distribuição aparece natural-
mente em problemas de inferência estatística, especialmente em testes de hipóteses, mas aqui
nos interessa sua definição e como simulá-la.
T ∼ tν .
A intuição da fórmula pode ser entendida a partir da estatística de teste para a média de uma
Normal. Se a variância σ2 fosse conhecida, teríamos
X−µ
Z= √ ∼ N (0, 1).
σ/ n
A ideia fundamental da simulação de Monte Carlo é usar amostras aleatórias para aproximar
quantidades numéricas que, de outra forma, seriam difíceis ou impossíveis de calcular analitica-
mente. Em sua forma mais simples, o método baseia-se na Lei dos Grandes Números.
µ = E [h( X )] .
1 n
n i∑
µ̂n = h ( Xi ) .
=1
Pela Lei dos Grandes Números, µ̂n → µ quase certamente quando n → ∞. Uma maneira
simples de justificar essa convergência é pela desigualdade de Chebyshev. Se Var [h( X )] = σ2 < ∞,
então
σ2
Var [µ̂n ] = .
n
Logo, para qualquer ε > 0,
Var [µ̂n ] σ2
P (|µ̂n − µ| > ε) ≤ = .
ε2 nε2
Portanto, P (|µ̂n − µ| > ε) → 0 quando n → ∞, mostrando que µ̂n converge para µ em probabi-
lidade, o que é precisamente a versão fraca da Lei dos Grandes Números.
5.1.1 Exemplos
63
64 CAPÍTULO 5. SIMULAÇÃO VIA MONTE CARLO
Geramos X1 , . . . , Xn ∼ Uniforme(0, 1) e usamos a identidade
h 2
i
I = E e− X .
Exemplo 10. Considere X ∼ N (0, 1) e o evento A = { X > 1}. Queremos estimar a probabilidade
p = P ( X > 1) .
1 n
1 { Xi > 1 } .
n i∑
p̂n =
=1
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 1 centrado na origem, definido por
x2 + y2 ≤ 1.
1 n
1{ Xi2 + Yi2 ≤ 1}.
n i∑
π̂n = 4 ×
=1
Exemplo 12. Considere X ∼ Uniforme(0, 1). Queremos estimar simultaneamente E [ X ] e Var [ X ] por
simulação.
Geramos X1 , . . . , Xn independentes e usamos
1 n 1 n
n i∑ ∑
µ̂n = Xi , σ̂n2 = ( Xi − µ̂n )2 .
=1
n − 1 i =1
5.1. ESTIMANDO MÉDIAS 65
Exemplo 13. Considere a integral em duas dimensões
Z 1Z 1
2 + y2 )
I= e−(x dx dy.
0 0
1 n −(X2 +Y2 )
n i∑
În = e i i .
=1
pendentes e calculamos
1 n
µ̂n = ∑ e−Xi .
n i =1
O valor verdadeiro é E e
−X
= 21 .
1 n 1
În = ∑
n i = 1 1 + Xi
.
Assim, a integral pode ser estimada por Monte Carlo sem precisar integrar diretamente uma função
oscilatória:
1 n sin( Xi )
În = ∑ , Xi ∼ Exponencial(1).
n i = 1 X i e − Xi
sn = 0.289.
sn 0.289
EP = √ = √ = 0.000914.
n 100000
5.2. INTERVALOS DE CONFIANÇA 67
(4) Pela regra empírica 68–95–99.7, sabemos que:
Assim, podemos construir um intervalo aproximado de 95% de confiança usando dois desvios padrão
em vez de 1.96.
2 × EP = 2 × 0.000914 = 0.00183.
Em outras palavras, esperamos que cerca de 95% das repetições do experimento de Monte Carlo produ-
zam valores de În dentro de dois erros padrão da média verdadeira. O valor teórico I = 0.7468 está de fato
dentro desse intervalo.
Redução de variância
E[ X ] = θ.
69
70 CAPÍTULO 6. REDUÇÃO DE VARIÂNCIA
1
• se X1 e X2 forem independentes, Cov( X1 , X2 ) = 0 e Var(θ̂ ) = 2 Var( X1 );
A questão, então, é: como gerar dois valores X1 e X2 com a mesma distribuição, mas negati-
vamente correlacionados?
X1 = h(U1 , . . . , Um ),
X2 = h(1 − U1 , . . . , 1 − Um ),
X1 = h ( U ) , X2 = h ( 1 − U ) .
X1 = U, X2 = 1 − U.
E[ X1 ]E[ X2 ] = 41 .
X1 = h(U1 , . . . , Um ), X2 = h(1 − U1 , . . . , 1 − Um ),
6.1. USO DE VARIÁVEIS ANTITÉTICAS 71
com h crescente em cada coordenada.
Nesse caso, X1 é uma função crescente do vetor (U1 , . . . , Um ), enquanto X2 é decrescente.
Considerando
g(U1 , . . . , Um ) = − X2 ,
o que implica
Cov( X1 , X2 ) ≤ 0.
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 inde-
pendentes X1 e X2 , construímos o par de forma que ambas tenham a mesma distribuição, mas
sejam negativamente correlacionadas. A variável antitética é dado por
X1 + X2
X′ = ,
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:
2. Calcule X1 = h(U1 , . . . , Um ).
3. Calcule X2 = h(1 − U1 , . . . , 1 − Um ).
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 − U també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 U1 e U2 para obter
duas amostras de X, basta gerar uma única variável uniforme U. Com ela, obtemos simultanea-
mente 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 6. 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
1
X=− 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
1 1
X1 = − log(U ), X2 = − log(1 − U ).
λ λ
Definimos, então, a variável final como a média
X1 + X2
Z= .
2
O algoritmo é:
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 = − 21 log U (1 − U ) .
6.1. USO DE VARIÁVEIS ANTITÉTICAS 73
Como U (1 − U ) ≤ 1/4, segue que
1
Z≥ 2 log 4 = log 2 ≈ 0.693.
Ou seja, a variável construída por antitéticos nunca assume valores menores que log 2.
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.
Observe que o termo e− x corresponde à densidade de uma variável X ∼ Exp(1). Assim, podemos reescrever
a integral como
I = E log(1 + X 2 ) , X ∼ Exp(1).
Portanto, a solução via Monte Carlo é imediata: basta gerar amostras Xi ∼ Exp(1), calcular log(1 +
Xi2 )e tirar a média. O algoritmo segue os passos:
2. Transformar em Xi = − log(Ui ).
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
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.
J = E[ e Z ], Z ∼ N (0, 1).
2. Calcular e Zi .
Note que a distribuição normal é simétrica em torno de zero, isto é, Z ∼ N (0, 1) implica que também
− Z ∼ N (0, 1). Assim, para cada Zi gerado, podemos considerar o par ( Zi , − Zi ) e formar o estimador
1 n 1 Zi
Ĵant = ∑
n i =1 2
e + e − Zi
.
Neste caso, como e x é uma função monotônica crescente, os valores e Zi e e−Zi tendem a se compen-
sar, gerando correlação negativa e uma variância muito menor na estimativa. O ponto essencial é que a
esperança se mantém inalterada, mas o uso da antitética torna o estimador mais eficiente.
onde X é o resultado de uma simulação. Agora suponha que exista outra variável Y cuja espe-
rança é conhecida, digamos
E [Y ] = µ Y .
Z = X + c (Y − µ Y )
Cov( X, Y )
c∗ = − .
Var(Y )
Cov( X, Y )2
Var( Z ) = Var( X ) − .
Var(Y )
Var( Z )
= 1 − ρ( X, Y )2 ,
Var( X )
onde
Cov( X, Y )
ρ( X, Y ) = p
Var( X ) Var(Y )
é a correlação entre X e Y. Isso mostra que a redução relativa de variância obtida é de 100 ·
ρ( X, Y )2 por cento, independentemente de a correlação ser positiva ou negativa.
76 CAPÍTULO 6. REDUÇÃO DE VARIÂNCIA
Na prática, Cov( X, Y ) e Var(Y ) não são conhecidos de antemão e precisam ser estimados a
partir dos dados simulados. Se n simulações são realizadas, gerando pares ( Xi , Yi ), podemos
calcular
n n
d ( X, Y ) = 1 ∑ ( Xi − X̄ )(Yi − Ȳ ), Var
Cov d (Y ) = 1 ∑ (Yi − Ȳ )2 ,
n − 1 i =1 n − 1 i =1
e então definir
d ( X, Y )
Cov
ĉ∗ = − .
d (Y )
Var
O estimador final com variável de controle é dado por
1 n
n i∑
∗
θ̂ctrl = X i + ĉ ( Yi − µ Y ) .
=1
2. Calcular X̄ e Ȳ.
d ( X, Y ) e Var
3. Estimar Cov d (Y ) .
θ = E[ eU ].
Assim, definindo X = eU , podemos estimar θ por Monte Carlo a partir da média amostral de X.
Agora considere a variável Y = U, cuja esperança é conhecida, µY = E[U ] = 0.5. Como X = eU
e Y = U são fortemente correlacionados positivamente, podemos utilizar Y como variável de controle. O
estimador controlado é dado por
Z = X + c ∗ (Y − µ Y ) ,
onde
Cov( X, Y )
c∗ = − .
Var(Y )
Na prática, basta gerar pares ( Xi , Yi ) a partir de Ui ∼ U (0, 1), calcular as estimativas amostrais de
Cov( X, Y ) e Var(Y ) para obter ĉ∗ , e então construir o estimador
1 n
n i∑
∗
θ̂ctrl = X i + ĉ ( Yi − µ Y ) .
=1
A correlação positiva entre X e Y faz com que o desvio de Y em relação à sua média indique a direção
do desvio de X em relação a θ. O termo de ajuste c∗ (Y − µY ) corrige esse efeito, reduzindo drasticamente
a variância do estimador. Na simulação, observamos que a variância caiu de aproximadamente 0.24 (sem
6.3. REDUÇÃO DE VARIÂNCIA POR CONDICIONAMENTO 77
controle) para 0.004 (com controle), com estimativas muito mais concentradas em torno do valor verdadeiro
e − 1 ≈ 1.718.
2. Calcular Xi = eUi e Yi = Ui .
d ( X, Y ) e Var
3. Estimar Cov d (Y ) a partir dos dados.
Portanto,
E E[ X 2 | Y ] = E Var( X | Y ) + E (E[ X | Y ])2 .
Substituindo de volta,
Var( X ) = E Var( X | Y ) + E[(E[ X | Y ])2 ] − (E[ X ])2 ,
78 CAPÍTULO 6. REDUÇÃO DE VARIÂNCIA
e o termo entre parênteses é precisamente Var(E[ X | Y ]).
Assim, chegamos à decomposição
Var( X ) ≥ Var E[ X | Y ] ,
Z = E[ X | Y ].
Se a pessoa sorteada for um homem, usamos a média de alturas dos homens (digamos, 175 cm);
se for uma mulher, usamos a média das mulheres (digamos, 162 cm). Note que a média global
continua correta: metade das vezes registramos 175, metade das vezes 162, o que resulta em
168,5 cm, exatamente a média real da população.
A identidade
Var( X ) = E[Var( X | Y )] + Var(E[ X | Y ])
mostra como essa substituição reduz a variância. O primeiro termo corresponde à variabilidade
dentro de cada grupo (diferenças entre indivíduos do mesmo sexo), enquanto o segundo termo
corresponde à variabilidade entre as médias dos grupos (diferença entre a média dos homens e a
das mulheres). Quando usamos diretamente X, ambos os termos estão presentes; quando usamos
E[ X | Y ], eliminamos o primeiro termo e ficamos apenas com a variabilidade entre grupos.
Assim, o estimador permanece não-viesado, mas com menor variância. Em outras palavras,
condicionar equivale a substituir um indivíduo ruidoso pela média de seu grupo, preservando a
esperança e reduzindo a dispersão.
como estimador de θ.
Esse algoritmo gera um estimador não-viesado de θ, mas com variância reduzida em compa-
ração ao estimador usual baseado diretamente em X.
Exemplo 22. Queremos estimar π via simulação. Podemos gerar dois números aleatórios U1 , U2 ∼
U (0, 1) e definir
Vi = 2Ui − 1, i = 1, 2,
então E[ I ] = π/4, pois a probabilidade de um ponto uniforme em [−1, 1]2 cair dentro do círculo unitário
é exatamente a razão entre a área do círculo e a área do quadrado. Assim, uma estimativa usual seria
4 n
n i∑
π̂ = Ii .
=1
Podemos, porém, melhorar esse estimador aplicando a técnica de condicionamento. Em vez de usar I
diretamente, consideremos E[ I | V1 ]. Temos
E[ I | V1 = v] = P V12 + V22 ≤ 1 V1 = v .
Isso equivale a
E[ I | V1 = v] = P V22 ≤ 1 − v2 .
80 CAPÍTULO 6. REDUÇÃO DE VARIÂNCIA
Como V2 ∼ U (−1, 1) e é independente de V1 , temos
p p
P V22 ≤ 1 − v2 = P − 1 − v2 ≤ V2 ≤ 1 − v2 .
Portanto,
Z √1− v2
1 p
E[ I | V1 = v] = √ dx = 1 − v2 .
− 1− v2 2
Assim, obtemos o estimador condicionado
q
Z= 1 − V12 ,
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 = 0 representa mulher e Y = 1 representa homem.
Condicionalmente a Y, a altura X tem distribuição normal
θ = E[ X ] = (1 − p)µ f + pµm .
1 n
n i∑
θ̂cond = Zi ,
=1
estritamente menor do que Var( X ), pois elimina a variabilidade interna de cada grupo (σ2f e σm2 ). 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 7
Uma alternativa é escolher uma outra densidade g( x ) tal que f ( x ) = 0 sempre que g( x ) = 0.
Nesse caso, podemos reescrever
f (x) f (X)
Z
θ = h( x ) g( x ) dx = Eg h( X ) ,
g( x ) g( X )
onde Eg denota esperança em relação à densidade g.
Assim, se gerarmos X1 , . . . , Xn ∼ g, um estimador natural é
1 n f ( Xi )
θ̂ = ∑
n i =1
h ( Xi )
g ( Xi )
.
81
82 CAPÍTULO 7. 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 probabilida-
des raras. Nesse caso, h( x ) é uma função indicadora de um conjunto A pouco provável sob f . Se
escolhermos g de modo que A seja mais frequente, então:
• 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.
etx f ( x )
f t (x) = .
M(t)
Intuitivamente, a densidade f t dá mais peso a valores grandes de X quando t > 0 e mais peso
a valores pequenos quando t < 0. Em muitos casos, f t pertence à mesma família paramétrica de
f , mas com parâmetros modificados.
Alguns exemplos:
• Exponencial. Se X ∼ Exp(λ), então f t é 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 f t é 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.
pet
pt = .
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.
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 esti-
mador. Uma alternativa é utilizar a técnica de amostragem por importância com densidades inclinadas.
Seja
etx f i ( x )
f i,t ( x ) = , Mi (t) = E[etXi ],
Mi ( t )
a densidade inclinada de Xi , onde t > 0 é um parâmetro comum a todas as variáveis. Ao simular cada Xi
segundo f i,t , obtemos que " #
n
θ = Et 1{S>a} exp(−tS) ∏ Mi (t) ,
i =1
( j) ( j)
onde S( j) = ∑in=1 Xi e cada Xi é simulado de f i,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.
2. Para j = 1, . . . , N:
( j) ( j)
(a) Gere X1 , . . . , Xn independentemente segundo as densidades inclinadas f i,t .
( j)
(b) Calcule S( j) = ∑in=1 Xi .
(c) Associe o peso
n
W ( j) = 1{S( j) >a} exp(−tS( j) ) ∏ Mi (t).
i =1
3. Estime θ por
N
1
θ̂ =
N ∑ W ( j) .
j =1
No caso particular em que cada Xi ∼ N (0, 1), temos que a soma S = ∑in=1 Xi é 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 é
etx f ( x ) 1 2
f t (x) = , M(t) = exp 2t ,
M(t)
o que implica que f t é 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 + 12 t2 ).
Nesse caso, a estimativa via amostragem por importância com N = 200,000 simulações foi
eλx f ( x )
f λ (x) = , Z (λ) = E f [eλX ].
Z (λ)
Podemos então expressar qualquer probabilidade como uma esperança sob essa nova medida:
f (x) f (X)
Z Z
P( X ≥ a ) = f ( x ) dx = f λ ( x ) dx = Eλ 1{ X ≥ a} .
x≥a x≥a f λ ( x ) fλ (X)
Usando a definição de f λ , a razão de verossimilhança é
f (X)
= e−λX +ψ(λ) ,
fλ (X)
e portanto h i h i
P( X ≥ a) = Eλ e−λX +ψ(λ) 1{ X ≥ a} = Z (λ)Eλ e−λX 1{ X ≥ a} .
Essa identidade expressa a probabilidade de um evento raro como uma esperança sob a me-
dida 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−λX por 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−λa Eλ [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 deter-
miní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:
ψ′ (λ⋆ ) = a.
86 CAPÍTULO 7. AMOSTRAGEM POR IMPORTÂNCIA
Para entender a condição para λ⋆ ótimo, calculemos a derivada da função log-partição. A
partir de Z
ψ(λ) = log eλx f ( x ) dx,
derivando em relação a λ obtemos:
R λx
′ xe f ( x ) dx
ψ (λ) = R λx .
e f ( x ) dx
Essa expressão pode ser reconhecida como a média de X sob a densidade inclinada f λ ( x ) ∝
eλx f ( x ):
ψ ′ ( λ ) = Eλ [ X ] .
Portanto, a derivada da função log-partição coincide com o valor esperado de X sob 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 proba-
bilí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.
Pela expansão de Taylor de segunda ordem de ψ, existe algum θ ∈ (0, λ) tal que
λ2 ′′ λ2 λ2
ψ ( λ ) = ψ (0) + ψ ′ (0) λ + ψ (θ ) = Varθ [ X ] ≤ sup Varθ [ X ] .
2 2 2 θ ∈(0,λ)
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 X crescem 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 qua-
drá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 7. AMOSTRAGEM POR IMPORTÂNCIA
7.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
ψ(λ) = log E[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 log-
partição verdadeira, ainda podemos controlar probabilidades de eventos raros e obter limites
exponenciais.
Por exemplo, suponha que sabemos que
σ 2 λ2
ψ(λ) ≤ , ∀λ ∈ R.
2
Isso significa que os momentos exponenciais de X crescem 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 h i h i
P( X ≥ a) = Eλ e−λX +ψ(λ) 1{ X ≥ a} = eψ(λ) Eλ e−λX 1{ X ≥ a} .
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
σ 2 λ2
P( X ≥ a) ≤ inf e−λa+ψ(λ) ≤ inf e−λa+ 2 .
λ >0 λ >0
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
( b − a )2
σ2 = .
4
Essa observação leva diretamente a um dos resultados mais fundamentais na teoria das desi-
gualdades 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] e E [ X ] = 0. Então,
para todo λ ∈ R,
h i λ2 ( b − a )2
ψ(λ) = log E eλX ≤ .
8
Em particular, X é sub-Gaussiana com parâmetro de variância σ2 = (b − a)2 /4.
7.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.
Essa desigualdade mostra que somas de variáveis aleatórias limitadas independentes exibem
2
concentração do tipo Gaussiana: suas caudas decaem tão rapidamente quanto e−ct , 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.
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
ψ∗ ( a) = sup{ λa − ψ(λ) }.
λ ∈R
Como essa desigualdade vale para qualquer λ, o lado esquerdo deve ser maior ou igual ao
maior valor possível do lado direito. Aplicando esse raciocínio, temos
Como isso vale para todo g satisfazendo Eg [ X ] = a, também vale para o menor valor possível
do lado esquerdo, isto é,
ψ∗ ( a) = sup{λa − ψ(λ)}
λ
E f XeλX
d ′
λa − ψ(λ) = a − ψ (λ) = a − .
dλ E f [eλX ]
E f λa [ X ] = a.
ψ∗ ( a) = aλ a − ψ(λ a ).
o que implica
inf DKL ( g∥ f ) ≥ ψ∗ ( a).
g:Eg [ X ]= a
ψ∗ ( a) = inf DKL ( g∥ f ).
g:Eg [ X ]= a
min DKL ( g∥ f ),
g:Eg [ X ]= a
mostrando que a inclinação exponencial não é uma escolha arbitrária, mas sim a escolha ótima
ditada pela dualidade convexa.
92 CAPÍTULO 7. AMOSTRAGEM POR IMPORTÂNCIA
Capítulo 8
n
1 1
E( X ) ≈
n
( X1 + · · · + X n ) = X n , Var( X ) ≈
n−1 ∑ ( X j − X n )2 .
j =1
A lei dos grandes números garante que essas aproximações serão boas se n for grande. Pode-
mos melhorar cada vez mais a aproximação aumentando n, bastando rodar o computador por
mais tempo (em vez de lidar com uma soma ou integral possivelmente intratável). Como dis-
cutido no Capítulo 10, essa abordagem, em que geramos valores aleatórios para aproximar uma
quantidade, é chamada de método de Monte Carlo.
Uma limitação importante da ideia de Monte Carlo é que precisamos saber como gerar as
amostras X1 , . . . , Xn (de preferência de forma eficiente, pois queremos n grande). Por exemplo,
suponha que desejamos simular de uma distribuição contínua com função densidade
Olhando apenas para a função de densidade, não é óbvio como obter uma variável aleatória com
essa distribuição. Reconhecemos, no entanto, que se trata de uma Beta(4.1, 5.2). Mesmo assim,
inverter a função de distribuição acumulada é difícil, e em distribuições mais complicadas nem
sequer conhecemos a constante de normalização da densidade.
Essas dificuldades motivam o uso de um conjunto poderoso de algoritmos que revolucio-
naram a estatística e a computação científica: os métodos de Monte Carlo via Cadeias de Markov
(MCMC). A ideia central é construir uma cadeia de Markov cuja distribuição estacionária seja
justamente a distribuição de interesse.
Antes de estudarmos o MCMC em si, precisamos entender melhor o objeto central que sus-
tenta esses métodos: as cadeias de Markov.
93
94 CAPÍTULO 8. CADEIAS DE MARKOV E MCMC
8.1 Cadeias de Markov (Resumo)
Cadeias de Markov “vivem” tanto no espaço quanto no tempo: o conjunto de valores possíveis
de Xn é chamado de espaço de estados, enquanto o índice n representa a evolução do processo ao
longo do tempo.
O espaço de estados de uma cadeia de Markov pode ser discreto ou contínuo, e o tempo
também pode ser discreto ou contínuo. Neste capítulo vamos focar exclusivamente em cadeias
de Markov com tempo discreto e espaço de estados finito. Em particular, assumiremos que Xn assume
valores em um conjunto finito, que usualmente denotamos por {1, 2, . . . , M} ou {0, 1, . . . , M }.
A condição acima é chamada de propriedade de Markov. Ela afirma que, dado o presente, o
passado e o futuro são condicionalmente independentes. Ou seja, para prever o próximo estado,
basta conhecer o estado atual.
Para descrever a dinâmica de uma cadeia de Markov, precisamos conhecer as probabilidades
de transição de qualquer estado para qualquer outro.
Definição 3 (Matriz de transição). Seja X0 , X1 , X2 , . . . uma cadeia de Markov com espaço de estados
{1, 2, . . . , M}. Definimos
qij = P( Xn+1 = j | Xn = i ),
M
como a probabilidade de transição do estado i para o estado j. A matriz Q = (qij )i,j =1 é chamada de matriz
de transição da cadeia.
A matriz de transição Q é não-negativa e cada linha soma 1, pois, dado um estado inicial i, a
cadeia deve transitar para algum estado do espaço.
Exemplo 25 (Cadeia chuva-sol). Suponha que em cada dia o clima possa ser chuvoso (R) ou ensolarado
(S). Se hoje está chuvoso, amanhã estará chuvoso com probabilidade 1/3 e ensolarado com probabilidade
2/3. Se hoje está ensolarado, amanhã estará chuvoso com probabilidade 1/2 e ensolarado com probabilidade
1/2.
Definindo Xn como o clima no dia n, temos uma cadeia de Markov com espaço de estados { R, S} e
matriz de transição !
1/3 2/3
Q= ,
1/2 1/2
onde a primeira linha corresponde ao estado R e a segunda ao estado S.
Uma vez conhecida a matriz de transição Q de uma cadeia de Markov, podemos calcular as
probabilidades de transição em horizontes de tempo maiores que um passo.
Exemplo 26 (Matriz de transição de uma cadeia de Markov com 4 estados). Considere a cadeia
de Markov com 4 estados representada na Figura abaixo. Quando não há probabilidades escritas sobre as
setas, isso significa que todas as transições saindo de um mesmo estado são igualmente prováveis.
Por exemplo, existem 3 setas saindo do estado 1, de modo que as transições 1 → 3, 1 → 2 e 1 → 1
ocorrem cada uma com probabilidade 1/3. Portanto, a matriz de transição da cadeia é
1/3 1/3 1/3 0
0 0 1/2 1/2
Q= .
0 1 0 0
1/2 0 0 1/2
Para calcular a probabilidade de que a cadeia esteja no estado 3 após 5 passos, partindo do estado 1,
basta olhar para o elemento (1, 3) da matriz Q5 .
Usando um computador, obtemos
853/3888 509/1944 52/243 395/1296
173/864 85/432 31/108 91/288
Q5 = 37/144
.
29/72 1/9 11/48
499/2592 395/1296 71/324 245/864
96 CAPÍTULO 8. CADEIAS DE MARKOV E MCMC
Assim,
(5) 52
q13 = .
243
A matriz de transição Q codifica a distribuição condicional de X1 dado o estado inicial da
cadeia. Especificamente, a i-ésima linha de Q é a PMF condicional de X1 dado X0 = i, escrita
como um vetor linha. De forma análoga, a i-ésima linha de Qn corresponde à PMF condicional
de Xn dado X0 = i.
Para obter as distribuições marginais de X0 , X1 , . . ., precisamos especificar não apenas a ma-
triz de transição, mas também as condições iniciais da cadeia. O estado inicial X0 pode ser
especificado de forma determinística, ou de forma aleatória segundo alguma distribuição. Seja
t = (t1 , t2 , . . . , t M ) a PMF de X0 , vista como vetor linha, em que ti = P( X0 = i ).
tQn .
P( Xn = j ).
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 no-
vamente a cadeia de Markov com 4 estados representada na Figura anterior. Suponha que as condições
iniciais sejam dadas por
t = 41 , 14 , 14 , 14 ,
isto é, a cadeia começa com igual probabilidade em cada um dos quatro estados.
Seja Xn a posição da cadeia no tempo n. A distribuição marginal de X1 é
1/3 1/3 1/3 0
0 0 1/2 1/2
= 5 , 1, 5 , 1 .
1 1 1 1
tQ = 4 4 4 4 0 24 3 24 4
1 0 0
0 1/2 0 0 1/2
A distribuição marginal de X5 é
853/3888 509/1944 52/243 395/1296
173/864 85/432 31/108 91/288
5 1 1 1 1 3379 2267 101 1469
tQ = 4 4 4 4 =
15552 , 7776 , 486 , 5184 .
37/144 29/72 1/9 11/48
499/2592 395/1296 71/324 245/864
Neste caso utilizamos o computador para realizar as multiplicações de matrizes.
8.1. CADEIAS DE MARKOV (RESUMO) 97
8.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 tran-
sientes, 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.
Definição 5 (Estados recorrentes e transientes). Um estado i de uma cadeia de Markov é dito recor-
rente 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 probabi-
lidade 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 i para sempre, a cadeia inevitavelmente deixará o estado i em algum
momento.
Além disso, é possível caracterizar a distribuição do número de retornos ao estado.
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 recor-
rentes 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 > 0 tal que a entrada (i, j) de Qn é
positiva.
Uma cadeia que não é irredutível é chamada de redutível.
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
(n)
i. Pela definição de irredutibilidade, existe algum n tal que q1i > 0.
Assim, toda vez que a cadeia visita o estado 1, há uma probabilidade positiva de que, após
n passos, 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 i infinitas vezes.
Como i foi 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.
8.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.
Na cadeia com 6 estados (à direita), partindo do estado 1, é possível retornar a ele após 3 passos, 6
passos, 9 passos e assim por diante. No entanto, não é possível retornar ao estado 1 em um número de
passos que não seja múltiplo de 3. Portanto, o estado 1 tem período 3. De forma análoga, os estados 2 e 3
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 8. CADEIAS DE MARKOV E MCMC
IMPORTANTE!!! Vale destacar uma diferença importante: algumas propriedades pertencem
a estados 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.
Esse sistema de equações lineares pode ser escrito de forma compacta como
sQ = s.
8.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 = s significa que, se X0 tem distribuição s, então X1 també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 s permanecerá 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 !
1 2
Q= 3 3 .
1 1
2 2
A distribuição estacionária é da forma s = (s, 1 − s). Devemos então resolver
!
1 2
(s, 1 − s) 3 3 = (s, 1 − s),
1 1
2 2
o que é equivalente ao sistema
13 s + 21 (1 − s) = s,
2 s + 1 (1 − s) = 1 − s.
3 2
as1 + (1 − b)s2 = s2 .
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.
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 distri-
buiçã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 Xn tende à distribuição estacionária s. O resultado a seguir formaliza
essa ideia.
P( Xn = i ) −→ si quando n → ∞.
Em termos matriciais, Qn converge para uma matriz cujas linhas são todas iguais a s.
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.
0 1 0 0 0
0 0 1 0 0
Q=
0 0 0 1 0.
0 0 0 0 1
1 0 0 0 0
8.1. CADEIAS DE MARKOV (RESUMO) 103
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.
Exemplo 31 (Comportamento de longo prazo de uma cadeia com dois estados). Considere nova-
mente a cadeia de dois estados discutida anteriormente, cuja matriz de transição é
!
1 2
Q= 3 3 .
1 1
2 2
A distribuição estacionária é
3 4
s= , .
7 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:
!n !
1 2 3 4
Qn = 3
1
3
1
−→ 7
3
7
4
, quando n → ∞.
2 2 7 7
8.1.3 Reversibilidade
Vimos que a distribuição estacionária de uma cadeia de Markov é extremamente útil para com-
preender seu comportamento de longo prazo. Entretanto, em muitos casos pode ser compu-
tacionalmente 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 , . . . , s M ), com si ≥ 0 e ∑i si = 1, se
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 i e ir para o estado j em um passo é si qij , e essa probabilidade
é exatamente igual à de sair de j e voltar para i, que é s j q ji . Em outras palavras, o fluxo de
probabilidade de i para j é o mesmo que o de j para i:
si qij = s j q ji .
8.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.
O lado esquerdo representa o número médio de partículas que sairão do estado j no 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 i e j,
nsi qij = ns j q ji .
O lado esquerdo corresponde ao número médio de partículas que vão de i para j, e o lado direito
ao número médio que vai de j para i. Assim, a reversibilidade garante que, par a par, o fluxo
entre dois estados é perfeitamente equilibrado.
Demonstração. Temos
∑ si qij = ∑ s j q ji = s j ∑ q ji = s j ,
i i i
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.
1 1 1
s= M, M, . . . , M .
106 CAPÍTULO 8. CADEIAS DE MARKOV E MCMC
De fato, se qij = q ji , a condição de reversibilidade si qij = s j q ji é satisfeita sempre que si = s j para
todos os pares (i, j).
Esse é um caso particular de um fato mais geral: quando cada coluna de Q também soma 1,
a distribuição uniforme continua sendo estacionária.
1 1 1
s= M, M, . . . , M ,
1 1 1 1 1 1
M, M, . . . , M Q= M, M, . . . , M ,
Uma matriz cujas linhas e colunas somam 1 é chamada de matriz duplamente estocástica. Toda
cadeia de Markov cuja matriz de transição é duplamente estocástica possui distribuição estacio-
nária uniforme.
Exemplo 32 (Caminhada aleatória em uma rede não direcionada). Uma rede é um conjunto de nós
conectados por arestas. A rede é dita não direcionada se as arestas puderem ser percorridas em ambos os
sentidos, isto é, se não houver “ruas de mão única”.
Considere um caminhante que percorre aleatoriamente as arestas de uma rede não direcionada. A partir
de um nó i, o caminhante escolhe uma das arestas conectadas a i com probabilidades iguais e então atravessa
a aresta escolhida.
O grau de um nó é o número de arestas conectadas a ele. A sequência de graus de uma rede com nós
1, 2, . . . , n é o vetor
d = ( d1 , d2 , . . . , d n ),
8.2. MARKOV CHAIN MONTE CARLO 107
onde d j é o grau do nó j. Arestas que ligam um nó a ele mesmo (self-loops) são permitidas e contam como
1 no grau desse nó.
Por exemplo, para a rede acima, a sequência de graus é
d = (4, 3, 2, 3, 2).
di qij = d j q ji ,
pois qij = 1/di se {i, j} é uma aresta e qij = 0 caso contrário (para i ̸= j). Logo, pela proposição anterior,
a distribuição estacionária é proporcional à sequência de graus:
si ∝ di .
De forma intuitiva, os nós com maior grau são mais “bem conectados”, e portanto o caminhante passa
mais tempo neles no longo prazo. No exemplo acima, isso resulta em
4 3 2 3 2
s = 14 , 14 , 14 , 14 , 14 ,
Pelo teorema da lei dos grandes números, essas aproximações se tornam mais precisas à medida
que n cresce. Podemos obter resultados cada vez melhores apenas aumentando o tempo de
simulação, sem precisar lidar com integrais ou somas intratáveis. Esse tipo de abordagem é
conhecido como método de Monte Carlo.
Um dos principais desafios dos métodos de Monte Carlo é a necessidade de saber gerar
amostras X1 , X2 , . . . , Xn da distribuição desejada. Em muitos casos, isso não é nada trivial. Por
exemplo, suponha que queremos gerar valores de uma variável contínua com densidade
O objetivo desta seção é introduzir o método Markov Chain Monte Carlo (MCMC), uma família
de algoritmos que permite gerar amostras de distribuições complexas a partir de cadeias de
Markov. O desenvolvimento do MCMC revolucionou a estatística e a computação científica, pois
tornou possível simular distribuições de alta dimensão ou com normalizadores desconhecidos.
A ideia central é inverter o problema estudado anteriormente. Antes, conhecíamos a matriz
de transição Q e procurávamos sua distribuição estacionária s. Agora, o processo é o oposto:
partimos de uma distribuição s que desejamos simular e construímos uma cadeia de Markov
cuja distribuição estacionária é exatamente s.
Se essa cadeia for executada por tempo suficiente, a distribuição de seus estados se aproxi-
mará da distribuição alvo s. O aspecto surpreendente é que isso pode ser feito sem conhecer
o termo de normalização de s, o que torna o método aplicável em uma enorme variedade de
contextos.
s j p ji
aij = min ,1 .
si pij
5. Com probabilidade aij , aceite a proposta e defina Xn+1 = j; caso contrário, rejeite a proposta
e mantenha Xn+1 = i.
A matriz P é chamada de matriz de proposta, pois serve apenas para gerar possíveis movimen-
tos da cadeia. As probabilidades de aceitação aij ajustam essas propostas para que o equilíbrio
da nova cadeia Q obedeça à distribuição s. Em notação compacta, a nova matriz de transição é
pij aij ,
i ̸= j,
qij =
1 − ∑ p a , i = j.
k ̸=i ik ik
• se pii > 0, pode ocorrer de a proposta coincidir com o estado atual, e a cadeia naturalmente
permanece em i.
O caso i = j é imediato, pois ambos os lados são iguais a si qii . Consideremos agora i ̸= j.
Se qij > 0, então pij > 0 — a cadeia só pode propor uma transição possível — e também p ji >
0, pois, do contrário, a probabilidade de aceitação seria nula. De forma análoga, se pij , p ji > 0,
então q ji > 0. Portanto, qij e q ji são simultaneamente nulos ou não nulos.
110 CAPÍTULO 8. CADEIAS DE MARKOV E MCMC
Para i ̸= j com qij > 0, temos
qij = pij aij ,
pois, partindo de i, a única forma de alcançar j é propor essa transição (com probabilidade pij ) e
aceitá-la (com probabilidade aij ).
Caso 2: se s j p ji > si pij , o mesmo raciocínio vale trocando i e j, de modo que novamente
si qij = s j q ji .
ns j p ji a ji .
A ideia central do algoritmo é ajustar as probabilidades de aceitação aij de modo que, no regime estaci-
onário, esses dois fluxos se equilibrem:
ou, equivalentemente,
si qij = s j q ji .
• Se o fluxo proposto é menor, isto é, há poucas partículas saindo de i em relação às que retornam
de j, o sistema já tem menos movimento em direção a j. Nesse caso, todas as propostas são aceitas
(aij = 1), pois não há risco de desequilíbrio: permitir todas as transições ajuda a compensar o déficit
de fluxo, mantendo o equilíbrio entre os dois estados.
Com esse ajuste, o sistema tende a um equilíbrio dinâmico, no qual a taxa média de partículas indo
de i para j é igual à taxa de partículas voltando de j para i. Não há acúmulo nem escoamento líquido de
probabilidade entre os estados, e a igualdade
si qij = s j q ji
u a −1 (1 − u ) b −1
a(w, u) = min ,1 .
w a −1 (1 − w ) b −1
4. Com probabilidade a(w, u), aceite a proposta e defina Wn+1 = u; caso contrário, rejeite a proposta e
mantenha Wn+1 = w.
Observe que o algoritmo não requer o conhecimento da constante de normalização da densidade Beta( a, b),
pois ela se cancela na razão entre as densidades no numerador e no denominador. Aqui, a densidade
Beta( a, b) desempenha o papel de s (a distribuição estacionária desejada), enquanto a densidade uniforme
Unif(0, 1) desempenha o papel de pij e p ji , já que as propostas são sempre independentes do estado atual.
112 CAPÍTULO 8. CADEIAS DE MARKOV E MCMC
s(u) p( x | u)
a( x, u) = min 1, ,
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.
p( x |u)
O termo p(u| x) está presente para corrigir possíveis assimetrias da distribuição de proposta. Por
p( x |u)
exemplo, se é mais fácil propor de x para u do que o contrário, o fator p(u| x) compensa essa diferença,
garantindo que o fluxo médio de partículas entre x e u permaneça equilibrado:
Consequentemente,
p( x | u)
p( x | u) = p( x ) = Unif(0, 1), e portanto = 1.
p(u | x )
O fator de correção se cancela, e a fórmula de aceitação se reduz para
s(u)
a( x, u) = min 1, ,
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 indepen-
dente 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.
8.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
Xn ainda 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.
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.
Desejamos construir uma cadeia de Markov ( Xn , Yn ) cuja distribuição estacionária seja p X,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:
x n +1 ∼ p ( x | Y = y n ).
Na versão aleatória, escolhe-se a cada iteração qual variável será atualizada, com probabili-
dades iguais. O procedimento segue:
3. Se X for escolhido:
Repetindo o processo, obtemos novamente uma cadeia cuja distribuição estacionária é p X,Y .
O algoritmo de Gibbs se estende naturalmente para o caso de d variáveis aleatórias. Nesse
(1) (d)
caso, o estado da cadeia é um vetor Wn = (Wn , . . . , Wn ). Em cada iteração:
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.
p( x ′ , y′ ) q(( x, y) | ( x ′ , y′ ))
′ ′
a(( x, y), ( x , y )) = min 1, .
p( x, y) q(( x ′ , y′ ) | ( x, y))
No caso do Gibbs aleatório, a proposta consiste em escolher aleatoriamente uma das coorde-
nadas e atualizá-la de acordo com sua distribuição condicional verdadeira. Mais precisamente:
Vamos considerar o segundo caso, em que apenas Y é atualizado (o caso simétrico em que X
é atualizado é análogo). Assim,
1 1
q(( x, y′ ) | ( x, y)) = p(y′ | x ) e q(( x, y) | ( x, y′ )) = p ( y | x ).
2 2
Substituindo esses termos na fórmula de aceitação, obtemos:
p( x, y′ ) p(y | x )
′
a(( x, y), ( x, y )) = min 1, .
p( x, y) p(y′ | x )
Logo,
a(( x, y), ( x, y′ )) = 1.
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] ,
X | p ∼ Poisson(λp).
Essa distribuição não tem forma fechada conhecida, o que dificulta a amostragem direta. Para con-
tornar isso, introduzimos a variável latente N, correspondente ao número total de ovos postos.
X | N = n, p ∼ Binomial(n, p).
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.
p | X = x, N = n ∼ Beta( x + a, n − x + b).
Y | p ∼ Poisson(λ(1 − p)).
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:
T
1
E( p | X = x ) ≈
T ∑ p(t) .
t =1
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
8.4. AMOSTRAGEM DE GIBBS 117
considerados diferentes valores observados de ovos chocados X ∈ {3, 5, 7, 9}, enquanto os demais
parâmetros permaneceram fixos.
X E [ 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 8. CADEIAS DE MARKOV E MCMC
Capítulo 9
Processos de Difusão
Var(Wt+s − Wt ) = sId .
Essa relação mostra que quanto maior o intervalo s, maior a variabilidade esperada do incre-
mento, refletindo a natureza difusiva do processo.
119
120 CAPÍTULO 9. 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.
3. Geração dos incrementos: Para cada passo k = 1, 2, . . . , n, gere um ruído gaussiano inde-
pendente
ε k ∼ N (0, 1),
e compute o incremento
√
∆Wk = s εk.
Wk = Wk−1 + ∆Wk .
O resultado é mostrado na Figura 9.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 9.1: Simulação de duas trajetórias de movimento browniano com diferentes tamanhos de
passo s.
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 Xt segundo
dXt = f ( Xt , t) dt + g( Xt , t) dWt .
9.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 Bt possui incrementos
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 ,
Xt+∆ ≈ Xt + f ( Xt , t) ∆ + g( Xt , t) ( Bt+∆ − Bt ).
Essa equação descreve como o processo Xt evolui 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 es-
tocá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 σ > 0 controla a intensidade do ruído.
f1 (x) = −x e f 2 ( x ) = x − x3 .
(b) Gere várias trajetórias para cada caso, mantendo o mesmo valor de σ e do passo temporal ∆.
(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.
122 CAPÍTULO 9. PROCESSOS DE DIFUSÃO
9.2 Equação de Fokker–Planck
dXt = f ( Xt , t) dt + g( Xt , t) dBt ,
onde f é o campo de drift e g o coeficiente de difusão. A função g controla como o ruído atua sobre
cada componente do sistema — por exemplo, se g é uma matriz, o ruído pode ter intensidades e
correlações diferentes em cada direção.
A densidade de probabilidade p( x, t) associada a Xt satisfaz a equação de Fokker–Planck
(também chamada de equação de Kolmogorov para frente):
∂p 1
= −∇ · f ( x, t) p( x, t) + ∇ · ∇ · ( g g⊤ )( x, t) p( x, t) .
∂t 2
O primeiro termo representa o transporte determinístico da densidade pelo campo de drift, en-
quanto o segundo termo descreve a difusão espacial causada pelo ruído multiplicativo g( x, t) dBt .
√
No caso em que o ruído é isotrópico e constante, isto é, g( x, t) = 2D Id , a equação se simplifica
para
∂p
= −∇ · f ( x ) p( x, t) + D ∆p( x, t),
∂t
onde D > 0 é o coeficiente de difusão e ∆ é o operador Laplaciano.
Essa forma mostra claramente a dualidade entre drift e difusão:
A equação de Fokker–Planck pode ainda ser escrita como uma equação de continuidade:
∂p
= −∇ · J, J = f ( x ) p( x, t) − D ∇ p( x, t),
∂t
onde J é o fluxo de probabilidade. O primeiro termo de J corresponde ao transporte devido ao drift,
e o segundo termo é o fluxo difusivo induzido pela variação espacial da densidade.
Aqui, o operador ∇ · J representa o divergente do campo de fluxo J, isto é, a taxa líquida de
probabilidade que sai (ou entra) de uma pequena região do espaço. De forma intuitiva:
Essa SDE é conhecida como dinâmica de Langevin (ou overdamped Langevin equation). O termo
determinístico ∇ log p( Xt ) empurra as partículas em direção às regiões de maior probabilidade
√
da distribuição, enquanto o ruído gaussiano 2 dBt garante que o processo explore todo o espaço
de estados.
Se tomarmos p∞ ( x ) = p( x ), queremos verificar que essa escolha satisfaz a condição estacio-
nária da equação de Fokker–Planck:
∇ · ∇ log p( x ) p( x ) = ∆p( x ).
∆p( x ) = ∇ · ∇ p( x ) .
d2 p
Em uma dimensão, o Laplaciano reduz-se simplesmente à segunda derivada dx2
.
∇ p( x )
∇ log p( x ) = .
p( x )
∇ p( x ) = p( x ) ∇ log p( x ).
∆p( x ) = ∇ · ∇ p( x ) = ∇ · p( x ) ∇ log p( x ) .
Essa relação entre a SDE de Langevin e a Fokker–Planck mostra que o processo preserva a
densidade p( x ) no equilíbrio. Na prática, a amostragem de Langevin consiste em discretizar
essa SDE (via o método de Euler–Maruyama) e usar as iterações resultantes para gerar amostras
aproximadamente distribuídas segundo p( x ).
124 CAPÍTULO 9. PROCESSOS DE DIFUSÃO
Exercício 29. Considere a densidade alvo p( x ) = N ( x; 0, 1), isto é,
1 2
p( x ) = √ e− x /2 .
2π
∇ log p( x ) = − x.
Observe que o termo constante − 12 log(2π ) desaparece na derivada. Discuta o motivo pelo qual o
termo normalizador da densidade não afeta o gradiente.
Use ∆ = 0.01 e inicialize X0 ∼ N (0, 1). Gere trajetórias e verifique que o histograma das amostras
converge para a densidade p( x ).
Xt+∆ = Xt − Xt ∆.
O processo converge para qual valor? Explique o papel do ruído na manutenção da variabilidade da
amostra.
4. Interpretação: Analise a dinâmica obtida. O termo − Xt ∆ atua como uma força de retorno à origem
√
(um poço de potencial), enquanto o termo aleatório 2∆ ε t impede que o processo colapse nesse ponto,
mantendo a distribuição estacionária N (0, 1).
Note que o termo normalizador √12π não influencia o gradiente ∇ log p( x ), pois é constante. Assim,
o comportamento da dinâmica de Langevin depende apenas da forma relativa de p( x ), e não da sua escala
global. Esse fato é fundamental: em métodos baseados em gradientes (como Langevin e Hamiltonian Monte
Carlo), apenas a derivada do logaritmo da densidade importa.
Exercício 30. Considere o problema de gerar amostras de uma mistura bimodal de normais:
p( x ) = 1
2 N ( x; −3, 12 ) + 12 N ( x; 3, 12 ).
1. Derive o gradiente de log p( x ). Expresse ∇ log p( x ) em termos das densidades de cada compo-
nente da mistura e de suas médias.
9.4. DENOISING SCORE MATCHING 125
2. Implemente a dinâmica de Langevin:
√
Xt+∆ = Xt + ∇ log p( Xt ) ∆ + 2∆ ε t , ε t ∼ N (0, 1).
Use ∆ = 0.02 e inicialize X0 ∼ N (0, 5). Gere trajetórias e compare o histograma das amostras
obtidas com a densidade-alvo p( x ).
Xt+∆ = Xt + ∇ log p( Xt ) ∆,
e observe o comportamento do processo. Ele converge para algum ponto específico? O que muda em
relação à versão estocástica?
p( x ) = 1
2 N ( x; −3, 12 ) + 12 N ( x; 3, 12 ),
e a dinâmica de Langevin
√
dXt = ∇ log p( Xt ) dt + 2 dBt .
1. Simule uma trajetória longa da dinâmica, iniciando em X0 ∼ N (0, 5). Plote o histograma das
amostras e o gráfico de posição t 7→ Xt .
Discuta o que está acontecendo no processo e como isso se liga com a matéria de MCMC.
2. Repita o experimento com várias trajetórias independentes cada uma com condição inicial dife-
rente, e compare os histogramas dos dois métodos.
possui como distribuição estacionária a própria densidade alvo p( x ). De fato, pela equação de
Fokker–Planck associada,
∂pt ( x ) 1
log p( x ) pt ( x ) + 12 ∆pt ( x ),
= −∇ x · 2 ∇x
∂t
vemos que pt ( x ) = p( x ) é uma solução estacionária.
O problema, entretanto, é que na prática não conhecemos p( x ) de forma explícita e portanto não
temos acesso ao seu gradiente ∇ x log p( x ), o score, que aparece diretamente no termo de drift da
dinâmica de Langevin. O que dispomos, em geral, são apenas amostras independentes
x1 , . . . , x n ∼ p ( x ),
126 CAPÍTULO 9. PROCESSOS DE DIFUSÃO
e queremos, a partir delas, construir um estimador para o score ou, de forma equivalente, um
campo vetorial sθ ( x ) ≈ ∇ x log p( x ) que possa ser usado em dinâmicas como a de Langevin para
simular a distribuição alvo.
O objetivo do Denoising Score Matching (DSM) é aprender o score de uma distribuição de
probabilidade desconhecida p( x ), definido por
s⋆ ( x ) = ∇ x log p( x ).
A ideia central do DSM é transformar o problema de estimar o score em uma tarefa de regres-
são supervisionada. Para isso, partimos de amostras x ∼ p( x ) e adicionamos ruído gaussiano
ε ∼ N (0, σ2 ), obtendo dados corrompidos
x̃ = x + ε.
x − x̃
t= ,
σ2
e treina um modelo sθ ( x̃, σ) para aproximar esse alvo a partir do dado corrompido. A função de
perda considerada é h i
2
L(θ ) = Ex∼ p Eε∼N (0,σ2 ) sθ ( x̃, σ) − t .
Para entender por que esse procedimento funciona, começamos observando que a variável
corrompida x̃ possui densidade
Z
qσ ( x̃ ) = p( x ) φσ ( x̃ − x ) dx,
e como
x̃ − x
∇ x̃ φσ ( x̃ − x ) = − φσ ( x̃ − x ),
σ2
temos
x − x̃
Z
∇ x̃ qσ ( x̃ ) = p( x ) φσ ( x̃ − x ) dx.
σ2
Dividindo por qσ ( x̃ ) e lembrando que a densidade conjunta de ( x, x̃ ) é dada por p( x, x̃ ) =
p( x ) φσ ( x̃ − x ), podemos reescrever o integrando em termos da distribuição condicional de x
dado x̃:
p( x ) φσ ( x̃ − x )
p( x | x̃ ) = .
qσ ( x̃ )
9.4. DENOISING SCORE MATCHING 127
Substituindo essa expressão na fórmula do gradiente, obtemos
x − x̃
Z
∇ x̃ log qσ ( x̃ ) = p( x | x̃ ) dx,
σ2
o que mostra que o score de qσ é precisamente a esperança condicional
x − x̃
∇ x̃ log qσ ( x̃ ) = E x̃ .
σ2
s⋆ ( x̃, σ) = ∇ x̃ log qσ ( x̃ ).
Note que, fixado o valor de σ, o DSM está aprendendo o score não da densidade original
p( x ), mas da versão suavizada qσ ( x ) = ( p ∗ φσ )( x ), obtida ao convoluir p com uma gaussiana de
variância σ2 . Ou seja, estamos removendo detalhes finos da distribuição original e focando nas
variações mais suaves de p.
O campo vetorial ∇ x log qσ ( x ) aponta na direção em que a densidade suavizada mais cresce:
ele indica como deveríamos mover uma partícula que foi corrompida pelo ruído para restaurar
a estrutura de p.
Assim, quando treinamos sθ ( x, σ) para aproximar E[( x − x̃ )/σ2 | x̃ = x ], estamos aprendendo
como “desfazer” o ruído gaussiano de nível σ. Por isso o nome denoising score matching.
Na dedução anterior vimos que, para um valor fixo de σ, o modelo sθ ( x̃, σ) treinado com o
Denoising Score Matching aprende o score da densidade suavizada qσ ( x ). Na prática, entretanto,
não queremos apenas um único valor de σ, mas uma coleção de níveis de ruído que permitam
capturar a estrutura de p( x ) em diferentes escalas.
Usualmente escolhemos uma escala geométrica de valores
k/(K−1)
σmax
σk = σmin , k = 0, . . . , K − 1,
σmin
128 CAPÍTULO 9. 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 xi e para cada σk , podemos gerar várias
versões corrompidas
x̃i,j,k = xi + ε i,j,k , ε i,j,k ∼ N (0, σk2 I ),
Dessa forma, podemos expandir o conjunto de dados criando várias amostras supervisionadas
( x̃i,j,k , σk , ti,j,k ), que descrevem, para diferentes níveis de ruído, a direção de denoising a ser apren-
dida.
O próximo passo é ajustar um modelo de regressão sθ que receba como entrada o ponto
corrompido x̃ e 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̃, σ) ≈ ∇ x̃ log 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 ).
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 xi do conjunto de dados e para cada nível de ruído σk , gerar algumas
versões corrompidas x̃i,j,k = xi + ε i,j,k , com ε i,j,k ∼ N (0, σk2 I ). Isso aumenta o tamanho do
conjunto de treinamento e ajuda o modelo a aprender a remover diferentes intensidades de
ruído.
9.4. DENOISING SCORE MATCHING 129
3. Calcular o alvo de regressão para cada par ( x̃i,j,k , σk ) como ti,j,k = −ε i,j,k /σk2 , que representa a
direção na qual a amostra corrompida deve ser movida para retornar à distribuição original.
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 t a partir de [ 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.
3. Para cada valor de σ, corrompa as amostras adicionando ruído gaussiano ε ∼ N (0, σ2 I ), 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, compa-
rando visualmente os resultados de diferentes modelos.
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
αi √
x t = x t −1 + sθ ( xt−1 , σi ) + αi zt , zt ∼ N (0, I ),
2
onde cada nível de ruído σi controla a escala das atualizações e αi é o passo de integração
proporcional a σi2 . 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 9. PROCESSOS DE DIFUSÃO
1. Inicialização: Inicie x0 como uma amostra de uma distribuição de ruído, por exemplo
2 I ).
x0 ∼ N (0, σmax
escala.
• Repita T vezes (por exemplo, T = 10):
αi √
x ← x+ 2 sθ ( x, σi ) + αi z, z ∼ N (0, I ).
3. Saída: Após percorrer todos os níveis de ruído, o vetor final x T é uma amostra aproximada
da distribuição de interesse p( x ).
Figura 9.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.
9.4. DENOISING SCORE MATCHING 131
1. Gere o conjunto de dados bidimensional usando a função make_moons da biblioteca [Link].
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.
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 caracte-
rí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 condici-
onal, 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 ∇ x̃ log 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
∇y log 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 yi com ruído gaussiano ε i,j,k ∼ N (0, σk2 I ), o alvo passa a ser
ti,j,k = −ε i,j,k /σk2 , e o modelo deve receber [ỹi,j,k , xi , σk ] como entrada.
4. Adapte o algoritmo de Annealed Langevin Dynamics para o caso condicional, mantendo x fixo e
atualizando apenas y:
αi √
yt = yt−1 + sθ (yt−1 , x, σi ) + αi zt , zt ∼ N (0, I ).
2
132 CAPÍTULO 9. PROCESSOS DE DIFUSÃO
5. Escolha um conjunto de dados simples para testar o modelo, por exemplo:
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 10
Bootstrap
1 n
n i∑
F̂n ( x ) = 1 { Xi ≤ x } .
=1
O bootstrap consiste em gerar amostras X1∗ , . . . , Xn∗ i.i.d. de F̂n (isto é, reamostrar com reposição
dos dados observados), e então computar
θ̂ ∗ = t( F̂n∗ ).
1 n
n i∑
θ̂ = X̄ = Xi .
=1
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 10. 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 X̄ sem 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 X1∗ , . . . , Xn∗ . Cada
amostra reamostrada define uma distribuição empírica F̂n∗ .
1 n ∗
n i∑
X̄ ∗ = Xi .
=1
O conjunto dessas médias forma uma aproximação empírica da distribuição de X̄. O desvio padrão das
médias reamostradas,
v
B B
u
u 1
¯∗ = 1
B − 1 b∑ ∑ X̄ ∗(b) ,
σ̂boot = t X̄ ¯ ∗ 2 ,
∗(b) − X̄ X̄
=1
B b =1
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.
10.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̂n−1 (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 F̂n (θ̂ ) = 1/2,
podemos relacionar θ̂ e θ por uma expansão local de F em torno de θ:
F (θ̂ ) ≈ F (θ ) + f (θ )(θ̂ − θ ),
0 ≈ F̂n (θ ) − F (θ ) + f (θ )(θ̂ − θ ).
√
Multiplicando por n,
√ 1 √
n(θ̂ − θ ) ≈ − n F̂n (θ ) − F (θ ) .
f (θ )
Pelo Teorema Central do Limite empírico,
√
n F̂n (θ ) − F (θ ) ⇒ N 0, F (θ )(1 − F (θ )) = N (0, 1/4).
Portanto,
1
Var θ̂ ≈ .
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) ∗(b)
θ̂ ∗(b) = median( X1 , . . . , Xn ), b = 1, . . . , B.
136 CAPÍTULO 10. BOOTSTRAP
A variabilidade entre as medianas reamostradas fornece uma estimativa do erro padrão de θ̂,
v
B
1 B ∗(b)
u
u 1
B − 1 b∑ B b∑
∗(b) − θ̄ ∗ 2 , ∗
σ̂boot = t θ̂ θ̄ = θ̂ .
=1 =1
Dessa forma, mesmo sem conhecer f (θ ) nem a forma de F, é possível avaliar empiricamente a incerteza da
mediana.
Antes de entender por que o bootstrap funciona, é preciso quantificar o quão bem a distribuição
empírica F̂n aproxima 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
2
P sup F̂n ( x ) − F ( x ) > ε ≤ Ce−nε , ∀ ε > 0.
x
Essa desigualdade mostra que o erro uniforme entre F e F̂n decai exponencialmente com n.
Em particular, F̂n converge quase certamente para F, o que é o conteúdo do teorema de Gli-
venko–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 F̂n está 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, tam-
bé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.
10.2. UMA VISÃO TEÓRICA DE BOOTSTRAP 137
Teorema 12 (Desigualdade de McDiarmid). Sejam X1 , . . . , Xn variáveis independentes assumindo
valores em um conjunto arbitrário X , e seja f : X n → R uma função tal que
| f ( x1 , . . . , xi , . . . , xn ) − f ( x1 , . . . , xi′ , . . . , xn )| ≤ ci
para todo i e para todos os valores possíveis das variáveis. Então, para todo ε > 0,
2ε2
P ( f ( X1 , . . . , Xn ) − E [ f ( X1 , . . . , Xn )] ≥ ε) ≤ exp − n 2 .
∑ i =1 c 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
Observe que alterar um único Xi muda no máximo um termo da soma que define F̂n ( x ), e
portanto o valor de f só pode variar em 1/n. Assim, podemos tomar ci = 1/n para todo i, o que
dá
n
1 1
∑ c2i = n × n2 = n .
i =1
Substituindo isso na desigualdade de McDiarmid, obtemos
2
P sup | F̂n ( x ) − F ( x )| − E sup | F̂n ( x ) − F ( x )| ≥ ε ≤ e−2nε .
x x
Esse é o passo essencial da prova da DKW: ele mostra que a distância uniforme entre F̂n e F
está fortemente concentrada em torno de sua média. O passo seguinte consiste em controlar o
valor esperado E supx | F̂n ( x ) − F ( x )| .
Para controlar o valor esperado de supx | F̂n ( 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.
138 CAPÍTULO 10. BOOTSTRAP
Lema 1. Sejam Z1 , . . . , Zk variáveis aleatórias independentes, centradas e subgaussianas com parâmetro
σ2 , ou seja, h i 2 2
E etZj ≤ et σ /2 , ∀t ∈ R, j = 1, . . . , k.
Então
E max Zj ≤ σ
p
2 log k.
1≤ j ≤ k
2 2
Pela hipótese subgaussiana, E etZj ≤ et σ /2 para todo j, logo
log k tσ2
1 2 2
E max Zj ≤ log ket σ /2 = + .
j t t 2
E max Zj ≤ σ 2 log k.
p
1≤ j ≤ k
1 n
∑
F̂n ( x ) − F ( x ) = 1 { Xi ≤ x } − F ( x )
n i =1
é a média de n variáveis independentes, centradas e limitadas no intervalo [−1, 1]. Portanto, cada
uma é subgaussiana com parâmetro 1/n, o que implica que o vetor ( F̂n ( x1 ) − F ( x1 ), . . . , F̂n ( xk ) −
√
F ( xk )) é formado por variáveis subgaussianas com desvio σ = 1/ n.
Com isso, temos todos os elementos da prova. A desigualdade de McDiarmid garante que
2
P sup | F̂n ( x ) − F ( x )| − E sup | F̂n ( x ) − F ( x )| > ε ≤ e−2nε .
x x
Pelo Teorema de Donsker, temos que Gn ⇒ B0 , onde B0 (t) é a ponte browniana padrão, um processo
gaussiano centrado com covariância
Lembrando que, por construção, o processo B0 (t) é centrado, isto é, E [ B0 (t)] = 0 para todo
t ∈ [0, 1]. Tomando então s = t,
o que implica que a variância é máxima no meio do intervalo (t = 1/2) e decresce até zero nos
extremos. Isso mostra que o processo necessariamente satisfaz B0 (0) = B0 (1) = 0 — o início e o
fim são fixos, mas entre eles o comportamento é aleatório e gaussiano.
Intuitivamente, portanto, uma ponte browniana é um movimento browniano “amarrado” nas
extremidades: começa em zero, flutua livremente no interior, e é progressivamente puxado de
volta até atingir zero novamente em t = 1.
Podemos agora perguntar: que tipo de equação estocástica seria capaz de gerar um pro-
cesso com esse comportamento? Um movimento browniano comum Wt é livre: suas flutua-
ções acumulam-se ao longo do tempo, e sua variância cresce linearmente, Var [Wt ] = t. Para
transformá-lo em uma ponte, precisamos forçá-lo a retornar a zero em t = 1. Isso pode ser feito
10.2. UMA VISÃO TEÓRICA DE BOOTSTRAP 141
introduzindo um termo de drift que “corrige” a trajetória conforme o tempo avança, puxando o
processo de volta à origem.
A ponte browniana padrão pode ser descrita como a solução da equação estocástica
B0 (t)
dB0 (t) = − dt + dWt , B0 (0) = 0, t ∈ [0, 1).
1−t
O termo de ruído dWt representa o movimento browniano livre, enquanto o termo de drift
− B0 (t)/(1 − t) atua como uma força restauradora que cresce quando t → 1. Essa força é justa-
mente o que garante que o processo retorne a zero no instante final.
mostrando que o processo termina praticamente em zero, restando apenas um ruído residual
que desaparece no limite ∆t → 0.
Assim, o termo de drift − B0 (t)/(1 − t) surge naturalmente como o mecanismo que, ao longo
do tempo, ajusta continuamente as trajetórias para que todas convirjam exatamente a B0 (1) = 0.
Em termos geométricos, trata-se de um movimento browniano sujeito a uma força elástica cada
vez mais intensa, que o “puxa” para a origem conforme o tempo restante se esgota.
Seguindo a mesma lógica do caso uniforme, o processo limite do Teorema de Donsker para
uma distribuição qualquer F também deve satisfazer uma condição de fronteira análoga. De fato,
vimos que a covariância do processo limite é
f (x)
q
dGF ( x ) = − GF ( x ) dx + f ( x ) dWx .
1 − F(x)
f ( xk )
q √
Gk+1 = Gk − G ∆x + f ( xk ) ξ k ∆x, ξ k ∼ N (0, 1), G0 = 0.
1 − F ( xk ) k
Tome n = 1000, gere várias trajetórias independentes de ambos os processos e trace o envelope de 90%
ponto a ponto para comparar visualmente. O envelope pode ser obtido tomando, para cada xk , os percentis
5 e 95 das trajetórias simuladas, por exemplo em Python:
O arquivo gerado pode ser inspecionado com o módulo pstats, que permite ordenar e filtrar
resultados:
1 python -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 )
Também é possível usar cProfile dentro do código, o que facilita em notebooks ou quando
queremos medir apenas um trecho específico:
1 import cProfile , pstats , io
2
3 pr = cProfile . Profile ()
145
146 CAPÍTULO 11. ESTRATÉGIAS PARA ACELERAR CÓDIGOS EM PYTHON
4 pr . enable ()
5
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 () )
• time: tempo gasto apenas dentro da função, sem contar chamadas internas.
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 loca-
lizar 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 pro-
filer em cada uma delas, exibindo as funções mais custosas. O script pode ser salvo como
profile_matmul.py.
1 import numpy as np
2 import math
3 import cProfile , pstats , io
4 import time
5
37 A = np . random . rand (n , n )
38 B = np . random . rand (n , n )
39
Esse script pode ser executado normalmente com python profile_matmul.py. Outra forma é
rodar o profiler diretamente no terminal, usando python -m cProfile -o [Link] profile_matmul.py.
Nesse caso o resultado fica salvo em [Link], 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 consu-
mido dentro de matmul_naive:
1 7 function calls in 8.532 seconds
2
Ao comparar com a versão vetorizada, vemos que a execução termina em milésimos de se-
gundo, com o tempo todo acumulado em matmul_numpy:
148 CAPÍTULO 11. ESTRATÉGIAS PARA ACELERAR CÓDIGOS EM PYTHON
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.
Receita de bolo
No exemplo acima, o parâmetro n_jobs define quantos trabalhadores serão usados (tipica-
mente 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 coor-
dena 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
11.2. PARALELIZAÇÃO COM [Link] 149
começar. Já com vários cozinheiros (workers), cada um recebe um prato e trabalha nele indepen-
dentemente, 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,
o joblib 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 pe-
sados. 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 .
3 res = Parallel ( n_jobs =8 , backend = " threading " ) (
4 delayed ( processa_io ) ( u ) for u in urls
5 )
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:
1 from joblib import Parallel , delayed
2 import math
3
4 def eh_primo ( n ) :
5 for i in range (2 , int ( math . sqrt ( n ) ) +1) :
6 if n % i == 0:
7 return False
8 return True
9
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 11. ESTRATÉGIAS PARA ACELERAR CÓDIGOS EM PYTHON
1 import requests
2 urls = [ " https :// httpbin . org / delay /1 " ] * 20
3
Exemplo I/O-bound
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.
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 11. ESTRATÉGIAS PARA ACELERAR CÓDIGOS EM PYTHON
51 out = func (* args )
52 t1 = time . perf_counter ()
53 print ( f " { label } [1 a chamada ]: { t1 - t0 :.3 f } s " )
54
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 ) " )
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 ver-
sã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.
1 import numpy as np
2 from numba import vectorize , float64
3
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++.
3 for i in $ ( seq 1 5)
4 do
5 echo " Iniciando tarefa $i "
6 sleep 3 &
7 done
8
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
3 seria executado em sequencia.
Para visualizar essa diferenca, vejamos primeiro a execucao sequencial:
1 # !/ bin / bash
154 CAPÍTULO 11. ESTRATÉGIAS PARA ACELERAR CÓDIGOS EM PYTHON
2
3 for i in $ ( seq 1 5)
4 do
5 echo " Rodando tarefa $i "
6 sleep 3
7 done
8
Execucao sequencial
E agora a versao em paralelo, onde o tempo total cai para cerca de 3 segundos:
1 # !/ bin / bash
2
3 for i in $ ( seq 1 5)
4 do
5 echo " Rodando tarefa $i "
6 sleep 3 &
7 done
8
9 wait
10 echo " Todas as tarefas terminaram ( paralelo ) "
Execucao em paralelo
Para garantir que o script so finalize depois que todas as tarefas concluirem, podemos usar
explicitamente o comando wait:
1 # !/ bin / bash
2
3 for i in $ ( seq 1 5)
4 do
5 sleep 3 &
6 done
7
8 wait
9 echo " Todas as tarefas terminaram ! "
Tambem é possivel limitar quantas tarefas rodam em paralelo. Uma tecnica simples é contro-
lar com um contador e usar wait -n para esperar pelo menos um job terminar antes de lancar o
proximo:
1 # !/ bin / bash
2
5 for i in $ ( seq 1 5)
6 do
11.4. PARALELISMO SIMPLES EM BASH 155
7 sleep 3 &
8 if (( $ ( jobs -r | wc -l ) >= N ) ) ; then
9 wait -n
10 fi
11 done
12
13 wait
14 echo " Fim das tarefas "
Esses exemplos usam apenas comandos nativos (sleep, echo), mas a ideia é exatamente a
mesma se quisermos chamar um script Python ou outro programa no lugar.
156 CAPÍTULO 11. ESTRATÉGIAS PARA ACELERAR CÓDIGOS EM PYTHON
Referências Bibliográficas
Efron, B. (1979). Bootstrap Methods: Another Look at the Jackknife. The Annals of Statistics, 7(1):1
– 26. 133
157