Apostila Inferência Bayesiana - Ricardo Ehlers, Notas de estudo de Estatística
fernanda-ribeiro-21
fernanda-ribeiro-21

Apostila Inferência Bayesiana - Ricardo Ehlers, Notas de estudo de Estatística

106 páginas
50Números de download
1000+Número de visitas
Descrição
Apostila Inferência Bayesiana - Ricardo Ehlers
80 pontos
Pontos de download necessários para baixar
este documento
Baixar o documento
Pré-visualização3 páginas / 106
Esta é apenas uma pré-visualização
3 mostrados em 106 páginas
Esta é apenas uma pré-visualização
3 mostrados em 106 páginas
Esta é apenas uma pré-visualização
3 mostrados em 106 páginas
Esta é apenas uma pré-visualização
3 mostrados em 106 páginas
bayes.dvi

INFERÊNCIA BAYESIANA

RICARDO S. EHLERS

Primeira publicaç ao em 2002

Segunda edição publicada em 2004

Terceira edição publicada em 2005

Quarta edição publicada em 2006

Quinta edição publicada em 2007

© RICARDO SANDES EHLERS 2003-2011

Sumário

1 Introdução 1

1.1 Teorema de Bayes . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.2 Prinćıpio da Verossimilhança . . . . . . . . . . . . . . . . . . . . . 11

1.3 Exerćıcios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

2 Distribuições a Priori 14

2.1 Prioris Conjugadas . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.2 Conjugação na Famı́lia Exponencial . . . . . . . . . . . . . . . . . 15

2.3 Principais Famı́lias Conjugadas . . . . . . . . . . . . . . . . . . . 19

2.3.1 Distribuição normal com variância conhecida . . . . . . . . 19

2.3.2 Distribuição de Poisson . . . . . . . . . . . . . . . . . . . . 20

2.3.3 Distribuição multinomial . . . . . . . . . . . . . . . . . . . 21

2.3.4 Distribuição normal com média conhecida e variância de-

sconhecida . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

2.3.5 Distribuição normal com média e variância desconhecidos . 23

2.4 Priori não Informativa . . . . . . . . . . . . . . . . . . . . . . . . 25

2.5 Prioris Hierárquicas . . . . . . . . . . . . . . . . . . . . . . . . . . 28

2.6 Problemas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

3 Estimação 35

3.1 Introdução à Teoria da Decisão . . . . . . . . . . . . . . . . . . . 35

3.2 Estimadores de Bayes . . . . . . . . . . . . . . . . . . . . . . . . . 36

3.3 Estimação por Intervalos . . . . . . . . . . . . . . . . . . . . . . . 38

3.4 Estimação no Modelo Normal . . . . . . . . . . . . . . . . . . . . 39

3.4.1 Variância Conhecida . . . . . . . . . . . . . . . . . . . . . 40

3.4.2 Média e Variância desconhecidas . . . . . . . . . . . . . . 41

3.4.3 O Caso de duas Amostras . . . . . . . . . . . . . . . . . . 42

3.4.4 Variâncias desiguais . . . . . . . . . . . . . . . . . . . . . . 45

3.5 Exerćıcios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

i

ii SUMÁRIO

4 Métodos Aproximados 48

4.1 Computação Bayesiana . . . . . . . . . . . . . . . . . . . . . . . . 48

4.2 Uma Palavra de Cautela . . . . . . . . . . . . . . . . . . . . . . . 48

4.3 O Problema Geral da Inferência Bayesiana . . . . . . . . . . . . . 49

4.4 Método de Monte Carlo Simples . . . . . . . . . . . . . . . . . . . 50

4.4.1 Monte Carlo via Função de Importância . . . . . . . . . . 54

4.5 Métodos de Reamostragem . . . . . . . . . . . . . . . . . . . . . . 57

4.5.1 Método de Rejeição . . . . . . . . . . . . . . . . . . . . . . 57

4.5.2 Reamostragem Ponderada . . . . . . . . . . . . . . . . . . 60

4.6 Monte Carlo via cadeias de Markov . . . . . . . . . . . . . . . . . 63

4.6.1 Cadeias de Markov . . . . . . . . . . . . . . . . . . . . . . 63

4.6.2 Acurácia Numérica . . . . . . . . . . . . . . . . . . . . . . 64

4.6.3 Algoritmo de Metropolis-Hastings . . . . . . . . . . . . . . 65

4.6.4 Casos Especiais . . . . . . . . . . . . . . . . . . . . . . . . 71

4.6.5 Amostrador de Gibbs . . . . . . . . . . . . . . . . . . . . . 72

4.7 Problemas de Dimensão Variável . . . . . . . . . . . . . . . . . . 78

4.7.1 MCMC com Saltos Reversiveis (RJMCMC) . . . . . . . . 81

4.8 Tópicos Relacionados . . . . . . . . . . . . . . . . . . . . . . . . . 86

4.8.1 Autocorrelação Amostral . . . . . . . . . . . . . . . . . . . 86

4.8.2 Monitorando a Convergência . . . . . . . . . . . . . . . . . 86

5 Modelos Lineares 88

5.1 Análise de Variância com 1 Fator de Classificação . . . . . . . . . 91

A Lista de Distribuições 93

A.1 Distribuição Normal . . . . . . . . . . . . . . . . . . . . . . . . . 93

A.2 Distribuição Log-Normal . . . . . . . . . . . . . . . . . . . . . . . 94

A.3 A Função Gama . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94

A.4 Distribuição Gama . . . . . . . . . . . . . . . . . . . . . . . . . . 94

A.5 Distribuição Wishart . . . . . . . . . . . . . . . . . . . . . . . . . 95

A.6 Distribuição Gama Inversa . . . . . . . . . . . . . . . . . . . . . . 95

A.7 Distribuição Wishart Invertida . . . . . . . . . . . . . . . . . . . . 95

A.8 Distribuição Beta . . . . . . . . . . . . . . . . . . . . . . . . . . . 96

A.9 Distribuição de Dirichlet . . . . . . . . . . . . . . . . . . . . . . . 96

A.10 Distribuição t de Student . . . . . . . . . . . . . . . . . . . . . . . 96

A.11 Distribuição F de Fisher . . . . . . . . . . . . . . . . . . . . . . . 97

A.12 Distribuição de Pareto . . . . . . . . . . . . . . . . . . . . . . . . 97

A.13 Distribuição Binomial . . . . . . . . . . . . . . . . . . . . . . . . . 97

A.14 Distribuição Multinomial . . . . . . . . . . . . . . . . . . . . . . . 97

A.15 Distribuição de Poisson . . . . . . . . . . . . . . . . . . . . . . . . 98

A.16 Distribuição Binomial Negativa . . . . . . . . . . . . . . . . . . . 98

SUMÁRIO iii

B Alguns Endereços Interessantes 99

References 101

Caṕıtulo 1

Introdução

A informação que se tem sobre uma quantidade de interesse θ é fundamental na

Estat́ıstica. O verdadeiro valor de θ é desconhecido e a idéia é tentar reduzir

este desconhecimento. Além disso, a intensidade da incerteza a respeito de θ

pode assumir diferentes graus. Do ponto de vista Bayesiano, estes diferentes

graus de incerteza são representados através de modelos probabiĺısticos para θ.

Neste contexto, é natural que diferentes pesquisadores possam ter diferentes graus

de incerteza sobre θ (especificando modelos distintos). Sendo assim, não existe

nenhuma distinção entre quantidades observáveis e os parâmetros de um modelo

estat́ıstico, todos são considerados quantidades aleatórias.

1.1 Teorema de Bayes

Considere uma quantidade de interesse desconhecida θ (tipicamente não ob-

servável). A informação de que dispomos sobre θ, resumida probabilisticamente

através de p(θ), pode ser aumentada observando-se uma quantidade aleatória X

relacionada com θ. A distribuição amostral p(x|θ) define esta relação. A idéia de que após observar X = x a quantidade de informação sobre θ aumenta é bastante

intuitiva e o teorema de Bayes é a regra de atualização utilizada para quantificar

este aumento de informação,

p(θ|x) = p(x, θ) p(x)

= p(x|θ)p(θ)

p(x) =

p(x|θ)p(θ) ∫

p(θ, x)dθ . (1.1)

Note que 1/p(x), que não depende de θ, funciona como uma constante norma-

lizadora de p(θ|x). Para um valor fixo de x, a função l(θ; x) = p(x|θ) fornece a plausibilidade ou

verossimilhança de cada um dos posśıveis valores de θ enquanto p(θ) é chamada

distribuição a priori de θ. Estas duas fontes de informação, priori e verossimi-

1

2 CAPÍTULO 1. INTRODUÇÃO

lhança, são combinadas levando à distribuição a posteriori de θ, p(θ|x). Assim, a forma usual do teorema de Bayes é

p(θ|x) ∝ l(θ; x)p(θ), (1.2)

(lê-se p(θ|x) é proporcional a l(θ; x)p(θ)). Em palavras temos que

distribuição a posteriori ∝ verossimilhança× distribuição a priori.

Note que, ao omitir o termo p(x), a igualdade em (1.1) foi substituida por

uma proporcionalidade. Esta forma simplificada do teorema de Bayes será útil em

problemas que envolvam estimação de parâmetros já que o denominador é apenas

uma constante normalizadora. Em outras situações, como seleção e comparação

de modelos, este termo tem um papel crucial.

É intuitivo também que a probabilidade a posteriori de um particular conjunto

de valores de θ será pequena se p(θ) ou l(θ; x) for pequena para este conjunto. Em

particular, se atribuirmos probabilidade a priori igual a zero para um conjunto

de valores de θ então a probabilidade a posteriori será zero qualquer que seja a

amostra observada.

A partir da forma (1.2) a constante normalizadora da posteriori em (1.1) é

recuperada como

p(x) =

p(x, θ)dθ =

p(x|θ)p(θ)dθ = Eθ[p(X|θ)]

que é chamada distribuição preditiva. Esta é a distribuição esperada para a

observação x dado θ. Assim,

ˆ Antes de observar X podemos checar a adequação da priori fazendo

predições via p(x).

ˆ Se X observado recebia pouca probabilidade preditiva então o modelo deve

ser questionado.

Em muitas aplicações (e.g. séries temporais e geoestatistica) o maior inter-

esse é na previsão do processo em pontos não observados do tempo ou espaço.

Suponha então que, após observar X = x, estamos interessados na previsão de

uma quantidade Y , também relacionada com θ, e descrita probabilisticamente

por p(y|x, θ). A distribuição preditiva de Y dado x é obtida por integração como

p(y|x) = ∫

p(y, θ|x)dθ = ∫

p(y|θ, x)p(θ|x)dθ. (1.3)

Em muitos problemas estatisticos a hipótese de independência condicional entre

1.1. TEOREMA DE BAYES 3

X e Y dado θ está presente e a distribuição preditiva fica

p(y|x) = ∫

p(y|θ)p(θ|x)dθ.

Note no entanto que esta não é uma hipótese razoável para dados espacialmente

distribuidos aonde estamos admitindo que exista alguma estrutura de correlação

no espaço. De qualquer modo, em muitas aplicações práticas a integral em (1.3)

não tem solução analitica e precisaá ser obtida por algum método de aproximação.

Note também que as previsões são sempre verificáveis uma vez que Y é uma

quantidade observável. Finalmente, segue da última equação que

p(y|x) = Eθ|x[p(Y |θ)].

Fica claro também que os conceitos de priori e posteriori são relativos àquela

observação que está sendo considerada no momento. Assim, p(θ|x) é a posteriori de θ em relação a X (que já foi observado) mas é a priori de θ em relação a Y (que

não foi observado ainda). Após observar Y = y uma nova posteriori (relativa a

X = x e Y = y) é obtida aplicando-se novamente o teorema de Bayes. Mas será

que esta posteriori final depende da ordem em que as observações x e y foram

processadas? Observando-se as quantidades x1, x2, · · · , xn, independentes dado θ e relacionadas a θ através de pi(xi|θ) segue que

p(θ|x1) ∝ l1(θ; x1)p(θ) p(θ|x2, x1) ∝ l2(θ; x2)p(θ|x1)

∝ l2(θ; x2)l1(θ; x1)p(θ) ...

...

p(θ|xn, xn−1, · · · , x1) ∝ [

n ∏

i=1

li(θ; xi)

]

p(θ)

∝ ln(θ; xn) p(θ|xn−1, · · · , x1).

Ou seja, a ordem em que as observações são processadas pelo teorema de Bayes

é irrelevante. Na verdade, elas podem até ser processadas em subgrupos.

Exemplo 1.1 : (Gamerman e Migon, 1993) Um médico, ao examinar uma pes-

soa, “desconfia” que ela possa ter uma certa doença. Baseado na sua experiência,

no seu conhecimento sobre esta doença e nas informações dadas pelo paciente ele

assume que a probabilidade do paciente ter a doença é 0,7. Aqui a quantidade

4 CAPÍTULO 1. INTRODUÇÃO

de interesse desconhecida é o indicador de doença

θ =

{

1, se o paciente tem a doença

0, se o paciente não tem a doença.

Para aumentar sua quantidade de informação sobre a doença o médico aplica um

teste X relacionado com θ através da distribuição

P (X = 1 | θ = 0) = 0, 40 e P (X = 1 | θ = 1) = 0, 95

e o resultado do teste foi positivo (X = 1).

É bem intuitivo que a probabilidade de doença deve ter aumentado após este

resultado e a questão aqui é quantificar este aumento. Usando o teorema de Bayes

segue que

P (θ = 1 | X = 1) ∝ l(θ = 1;X = 1) p(θ = 1) = (0, 95)(0, 7) = 0, 665

P (θ = 0 | X = 1) ∝ l(θ = 0;X = 1) p(θ = 0) = (0, 40)(0, 3) = 0, 120.

Uma vez que as probabilidades a posteriori somam 1, i.e.

P (θ = 0 | X = 1) + P (θ = 1 | X = 1) = 1,

a constante normalizadora é obtida fazendo-se 0, 665k + 0, 120k = 1 e então

k = 1/0, 785. Portanto, a distribuição a posteriori de θ é

P (θ = 1 | X = 1) = 0, 665/0, 785 = 0, 847

P (θ = 0 | X = 1) = 0, 120/0, 785 = 0, 153.

O aumento na probabilidade de doença não foi muito grande porque a verossimil-

hança l(θ = 0;X = 1) também era grande (o modelo atribuia uma plausibilidade

grande para θ = 0 mesmo quando X = 1).

Agora o médico aplica outro teste Y cujo resultado está relacionado a θ através

da seguinte distribuição

P (Y = 1 | θ = 0) = 0, 04 e P (Y = 1 | θ = 1) = 0, 99.

Mas antes de observar o resultado deste teste é interessante obter sua distribuição

preditiva. Como θ é uma quantidade discreta segue que

p(y|x) = 1

θ=0

p(y|θ)p(θ|x)

1.1. TEOREMA DE BAYES 5

e note que p(θ|x) é a priori em relação a Y . Assim,

P (Y = 1 | X = 1) = P (Y = 1 | θ = 0)P (θ = 0 | X = 1) + P (Y = 1 | θ = 1)P (θ = 1 | X = 1) = (0, 04)(0, 153) + (0, 99)(0, 847) = 0, 845

P (Y = 0 | X = 1) = 1− P (Y = 1 | X = 1) = 0, 155.

O resultado deste teste foi negativo (Y = 0). Neste caso, é também intuitivo

que a probabilidade de doença deve ter diminuido e esta redução será quantificada

por uma nova aplicação do teorema de Bayes,

P (θ = 1 | X = 1, Y = 0) ∝ l(θ = 1;Y = 0)P (θ = 1 | X = 1) ∝ (0, 01)(0, 847) = 0, 0085

P (θ = 0 | X = 1, Y = 0) ∝ l(θ = 0;Y = 0)P (θ = 0 | X = 1) ∝ (0, 96)(0, 153) = 0, 1469.

A constante normalizadora é 1/(0,0085+0,1469)=1/0,1554 e assim a distribuição

a posteriori de θ é

P (θ = 1 | X = 1, Y = 0) = 0, 0085/0, 1554 = 0, 055

P (θ = 0 | X = 1, Y = 0) = 0, 1469/0, 1554 = 0, 945.

Verifique como a probabilidade de doença se alterou ao longo do experimento

P (θ = 1) =

0, 7, antes dos testes

0, 847, após o teste X

0, 055, após X e Y .

Note também que o valor observado de Y recebia pouca probabilidade preditiva.

Isto pode levar o médico a repensar o modelo, i.e.,

(i) Será que P (θ = 1) = 0, 7 é uma priori adequada?

(ii) Será que as distribuições amostrais de X e Y estão corretas ? O teste X é

tão inexpressivo e Y é realmente tão poderoso?

6 CAPÍTULO 1. INTRODUÇÃO

Exemplo 1.2 : Seja Y ∼ Binomial(12, θ) e em um experimento observou-se Y = 9. A função de verossimilhança de θ é dada por

l(θ) =

(

12

9

)

θ 9(1− θ)3, θ ∈ (0, 1).

Que distribuição poderia ser usada para resumir probabilisticamente nosso

conhecimento sobre o parâmetro θ? Note que, como 0 < θ < 1 queremos que,

p(θ) = 0 ⇒ p(θ|y) = 0, ∀θ ∋ (0, 1).

Podemos por exemplo assumir que θ ∼ N(µ, σ2) truncada no intervalo (0,1). Neste caso, denotando por fN(·|µ, σ2) a função de densidade da distribuição N(µ, σ2) segue que a função de densidade a priori de θ é dada por

p(θ) = fN(θ|µ, σ2)

∫ 1

0

fN(θ|µ, σ2)dθ .

Na Figura 1.1 esta função de densidade está representada para alguns valores de

µ e σ2. Os comandos do R abaixo podem ser utilizados para gerar as curvas. Note

como informações a priori bastante diferentes podem ser representadas.

> dnorm.t <- function(x, mean = 0, sd = 1) {

+ aux = pnorm(1, mean, sd) - pnorm(0, mean, sd)

+ dnorm(x, mean, sd)/aux

+ }

Outra possibilidade é através de uma reparametrização. Assumindo-se que

δ ∼ N(µ, σ2) e fazendo a transformação

θ = exp(δ)

1 + exp(δ)

segue que a transformação inversa é simplesmente

δ = log

(

θ

1− θ

)

= logito(θ).

Portanto a densidade a priori de θ fica

p(θ) = fN(δ(θ)|µ, σ2) ∣

= (2πσ2)−1/2 exp

{

− 1 2σ2

(

log

(

θ

1− θ

)

− µ )2

}

1

θ(1− θ)

1.1. TEOREMA DE BAYES 7

0.0 0.2 0.4 0.6 0.8 1.0

0. 0

0. 5

1. 0

1. 5

2. 0

2. 5

3. 0

θ

p( θ)

N(0.5,0.5) N(0,0.5) N(1,0.5) N(2,0.5)

Figura 1.1: Densidades a priori normais truncadas para o parametro θ no Exemplo 1.2.

e é chamada de normal-logistica. Na Figura 1.2 esta função de densidade está

representada para alguns valores de µ e σ2. Os comandos do R abaixo foram

utilizados. Novamente note como informações a priori bastante diferentes podem

ser representadas. Em particular a função de densidade de θ será sempre unimodal

quando σ2 ≤ 2 e bimodal quando σ2 > 2.

> dlogist = function(x, mean, sd) {

+ z = log(x/(1 - x))

+ dnorm(z, mean, sd)/(x - x^2)

+ }

Finalmente, podemos atribuir uma distribuição a priori θ ∼ Beta(a, b) (ver Apêndice A),

p(θ) = Γ(a+ b)

Γ(a)Γ(b) θa−1(1− θ)b−1, a > 0, b > 0, θ ∈ (0, 1).

Esta distribuição é simétrica em torno de 0,5 quando a = b e assimétrica quando

a 6= b. Variando os valores de a e b podemos definir uma rica familia de dis- tribuições a priori para θ, incluindo a distribuição Uniforme no intervalo (0,1) se

a = b = 1. Algumas possibilidades estão representadas na Figura 1.3.

Um outro resultado importante ocorre quando se tem uma única observação

da distribuição normal com média desconhecida. Se a média tiver priori normal

8 CAPÍTULO 1. INTRODUÇÃO

0.0 0.2 0.4 0.6 0.8 1.0

0 1

2 3

4

θ

p( θ)

N(−1,0.25) N(1,1) N(0,4)

Figura 1.2: Densidades a priori tipo logisticas para o parâmetro θ no Exemplo 1.2.

então os parâmetros da posteriori são obtidos de uma forma bastante intuitiva

como visto no teorema a seguir.

Teorema 1.1 Se X|θ ∼ N(θ, σ2) sendo σ2 conhecido e θ ∼ N(µ0, τ 20 ) então θ|x ∼ N(µ1, τ 21 ) sendo

µ1 = τ−20 µ0 + σ

−2x

τ−20 + σ −2

e τ−21 = τ −2 0 + σ

−2.

Prova. Temos que

p(x|θ) ∝ exp{−σ−2(x− θ)2/2} e p(θ) ∝ exp{−τ−20 (θ − µ0)/2}

e portanto

p(θ|x) ∝ exp {

−1 2 [σ−2(θ2 − 2xθ) + τ−20 (θ2 − 2µ0θ)]

}

∝ exp {

−1 2 [θ2(σ−2 + τ−20 )− 2θ(σ−2x+ τ−20 µ0)]

}

.

sendo que os termos que não dependem de θ foram incorporados à constante de

proporcionalidade. Definindo τ−21 = σ −2+ τ−20 e τ

−2 1 µ1 = σ

−2x− τ−20 µ0 segue que

p(θ|x) ∝ exp {

−τ −2 1

2 (θ2 − 2θµ1)

}

∝ exp {

−τ −2 1

2 (θ − µ1)2

}

pois µ1 não depende de θ. Portanto, a função de densidade a posteriori (a menos

1.1. TEOREMA DE BAYES 9

0.0 0.2 0.4 0.6 0.8 1.0

0 1

2 3

4 5

θ

p( θ)

Beta(1.5,4) Beta(2,0.5) Beta(7,1.5) Beta(3,3)

Figura 1.3: Densidades a priori Beta para o parâmetro θ no Exemplo 1.2.

de uma constante) tem a mesma forma de uma normal com média µ1 e variância

τ 21 .

Note que, definindo precisão como o inverso da variância, segue do teorema

que a precisão a posteriori é a soma das precisões a priori e da verossimilhança

e não depende de x. Interpretando precisão como uma medida de informação

e definindo w = τ−20 /(τ −2 0 + σ

−2) ∈ (0, 1) então w mede a informação relativa contida na priori com respeito à informação total. Podemos escrever então que

µ1 = wµ0 + (1− w)x

ou seja, µ1 é uma combinação linear convexa de µ0 e x e portanto

min{µ0, x} ≤ µ1 ≤ max{µ0, x}.

A distribuição preditiva de X também é facilmente obtida notando que pode-

mos reescrever as informações na forma de equações com erros não correlaciona-

dos. Assim,

X = θ + ǫ, ǫ ∼ N(0, σ2) θ = µ0 + w, w ∼ N(0, τ 20 )

tal que Cov(θ, ǫ) = Cov(µ0, w) = 0. Portanto a distribuição (incondicional) de

X é normal pois ele resulta de uma soma de variáveis aleatórias com distribuição

10 CAPÍTULO 1. INTRODUÇÃO

normal. Além disso,

E(X) = E(θ) + E(ǫ) = µ0

V ar(X) = V ar(θ) + V ar(ǫ) = τ 20 + σ 2

Conclusão, X ∼ N(µ0, τ 20 + σ2).

Exemplo 1.3 : (Box & Tiao, 1992) Os f́ısicos A e B desejam determinar uma

constante f́ısica θ. O f́ısico A tem mais experiência nesta área e especifica sua

priori como θ ∼ N(900, 202). O f́ısico B tem pouca experiência e especifica uma priori muito mais incerta em relação à posição de θ, θ ∼ N(800, 802). Assim, não é dificil verificar que

para o fisico A: P (860 < θ < 940) ≈ 0, 95

para o fisico B: P (640 < θ < 960) ≈ 0, 95.

Faz-se então uma medição X de θ em laboratório com um aparelho calibrado

com distribuição amostral X|θ ∼ N(θ, 402) e observou-se X = 850. Aplicando o teorema 1.1 segue que

(θ|X = 850) ∼ N(890, 17, 92) para o f́ısico A

(θ|X = 850) ∼ N(840, 35, 72) para o f́ısico B.

Note também que os aumentos nas precisões a posteriori em relação às precisões

a priori foram,

ˆ para o f́ısico A: precisão(θ) passou de τ−20 = 0, 0025 para τ −2 1 = 0, 00312

(aumento de 25%).

ˆ para o f́ısico B: precisão(θ) passou de τ−20 = 0, 000156 para τ −2 1 = 0, 000781

(aumento de 400%).

A situação está representada graficamente na Figura 1.4 a seguir. Note como a

distribuição a posteriori representa um compromisso entre a distribuição a priori

e a verossimilhança. Além disso, como as incertezas iniciais são bem diferentes

o mesmo experimento fornece muito pouca informação adicional para o fisico A

enquanto que a incerteza do fisico B foi bastante reduzida. Os comandos do R

abaixo podem ser usados nos cálculos.

1.2. PRINCÍPIO DA VEROSSIMILHANÇA 11

> norm.norm <- function(x, mu0, tau0, s0) {

+ precisao = 1/tau0 + length(x)/s0

+ tau1 = 1/precisao

+ w = (1/tau0)/precisao

+ mu1 = w * mu0 + (1 - w) * mean(x)

+ return(list(m = mu1, tau = tau1))

+ }

700 750 800 850 900 950 1000

0. 00

0 0.

00 5

0. 01

0 0.

01 5

0. 02

0

θ

priori posteriori verossimilhanca Fisico A

Fisico B

Figura 1.4: Densidades a priori e a posteriori e função de verossimilhança para o Exemplo 1.3.

1.2 Prinćıpio da Verossimilhança

O exemplo a seguir (DeGroot, 1970, páginas 165 e 166) ilustra esta propriedade.

Imagine que cada item de uma população de itens manufaturados pode ser clas-

sificado como defeituoso ou não defeituoso. A proporção θ de itens defeituosos

na população é desconhecida e uma amostra de itens será selecionada de acordo

com um dos seguintes métodos:

(i) n itens serão selecionados ao acaso.

(ii) Itens serão selecionados ao acaso até que y defeituosos sejam obtidos.

(iii) Itens serão selecionados ao acaso até que o inspetor seja chamado para

resolver um outro problema.

12 CAPÍTULO 1. INTRODUÇÃO

(iv) Itens serão selecionados ao acaso até que o inspetor decida que já acumulou

informação suficiente sobre θ.

Qualquer que tenha sido o esquema amostral, se foram inspecionados n itens

x1, · · · , xn dos quais y eram defeituosos então

l(θ; x) ∝ θy(1− θ)n−y.

O Prinćıpio da Verossimilhança postula que para fazer inferência sobre uma

quantidade de interesse θ só importa aquilo que foi realmente observado e não

aquilo que “poderia” ter ocorrido mas efetivamente não ocorreu.

1.3 Exerćıcios

1. No Exemplo 1.3, obtenha também a distribuição preditiva de X e compare

o valor observado com a média desta preditiva para os 2 f́ısicos. Faça uma

previsão para uma 2a medição Y feita com o mesmo aparelho.

2. Uma máquina produz 5% de itens defeituosos. Cada item produzido passa

por um teste de qualidade que o classifica como “bom”, “defeituoso” ou

“suspeito”. Este teste classifica 20% dos itens defeituosos como bons e 30%

como suspeitos. Ele também classifica 15% dos itens bons como defeituosos

e 25% como suspeitos.

(a) Que proporção dos itens serão classificados como suspeitos ?

(b) Qual a probabilidade de um item classificado como suspeito ser de-

feituoso ?

(c) Outro teste, que classifica 95% dos itens defeituosos e 1% dos itens

bons como defeituosos, é aplicado somente aos itens suspeitos.

(d) Que proporção de itens terão a suspeita de defeito confirmada ?

(e) Qual a probabilidade de um item reprovado neste 2o teste ser defeituoso

?

3. Uma empresa de crédito precisa saber como a inadimplência está distribuida

entre seus clentes. Sabe-se que um cliente pode pertencer às classes A, B,

C ou D com probabilidades 0,50, 0,20, 0,20 e 0,10 respectivamente. Um

cliente da classe A tem probabilidade 0,30 de estar inadimplente, um da

classe B tem probabilidade 0,10 de estar inadimplente, um da classe C tem

probabilidade 0,05 de estar inadimplente e um da classe D tem probabili-

dade 0,05 de estar inadimplente. Um cliente é sorteado aleatoriamente.

(a) Defina os eventos e enumere as probabilidades fornecidas no problema.

1.3. EXERCÍCIOS 13

(b) Qual a probabilidade dele estar inadimplente ?

(c) Sabendo que ele está inadimplente, qual a sua classe mais provável?

4. Suponha que seus dados x1, . . . , xn são processados sequencialmente, i.e. x1 é observado antes de x2 e assim por diante. Escreva um programa que aplica

o Teorema 1.1 obtendo a média e a variância a posteriori dado x1, use esta

distribuição como priori para obter a média e a variância a posteriori dados

x1, x2 e repita o procedimento sequencialmente até obter a posteriori dados

x1, . . . , xn. Faça um gráfico com as médias a posteriori mais ou menos 2

desvios padrão a posteriori.

Caṕıtulo 2

Distribuições a Priori

A utilização de informação a priori em inferência Bayesiana requer a especificação

de uma distribuição a priori para a quantidade de interesse θ. Esta distribuição

deve representar (probabilisticamente) o conhecimento que se tem sobre θ antes

da realização do experimento. Neste capitulo serão discutidas diferentes formas

de especificação da distribuição a priori.

2.1 Prioris Conjugadas

A partir do conhecimento que se tem sobre θ, pode-se definir uma famı́lia

paramétrica de densidades. Neste caso, a distribuição a priori é representada

por uma forma funcional, cujos parâmetros devem ser especificados de acordo

com este conhecimento. Estes parâmetros indexadores da familia de distribuições

a priori são chamados de hiperparâmetros para distingui-los dos parâmetros de

interesse θ.

Esta abordagem em geral facilita a análise e o caso mais importante é o de

prioris conjugadas. A idéia é que as distribuições a priori e a posteriori pertençam

a mesma classe de distribuições e assim a atualização do conhecimento que se tem

de θ envolve apenas uma mudança nos hiperparâmetros. Neste caso, o aspecto

sequencial do método Bayesiano pode ser explorado definindo-se apenas a regra de

atualização dos hiperparâmetros já que as distribuições permanecem as mesmas.

Definição 2.1 Se F = {p(x|θ), θ ∈ Θ} é uma classe de distribuições amostrais então uma classe de distribuições P é conjugada a F se

∀ p(x|θ) ∈ F e p(θ) ∈ P ⇒ p(θ|x) ∈ P.

Gamerman (1996, 1997 Cap. 2) alerta para o cuidado com a utilização in-

discriminada de prioris conjugadas. Essencialmente, o problema é que a priori

14

2.2. CONJUGAÇÃO NA FAMÍLIA EXPONENCIAL 15

conjugada nem sempre é uma representação adequada da incerteza a priori. Sua

utilização está muitas vezes associada à tratabilidade anaĺıtica decorrente.

Uma vez entendidas suas vantagens e desvantagens a questão que se coloca

agora é “como” obter uma famı́lia de distribuições conjugadas.

(i) Identifique a classe P de distribuições para θ tal que l(θ; x) seja proporcional

a um membro desta classe.

(ii) Verifique se P é fechada por amostragem, i.e., se ∀ p1, p2 ∈ P ∃ k tal que kp1p2 ∈ P .

Se, além disso, existe uma constante k tal que k−1 = ∫

l(θ; x)dθ < ∞ e todo p ∈ P é definido como p(θ) = k l(θ; x) então P é a famı́lia conjugada natural ao modelo amostral gerador de l(θ; x).

Exemplo 2.1 : Sejam X1, . . . , Xn ∼ Bernoulli(θ). Então a densidade amostral conjunta é

p(x|θ) = θt(1− θ)n−t, 0 < θ < 1 sendo t = n

i=1

xi

e pelo teorema de Bayes segue que

p(θ|x) ∝ θt(1− θ)n−tp(θ).

Note que l(θ; x) é proporcional à densidade de uma distribuição

Beta(t + 1, n − t + 1). Além disso, se p1 e p2 são as densidades das dis- tribuições Beta(a1, b1) e Beta(a2, b2) então

p1p2 ∝ θa1+a2−2(1− θ)b1+b2−2,

ou seja p1p2 é proporcional a densidade da distribuição Beta(a1 + a2 − 1, b1 + b2 − 1). Conclui-se que a famı́lia de distribuições Beta com parâmetros inteiros é conjugada natural à famı́lia Bernoulli. Na prática esta classe pode ser ampliada

para incluir todas as distribuições Beta, i.e. incluindo todos os valores positivos

dos parâmetros.

2.2 Conjugação na Famı́lia Exponencial

A famı́lia exponencial inclui muitas das distribuições de probabilidade mais comu-

mente utilizadas em Estatistica, tanto continuas quanto discretas. Uma caracter-

istica essencial desta familia é que existe uma estatistica suficiente com dimensão

16 CAPÍTULO 2. DISTRIBUIÇÕES A PRIORI

fixa. Veremos adiante que a classe conjugada de distribuições é muito fácil de

caracterizar.

Definição 2.2 A familia de distribuições com função de (densidade) de probabil-

idade p(x|θ) pertence à familia exponencial a um parâmetro se podemos escrever

p(x|θ) = a(x) exp{u(x)φ(θ) + b(θ)}.

Note que pelo critério de fatoração de Neyman U(x) é uma estatistica suficiente

para θ.

Neste caso, a classe conjugada é facilmente identificada como,

p(θ) = k(α, β) exp{αφ(θ) + βb(θ)}.

e aplicando o teorema de Bayes segue que

p(θ|x) = k(α + u(x), β + 1) exp{[α + u(x)]φ(θ) + [β + 1]b(θ)}.

Agora, usando a constante k, a distribuição preditiva pode ser facilmente obtida

sem necessidade de qualquer integração. A partir da equação p(x)p(θ|x) = p(x|θ)p(θ) e após alguma simplificação segue que

p(x) = p(x|θ)p(θ) p(θ|x) =

a(x)k(α, β)

k(α + u(x), β + 1) .

Exemplo 2.2 : Uma extensão direta do Exemplo 2.1 é o modelo binomial, i.e.

X|θ ∼ Binomial(n, θ). Neste caso,

p(x|θ) = (

n

x

)

exp

{

x log

(

θ

1− θ

)

+ n log(1− θ) }

e a famı́lia conjugada natural é Beta(r, s). Podemos escrever então

p(θ) ∝ θr−1(1− θ)s−1

∝ exp {

(r − 1) log (

θ

1− θ

)

+

(

s+ r − 2 n

)

n log(1− θ) }

∝ exp {αφ(θ) + βb(θ)} .

A posteriori também é Beta com parâmetros α + x e β + 1 ou equivalentemente

2.2. CONJUGAÇÃO NA FAMÍLIA EXPONENCIAL 17

r + x e s+ n− x, i.e.

p(θ|x) ∝ exp {

(r + x− 1)φ(θ) + [

s+ r − 2 + n n

]

b(θ)

}

∝ θr+x−1(1− θ)s+n−x−1.

Como ilustração, no Exemplo 2.2 suponha que n = 12, X = 9 e usamos pri-

oris conjugadas Beta(1,1), Beta(2,2) e Beta(1,3). As funções de densidade destas

distribuições juntamente com a função de verossimilhança normalizada e as re-

spectivas densidades a posteriori estão na Figura 2.1. A distribuição preditiva é

dada por

p(x) =

(

n

x

)

B(r + x, s+ n− x) B(r, s)

, x = 0, 1, . . . , n, n ≥ 1,

onde B−1 é a constante normalizadora da distribuição Beta, i.e. (ver Apêndice

A)

B−1(a, b) = Γ(a+ b)

Γ(a)Γ(b) .

Esta distribuição é denominada Beta-Binomial.

0.0 0.2 0.4 0.6 0.8 1.0

0. 0

1. 0

2. 0

3. 0

θ

veross priori posteriori

0.0 0.2 0.4 0.6 0.8 1.0

0. 0

1. 0

2. 0

3. 0

θ

veross priori posteriori

0.0 0.2 0.4 0.6 0.8 1.0

0. 0

1. 0

2. 0

3. 0

θ

veross priori posteriori

Figura 2.1: Densidades a priori, a posteriori e função de verossimilhança normalizada para o Exemplo 2.2.

18 CAPÍTULO 2. DISTRIBUIÇÕES A PRIORI

No Exemplo 2.2 suponha novamente que n = 12, X = 9 e usamos as prioris

conjugadas Beta(1,1), Beta(2,2) e Beta(1,3). Na Tabela 2.1 estão listadas as

probabilidades preditivas P (X = k) associadas a estas prioris. Os comandos do

R a seguir podem ser usados no cálculo destas probabilidades.

> beta.binomial = function(n, a, b) {

+ m = matrix(0, n + 1, 2)

+ m[, 1] = 0:n

+ for (x in 0:n) m[x, 2] = round(choose(n, x) * beta(a + x,

+ b + n - x)/beta(a, b), 4)

+ return(list(m = m))

+ }

Tabela 2.1: Probabilidades preditivas da Beta-Binomial para o Exemplo 2.2

k Beta(1,1) Beta(2,2) Beta(1,3) 0 0.0769 0.0527 0.1714 1 0.0769 0.0725 0.1451 2 0.0769 0.0879 0.1209 3 0.0769 0.0989 0.0989 4 0.0769 0.1055 0.0791 5 0.0769 0.1077 0.0615 6 0.0769 0.1055 0.0462 7 0.0769 0.0989 0.0330 8 0.0769 0.0879 0.0220 9 0.0769 0.0725 0.0132 10 0.0769 0.0527 0.0066 11 0.0769 0.0286 0.0022 12 0.0000 0.0000 0.0000

No caso geral em que se tem uma amostra X1, . . . , Xn da famı́lia exponencial

a natureza sequencial do teorema de Bayes permite que a análise seja feita por

replicações sucessivas. Assim a cada observação xi os parâmetros da distribuição

a posteriori são atualizados via

αi = αi−1 + u(xi)

βi = βi−1 + 1

2.3. PRINCIPAIS FAMÍLIAS CONJUGADAS 19

com α0 = α e β0 = β. Após n observações temos que

αn = α + n

i=1

u(xi)

βn = β + n

e a distribuição preditiva é dada por

p(x) =

[

n ∏

i=1

a(xi)

]

k(α, β)

k(α + ∑

u(xi), β + n) .

Finalmente, a definição de famı́lia exponencial pode ser extendida ao caso

multiparamétrico, i.e.

p(x|θ) = [

n ∏

i=1

a(xi)

]

exp

{

r ∑

j=1

[

n ∑

i=1

uj(xi)

]

φj(θ) + nb(θ)

}

com θ = (θ1, . . . , θr). Neste caso, pelo critério de fatoração, temos que ∑

U1(xi), . . . , ∑

Ur(xi) é uma estat́ıstica conjuntamente suficiente para o vetor

de parâmetros θ.

2.3 Principais Famı́lias Conjugadas

Já vimos que a famı́lia de distribuições Beta é conjugada ao modelo Bernoulli e

binomial. Não é dif́ıcil mostrar que o mesmo vale para as distribuições amostrais

geométrica e binomial-negativa (ver Exerćıcio 1). A seguir veremos resultados

para outros membros importantes da famı́lia exponencial.

2.3.1 Distribuição normal com variância conhecida

Para uma única observação vimos pelo Teorema 1.1 que a famı́lia de distribuições

normais é conjugada ao modelo normal. Para uma amostra de tamanho n, a

função de verossimilhança pode ser escrita como

l(θ; x) = (2πσ2)−n/2 exp

{

− 1 2σ2

n ∑

i=1

(xi − θ)2 }

∝ exp {

− n 2σ2

(x− θ)2 }

onde os termos que não dependem de θ foram incorporados à constante de pro-

porcionalidade. Portanto, a verossimilhança tem a mesma forma daquela baseada

em uma única observação bastando substituir x por x e σ2 por σ2/n. Logo vale

20 CAPÍTULO 2. DISTRIBUIÇÕES A PRIORI

o Teorema 1.1 com as devidas substituições, i.e. a distribuição a posteriori de θ

dado x é N(µ1, τ 2 1 ) sendo

µ1 = τ−20 µ0 + nσ

−2x

τ−20 + nσ −2

e τ−21 = τ −2 0 + nσ

−2.

Note que a média a posteriori pode ser reescrita como wµ0 + (1 − w)x sendo w = τ−20 /(τ

−2 0 + nσ

−2).

Uma função geral pode ser escrita no R para calcular estes parâmetros e op-

cionalmente fazer os gráficos das densidades.

> norm.norm <- function(x, sigma, mu0, tau0, plot = F) {

+ n = length(x)

+ xbar = mean(x)

+ ep = sigma/sqrt(n)

+ sigma2 = sigma^2

+ tau1 = n * (1/sigma2) + (1/tau0)

+ mu1 = (n * (1/sigma2) * xbar + (1/tau0) * mu0)/tau1

+ if (plot) {

+ curve(dnorm(x, xbar, ep), xbar - 3 * ep, xbar + 3 * ep)

+ curve(dnorm(x, mu0, sqrt(tau0)), add = T, col = 2)

+ curve(dnorm(x, mu1, 1/sqrt(tau1)), add = T, col = 3)

+ legend(-0.5, 1.2, legend = c("veross.", "priori", "posteriori"),

+ col = 1:3, lty = c(1, 1, 1))

+ }

+ return(list(mu1 = mu1, tau1 = tau1))

+ }

2.3.2 Distribuição de Poisson

Seja X1, . . . , Xn uma amostra aleatória da distribuição de Poisson com parâmetro

θ. Sua função de probabilidade conjunta é dada por

p(x|θ) = e −nθθt ∏

xi! ∝ e−nθθt, θ > 0, t =

n ∑

i=1

xi.

O núcleo da verossimilhança é da forma θae−bθ que caracteriza a famı́lia de

distribuições Gama a qual é fechada por amostragem (verifique!). Assim, dis-

tribuição a priori conjugada natural de θ é Gama com parâmetros positivos α e

β, i.e.

p(θ) = βα

Γ(α) θα−1e−βθ, α > 0, β > 0, θ > 0.

2.3. PRINCIPAIS FAMÍLIAS CONJUGADAS 21

A densidade a posteriori fica

p(θ|x) ∝ θα+t−1 exp {−(β + n)θ}

que corresponde à densidade Gama(α+ t, β + n). Note que a média a posteriori

pode ser reescrita como uma combinação linear da média a priori e da média

amostral (ver exerćıcio 6). A distribuição preditiva também é facilmente obtida

pois

p(x|θ) = [

n ∏

i=1

1

xi!

]

exp {t log θ − nθ}

e portanto

p(x) =

[

n ∏

i=1

1

xi!

]

βα

Γ(α)

Γ(α + t)

(β + n)α+t .

Para uma única observação x segue então que

p(x) = 1

x!

βα Γ(α + x)

Γ(α) (β + 1)α+x =

1

x!

(

β

β + 1

)α ( 1

β + 1

)x (α + x− 1)! (α− 1)!

=

(

α + x− 1 x

)(

β

β + 1

)α ( 1

β + 1

)x

.

Esta distribuição é chamada de Binomial-Negativa com parâmetros α e β e sua

média e variância são facilmente obtidos como

E(X) = E[E(X|θ)] = E(θ) = α/β

V ar(X) = E[V ar(X|θ)] + V ar[E(X|θ)] = E(θ) + V ar(θ) = α(β + 1) β2

.

2.3.3 Distribuição multinomial

Denotando por X = (X1, . . . , Xp) o número de ocorrências em cada uma de p

categorias em n ensaios independentes e por θ = (θ1, . . . , θp) as probabilidades

associadas, deseja-se fazer inferência sobre estes p parâmetros. No entanto, note

que existem efetivamente p − 1 parâmetros já que temos a seguinte restrição ∑p

i=1 θi = 1. Além disso, a restrição ∑p

i=1 Xi = n obviamente também se aplica.

Dizemos que X tem distribuição multinomial com parâmetros n e θ e função de

probabilidade conjunta das p contagens X é dada por

p(x|θ) = n!∏p i=1 xi!

p ∏

i=1

θxii .

22 CAPÍTULO 2. DISTRIBUIÇÕES A PRIORI

Note que esta é uma generalização da distribuição binomial que tem apenas duas

categorias. Não é dif́ıcil mostrar que esta distribuição também pertence à famı́lia

exponencial. A função de verossimilhança para θ é

l(θ;x) ∝ p ∏

i=1

θxii

que tem o mesmo núcleo da função de densidade de uma distribuição de Dirichlet.

A famı́lia Dirichlet com parâmetros inteiros a1, . . . , ap é a conjugada natural do

modelo multinomial, porém na prática a conjugação é extendida para parâmetros

não inteiros. A distribuição a posteriori é dada por

p(θ|x) ∝ p ∏

i=1

θxii

p ∏

i=1

θai−1i =

p ∏

i=1

θxi+ai−1i .

Note que estamos generalizando a análise conjugada para amostras binomiais com

priori beta.

2.3.4 Distribuição normal com média conhecida e variân-

cia desconhecida

Seja X1, . . . , Xn uma amostra aleatória da distribuição N(θ, σ 2), com θ conhecido

e φ = σ−2 desconhecido. Neste caso a função de densidade conjunta é dada por

p(x|θ, φ) ∝ φn/2 exp{−φ 2

n ∑

i=1

(xi − θ)2}.

Note que o núcleo desta verossimilhança tem a mesma forma daquele de uma

distribuição Gama. Como sabemos que a famı́lia Gama é fechada por amostragem

podemos considerar uma distribuição a priori Gama com parâmetros n0/2 e

n0σ 2 0/2, i.e.

φ ∼ Gama (

n0 2 , n0σ

2 0

2

)

.

Equivalentemente, podemos atribuir uma distribuição a priori qui-quadrado com

n0 graus de liberdade para n0σ 2 0φ. A forma funcional dos parâmetros da dis-

tribuição a priori é apenas uma conveniência matemática como veremos a seguir.

Definindo ns20 = ∑n

i=1(xi − θ)2 e aplicando o teorema de Bayes obtemos a

2.3. PRINCIPAIS FAMÍLIAS CONJUGADAS 23

distribuição a posteriori de φ,

p(φ|x) ∝ φn/2 exp {

−φ 2 ns20

}

φn0/2−1 exp

{

−φ 2 n0σ

2 0

}

= φ(n0+n)/2−1 exp

{

−φ 2 (n0σ

2 0 + ns

2 0)

}

.

Note que esta expressão corresponde ao núcleo da distribuição Gama, como

era esperado devido à conjugação. Portanto,

φ|x ∼ Gama (

n0 + n

2 , n0σ

2 0 + ns

2 0

2

)

.

Equivalentemente podemos dizer que (n0σ 2 0 + ns

2 0)φ | x ∼ χ2n0+n.

2.3.5 Distribuição normal com média e variância descon-

hecidos

Seja X1, . . . , Xn uma amostra aleatória da distribuição N(θ, σ 2), com ambos θ

e φ=σ−2 desconhecidos. Precisamos então especificar uma distribuição a priori

conjunta para θ e φ. Uma possibilidade é fazer a especificação em dois estágios

já que podemos sempre escrever p(θ, φ) = p(θ|φ)p(φ). No primeiro estágio,

θ|φ ∼ N(µ0, (c0φ)−1), φ = σ−2

e a distribuição a priori marginal de φ é a mesma do caso anterior, i.e.

φ ∼ Gama (

n0 2 , n0σ

2 0

2

)

.

A distribuição conjunta de (θ, φ) é geralmente chamada de Normal-Gama com

parâmetros (µ0, c0, n0, σ 2 0) e sua função de densidade conjunta é dada por,

p(θ, φ) = p(θ|φ)p(φ)

∝ φ1/2 exp {

−c0φ 2

(θ − µ0)2 }

φn0/2−1 exp

{

−n0σ 2 0φ

2

}

∝ φ(n0+1)/2−1 exp {

−φ 2 (n0σ

2 0 + c0(θ − µ0)2)

}

.

A partir desta densidade conjunta podemos obter a distribuição marginal de

24 CAPÍTULO 2. DISTRIBUIÇÕES A PRIORI

θ por integração

p(θ) =

p(θ|φ)p(φ)dφ

∝ ∫ ∞

0

φ(n0+1)/2−1 exp

{

−φ 2 [n0σ

2 0 + c0(θ − µ0)2]

}

∝ [

n0σ 2 0 + c0(θ − µ0)2

2

]− n0+1

2

∝ [

1 + (θ − µ0)2 n0(σ20/c0)

]− n0+1

2

,

que é o núcleo da distribuição t de Student com n0 graus de liberdade, parâmetro

de locação µ0 e parâmetro de escala σ 2 0/c0 (ver Apêndice A). Denotamos θ ∼

tn0(µ0, σ 2 0/c0). A distribuição condicional de φ dado θ também é facilmente obtida

como

p(φ|θ) ∝ p(θ|φ)p(φ)

∝ φ(n0+1)/2−1 exp {

−φ 2 [n0σ

2 0 + c0(θ − µ0)2]

}

,

e portanto,

φ|θ ∼ Gama (

n0 + 1

2 , n0σ

2 0 + c0(θ − µ0)2

2

)

.

A posteriori conjunta de (θ, φ) também é obtida em 2 etapas como segue.

Primeiro, para φ fixo podemos usar o resultado da Seção 2.3.1 de modo que a

distribuição a posteriori de θ dado φ fica

θ|φ,x ∼ N(µ1, (c1φ)−1)

sendo

µ1 = c0φµ0 + nφx

c0φ+ nφ =

c0µ0 + nx

c0 + n e c1 = c0 + n.

Na segunda etapa, combinando a verossimilhança com a priori de φ obtemos que

φ|x ∼ Gama (

n1 2 , n1σ

2 1

2

)

sendo

n1 = n0 + n e n1σ 2 1 = n0σ

2 0 +

(xi − x)2 + c0n(µ0 − x)2/(c0 + n).

Equivalentemente, podemos escrever a posteriori de φ como n1σ 2 1φ ∼ χ2n1 . As-

sim, a posteriori conjunta é (θ, φ|x) ∼ Normal-Gama(µ1, c1, n1, σ21) e portanto a

2.4. PRIORI NÃO INFORMATIVA 25

posteriori marginal de θ fica

θ | x ∼ tn1(µ1, σ21/c1).

Em muitas situações é mais fácil pensar em termos de algumas caracteŕısticas

da distribuição a priori do que em termos de seus hiperparâmetros. Por exemplo,

se E(θ) = 2, V ar(θ) = 5, E(φ) = 3 e V ar(φ) = 3 então

(i) µ0 = 2 pois E(θ) = µ0.

(ii) σ20 = 1/3 pois E(φ) = 1/σ 2 0.

(iii) n0 = 6 pois V ar(φ) = 2/(n0σ 4 0) = 18/n0.

(iv) c0 = 1/10 pois V ar(θ) =

(

n0 n0 − 2

)

σ20 c0

= 1

2c0

2.4 Priori não Informativa

Esta seção refere-se a especificação de distribuições a priori quando se espera que

a informação dos dados seja dominante, no sentido de que a nossa informação

a priori é vaga. Os conceitos de “conhecimento vago”, “não informação”, ou “ig-

norância a priori” claramente não são únicos e o problema de caracterizar prioris

com tais caracteŕısticas pode se tornar bastante complexo.

Por outro lado, reconhece-se a necessidade de alguma forma de análise que,

em algum sentido, consiga captar esta noção de uma priori que tenha um efeito

mı́nimo, relativamente aos dados, na inferência final. Tal análise pode ser pen-

sada como um ponto de partida quando não se consegue fazer uma elicitação

detalhada do “verdadeiro” conhecimento a priori. Neste sentido, serão apresen-

tadas aqui algumas formas de “como” fazer enquanto discussões mais detalhadas

são encontradas em Berger (1985), Box & Tiao (1992), Bernardo & Smith (1994)

e O’Hagan (1994).

A primeira idéia de “não informação” a priori que se pode ter é pensar em

todos os posśıveis valores de θ como igualmente prováveis, i.e. com uma dis-

tribuição a priori uniforme. Neste caso, fazendo p(θ) ∝ k para θ variando em um subconjunto da reta significa que nenhum valor particular tem preferência (Bayes,

1763). Porém esta escolha de priori pode trazer algumas dificuldades técnicas,

(i) Se o intervalo de variação de θ for ilimitado então a distribuição a priori é

imprópria, i.e. ∫

p(θ)dθ = ∞.

26 CAPÍTULO 2. DISTRIBUIÇÕES A PRIORI

(ii) Se φ = g(θ) é uma reparametrização não linear monótona de θ então p(φ) é

não uniforme já que pelo teorema de transformação de variáveis

p(φ) = p(θ(φ))

∝ ∣

.

Na prática, como estaremos interessados na distribuição a posteriori não dare-

mos muita importância à impropriedade da distribuição a priori. No entanto de-

vemos sempre nos certificar de que a posterior é própria antes de fazer qualquer

inferência.

A classe de prioris não informativas proposta por Jeffreys (1961) é invariante

a transformações 1 a 1, embora em geral seja imprópria e será definida a seguir.

Antes porém precisamos da definição da medida de informação de Fisher.

Definição 2.3 Considere uma única observação X com função de (densidade)

de probabilidade p(x|θ). A medida de informação esperada de Fisher de θ através de X é definida como

I(θ) = E

[

−∂ 2 log p(x|θ)

∂θ2

]

.

Se θ for um vetor paramétrico define-se então a matriz de informação esperada

de Fisher de θ através de X como

I(θ) = E

[

−∂ 2 log p(x|θ) ∂θ∂θ′

]

.

Note que o conceito de informação aqui está sendo associado a uma espécie de

curvatura média da função de verossimilhança no sentido de que quanto maior a

curvatura mais precisa é a informação contida na verossimilhança, ou equivalen-

temente maior o valor de I(θ). Em geral espera-se que a curvatura seja negativa

e por isso seu valor é tomado com sinal trocado. Note também que a esperança

matemática é tomada em relação à distribuição amostral p(x|θ). Podemos considerar então I(θ) uma medida de informação global enquanto

que uma medida de informação local é obtida quando não se toma o valor esperado

na definição acima. A medida de informação observada de Fisher J(θ) fica então

definida como

J(θ) = −∂ 2 log p(x|θ)

∂θ2

e que será utilizada mais adiante quando falarmos sobre estimação.

Definição 2.4 Seja uma observação X com função de (densidade) de probabili-

dade p(x|θ). A priori não informativa de Jeffreys tem função de densidade dada por

p(θ) ∝ [I(θ)]1/2.

2.4. PRIORI NÃO INFORMATIVA 27

Se θ for um vetor paramétrico então p(θ) ∝ | det I(θ)|1/2.

Exemplo 2.3 : Seja X1, . . . , Xn ∼ Poisson(θ). Então o logaritmo da função de probabilidade conjunta é dado por

log p(x|θ) = −nθ + n

i=1

xi log θ − log n ∏

i=1

xi!

e tomando-se a segunda derivada segue que

∂2 log p(x|θ) ∂θ2

= ∂

∂θ

[

−n+ ∑n

i=1 xi θ

]

= − ∑n

i=1 xi θ2

e assim,

I(θ) = 1

θ2 E

[

n ∑

i=1

xi

]

= n/θ ∝ θ−1.

Portanto, a priori não informativa de Jeffreys para θ no modelo Poisson é p(θ) ∝ θ−1/2. Note que esta priori é obtida tomando-se a conjugada natural Gama(α, β)

e fazendo-se α = 1/2 e β → 0. Em geral a priori não informativa é obtida fazendo-se o parâmetro de escala

da distribuição conjugada tender a zero e fixando-se os demais parâmetros conve-

nientemente. Além disso, a priori de Jeffreys assume formas espećıficas em alguns

modelos que são frequentemente utilizados como veremos a seguir.

Definição 2.5 X tem um modelo de locação se existem uma função f e uma

quantidade θ tais que p(x|θ) = f(x − θ). Neste caso θ é chamado de parâmetro de locação.

A definição vale também quando θ é um vetor de parâmetros. Alguns exem-

plos importantes são a distribuição normal com variância conhecida, e a dis-

tribuição normal multivariada com matriz de variância-covariância conhecida.

Pode-se mostrar que para o modelo de locação a priori de Jeffreys é dada por

p(θ) ∝ constante.

Definição 2.6 X tem um modelo de escala se existem uma função f e uma

quantidade σ tais que p(x|σ) = (1/σ)f(x/σ). Neste caso σ é chamado de parâmetro de escala.

Alguns exemplos são a distribuição exponencial com parâmetro θ, com parâmetro

de escala σ = 1/θ, e a distribuição N(θ, σ2) com média conhecida e escala σ.

Pode-se mostrar que para o modelo de escala a priori de Jeffreys é dada por

p(σ) ∝ σ−1.

28 CAPÍTULO 2. DISTRIBUIÇÕES A PRIORI

Definição 2.7 X tem um modelo de locação e escala se existem uma função f

e as quantidades θ e σ tais que

p(x|θ, σ) = 1 σ f

(

x− θ σ

)

.

Neste caso θ é chamado de parâmetro de locação e σ de parâmetro de escala.

Alguns exemplos são a distribuição normal (uni e multivariada) e a distribuição

de Cauchy. Em modelos de locação e escala, a priori não informativa pode ser

obtida assumindo-se independência a priori entre θ e σ de modo que p(θ, σ) =

p(θ)p(σ) ∝ σ−1.

Exemplo 2.4 : Seja X1, . . . , Xn ∼ N(µ, σ2) com µ e σ2 desconhecidos. Neste caso,

p(x|µ, σ2) ∝ 1 σ exp

{

−1 2

(

x− µ σ

)2 }

,

portanto (µ, σ) é parâmetro de locação-escala e p(µ, σ) ∝ σ−1 é a priori não informativa. Então, pela propriedade da invariância, a priori não informativa

para (µ, σ2) no modelo normal é p(µ, σ2) ∝ σ−2.

Vale notar entretanto que a priori não informativa de Jeffreys viola o prinćı-

pio da verossimilhança, já que a informação de Fisher depende da distribuição

amostral.

2.5 Prioris Hierárquicas

A idéia aqui é dividir a especificação da distribuição a priori em estágios. Além

de facilitar a especificação esta abordagem é natural em determinadas situações

experimentais.

A distribuição a priori de θ depende dos valores dos hiperparâmetros φ e pode-

mos escrever p(θ|φ) ao invés de p(θ). Além disso, ao invés de fixar valores para os hiperparâmetros podemos especificar uma distribuição a priori p(φ) completando

assim o segundo estágio na hierarquia. Assim, a distribuição a priori conjunta é

simplesmente p(θ, φ) = p(θ|φ)p(φ) e a distribuição a priori marginal de θ pode ser então obtida por integração como

p(θ) =

p(θ, φ)dφ =

p(θ|φ)p(φ)dφ.

2.5. PRIORIS HIERÁRQUICAS 29

A distribuição a posteriori conjunta fica

p(θ, φ|x) ∝ p(x|θ, φ)p(θ|φ)p(φ) ∝ p(x|θ)p(θ|φ)p(φ)

pois a distribuição dos dados depende somente de θ. Em outras palavras, dado

θ, x e φ são independentes.

Exemplo 2.5 : Sejam X1, . . . , Xn tais que Xi ∼ N(θi, σ2) com σ2 conhecido e queremos especificar uma distribuição a priori para o vetor de parâmetros θ =

(θ1, . . . , θn). Suponha que no primeiro estágio assumimos que θi ∼ N(µ, τ 2), i = 1, . . . , n. Neste caso, se fixarmos o valor de τ 2 = τ 20 e assumirmos que µ tem

distribuição normal então θ terá distribuição normal multivariada. Por outro

lado, fixando um valor para µ = µ0 e assumindo que τ −2 tem distribuição Gama

implicará em uma distribuição t de Student multivariada para θ.

Teoricamente, não há limitação quanto ao número de estágios, mas devido às

complexidades resultantes as prioris hierárquicas são especificadas em geral em 2

ou 3 estágios. Além disso, devido à dificuldade de interpretação dos hiperparâmet-

ros em estágios mais altos é prática comum especificar prioris não informativas

para este ńıveis.

Uma aplicação interessante do conceito de hierarquia é quando a informação a

priori dispońıvel só pode ser convenientemente resumida através de uma mistura

de distribuições. Isto implica em considerar uma distribuição discreta para φ de

modo que, se φ assume os posśıveis valores φ1, . . . , φk então

p(θ) = k

i=1

p(θ|φi)p(φi).

Não é dif́ıcil verificar que a distribuição a posteriori de θ é também uma mistura

com veremos a seguir. Aplicando o teorema de Bayes temos que,

p(θ|x) = p(θ)p(x|θ)∫ p(θ)p(x|θ)dθ

=

k ∑

i=1

p(x|θ)p(θ|φi)p(φi)

k ∑

i=1

p(φi)

p(x|θ)p(θ|φi)dθ .

Mas note que a distribuição a posteriori condicional de θ dado φi é obtida via

teorema de Bayes como

p(θ|x, φi) = p(x|θ)p(θ|φi)

p(x|θ)p(θ|φi)dθ =

p(x|θ)p(θ|φi) m(x|φi)

30 CAPÍTULO 2. DISTRIBUIÇÕES A PRIORI

e a distribuição a posteriori de φi é obtida como

p(φi) = m(x|φi)p(φ)

p(x) .

Portanto p(x|θ)p(θ|φi)=p(θ|x, φi)m(x|φi). Assim, podemos escrever a posteriori de θ como

p(θ |x) =

k ∑

i=1

p(θ|x, φi)m(x|φi)p(φi)

k ∑

i=1

m(x|φi)p(φi) =

k ∑

i=1

p(θ|x, φi)p(φi|x)

Note também que p(x) = ∑

m(x|φi)p(φi), isto é a distribuição preditiva, é uma mistura de preditivas condicionais.

Exemplo 2.6 : Se θ ∈ (0, 1), a famı́lia de distribuições a priori Beta(a, b) é con- veniente. Mas estas são sempre unimodais e (se a 6= b) assimétricas à esquerda ou à direita. Outras formas interessantes, e mais de acordo com a nossa informação

a priori, podem ser obtidas misturando-se 2 ou 3 elementos desta famı́lia. Por

exemplo,

θ ∼ 0, 25Beta(3, 8) + 0, 75Beta(8, 3)

representa a informação a priori de que θ ∈ (0, 5; 0, 95) com alta probabilidade (0,71) mas também que θ ∈ (0, 1; 0, 4) com probabilidade moderada (0,20). As modas desta distribuição são 0,23 e 0,78. Por outro lado

θ ∼ 0, 33Beta(4, 10) + 0, 33Beta(15, 28) + 0, 33Beta(50, 70)

representa a informação a priori de que θ > 0, 6 com probabilidade despreźıvel.

Estas densidades estão representadas graficamente na Figura 2.2 a seguir. Note

que a primeira mistura deu origem a uma distribuição a priori bimodal enquanto

a segunda originou uma priori assimétrica à esquerda com média igual a 0,35.

Para outros exemplos de misturas de prioris ver O’Hagan (1994). Para um

excelente material sobre modelos hierárquicos ver (Gelman et al. 2004).

2.6 Problemas

1. Mostre que a famı́lia de distribuições Beta é conjugada em relação às dis-

tribuições amostrais binomial, geométrica e binomial negativa.

2.6. PROBLEMAS 31

0.0 0.2 0.4 0.6 0.8 1.0

0 1

2 3

4

θ

.33B(4,10)+.33B(15,28)+.33B(50,70)

.25 B(3,8)+.75 B(8,3)

Figura 2.2: Misturas de funções de densidade Beta(3,8) e Beta(8,3) com pesos 0,25 e 0,75 e Beta(4,10), Beta(15,28) e Beta(50,70) com pesos iguais a 0,33.

2. Para uma amostra aleatória de 100 observações da distribuição normal com

média θ e desvio-padrão 2 foi especificada uma priori normal para θ.

(a) Mostre que o desvio-padrão a posteriori será sempre menor do que 1/5.

Interprete este resultado.

(b) Se o desvio-padrão a priori for igual a 1 qual deve ser o menor número

de observações para que o desvio-padrão a posteriori seja 0,1?

3. Seja X1, . . . , Xn uma amostra aleatória da distribuição N(θ, σ 2), com θ con-

hecido. Utilizando uma distribuição a priori Gama para σ−2 com coeficiente

de variação 0,5, qual deve ser o tamanho amostral para que o coeficiente de

variação a posteriori diminua para 0,1?

4. Seja X1, . . . , Xn uma amostra aleatória da distribuição N(θ, σ 2), com θ e

σ2 desconhecidos, e considere a priori conjugada de (θ, φ).

(a) Determine os parâmetros (µ0, c0, n0, σ 2 0) utilizando as seguintes infor-

mações a priori: E(θ) = 0, P (|θ| < 1, 412) = 0, 5, E(φ) = 2 e E(φ2) = 5.

32 CAPÍTULO 2. DISTRIBUIÇÕES A PRIORI

(b) Em uma amostra de tamanho n = 10 foi observado X = 1 e ∑n

i=1(Xi − X)2 = 8. Obtenha a distribuição a posteriori de θ e es- boce os gráficos das distribuições a priori, a posteriori e da função de

verossimilhança, com φ fixo.

(c) Calcule P (|Y | > 1|x) onde Y é uma observação tomada da mesma população.

5. Suponha que o tempo, em minutos, para atendimento a clientes segue uma

distribuição exponencial com parâmetro θ desconhecido. Com base na ex-

periência anterior assume-se uma distribuição a priori Gama com média 0,2

e desvio-padrão 1 para θ.

(a) Se o tempo médio para atender uma amostra aleatória de 20 clientes

foi de 3,8 minutos, qual a distribuição a posteriori de θ.

(b) Qual o menor número de clientes que precisam ser observados para

que o coeficiente de variação a posteriori se reduza para 0,1?

6. Seja X1, . . . , Xn uma amostra aleatória da distribuição de Poisson com

parâmetro θ.

(a) Determine os parâmetros da priori conjugada de θ sabendo que E(θ) =

4 e o coeficiente de variação a priori é 0,5.

(b) Quantas observações devem ser tomadas até que a variância a poste-

riori se reduza para 0,01 ou menos?

(c) Mostre que a média a posteriori é da forma γnx + (1 − γn)µ0, onde µ0 = E(θ) e γn → 1 quando n → ∞. Interprete este resultado.

7. O número médio de defeitos por 100 metros de uma fita magnética é descon-

hecido e denotado por θ. Atribui-se uma distribuição a priori Gama(2,10)

para θ. Se um rolo de 1200 metros desta fita foi inspecionado e encontrou-se

4 defeitos qual a distribuição a posteriori de θ?

8. Seja X1, . . . , Xn uma amostra aleatória da distribuição Bernoulli com

parâmetro θ e usamos a priori conjugada Beta(a, b). Mostre que a mé-

dia a posteriori é da forma γnx + (1 − γn)µ0, onde µ0 = E(θ) e γn → 1 quando n → ∞. Interprete este resultado.

9. Para uma amostra aleatória X1, . . . , Xn tomada da distribuição U(0, θ),

mostre que a famı́lia de distribuições de Pareto com parâmetros a e b, cuja

função de densidade é p(θ) = aba/θa+1, é conjugada à uniforme.

2.6. PROBLEMAS 33

10. Para uma variável aleatória θ > 0 a famı́lia de distribuições Gama-invertida

tem função de densidade de probabilidade dada por

p(θ) = βα

Γ(α) θ−(α+1)e−β/θ, α, β > 0.

Mostre que esta famı́lia é conjugada ao modelo normal com média µ con-

hecida e variância θ desconhecida.

11. Suponha que X = (X1, X2, X3) tenha distribuição trinomial com parâmet-

ros n (conhecido) e π = (π1, π2, π3) com π1 + π2 + π3 = 1. Mostre que a

priori não informativa de Jeffreys para π é p(π) ∝ [π1π2(1− π1 − π2)]−1/2.

12. Para cada uma das distribuições abaixo verifique se o modelo é de locação,

escala ou locação-escala e obtenha a priori não informativa para os parâmet-

ros desconhecidos.

(a) Cauchy(0,β).

(b) tν(µ, σ 2), ν conhecido.

(c) Pareto(a, b), b conhecido.

(d) Uniforme (θ − 1, θ + 1). (e) Uniforme (−θ, θ).

13. Seja uma coleção de variáveis aleatórias independentes Xi com distribuições

p(xi|θi) e seja pi(θi) a priori não informativa de θi, i = 1, . . . , k. Mostre que a priori não informativa de Jeffreys para o vetor paramétrico θ = (θ1, . . . , θk)

é dada por ∏k

i=1 pi(θi).

14. Se θ tem priori não informativa p(θ) ∝ k, θ > 0 mostre que a priori de φ = aθ + b, a 6= 0 também é p(φ) ∝ k.

15. Se θ tem priori não informativa p(θ) ∝ θ−1 mostre que a priori de φ = θa, a 6= 0 também é p(φ) ∝ φ−1 e que a priori de ψ = log θ é p(ψ) ∝ k.

16. No Exemplo 1.3, sejam φi = (µi, τ 2 i ), i = 1, 2, as médias e variâncias a

priori dos f́ısicos A e B respectivamente. As prioris condicionais foram

então combinadas como

p(θ) = p(φ1)p(θ|φ1) + p(φ2)p(θ|φ2)

com p(φ1) = 0, 25 e p(φ2) = 0, 75. Usando as posterioris condicionais obti-

das naquele exemplo obtenha a distribuição a posteriori de θ (incondicional).

Esboce e comente os gráficos das densidades a priori e posteriori.

34 CAPÍTULO 2. DISTRIBUIÇÕES A PRIORI

17. Se X ∼ Binomial Negativa(v, θ) obtenha a priori de Jeffreys para θ.

18. Se X ∼ Geometrica(θ) obtenha a priori de Jeffreys para θ.

Caṕıtulo 3

Estimação

A distribuição a posteriori de um parâmetro θ contém toda a informação prob-

abiĺıstica a respeito deste parâmetro e um gráfico da sua função de densidade a

posteriori é a melhor descrição do processo de inferência. No entanto, algumas

vezes é necessário resumir a informação contida na posteriori através de alguns

poucos valores numéricos. O caso mais simples é a estimação pontual de θ onde se

resume a distribuição a posteriori através de um único número, θ̂. Como veremos

a seguir, será mais fácil entender a escolha de θ̂ no contexto de teoria da decisão.

3.1 Introdução à Teoria da Decisão

Um problema de decisão fica completamente especificado pela descrição dos

seguintes espaços:

(i) Espaço do parâmetro ou estados da natureza, Θ.

(ii) Espaço dos resultados posśıveis de um experimento, Ω.

(iii) Espaço de posśıveis ações, A.

Uma regra de decisão δ é uma função definida em Ω que assume valores em A,

i.e. δ : Ω → A. A cada decisão δ e a cada posśıvel valor do parâmetro θ podemos associar uma perda L(δ, θ) assumindo valores positivos. Definimos assim uma

função de perda.

Definição 3.1 O risco de uma regra de decisão, denotado por R(δ), é a perda

esperada a posteriori, i.e. R(δ) = Eθ|x[L(δ, θ)].

Definição 3.2 Uma regra de decisão δ∗ é ótima se tem risco mı́nimo, i.e.

R(δ∗) < R(δ), ∀δ. Esta regra será denominada regra de Bayes e seu risco, risco de Bayes.

35

36 CAPÍTULO 3. ESTIMAÇÃO

Exemplo 3.1 : Um laboratório farmaceutico deve decidir pelo lançamento ou

não de uma nova droga no mercado. É claro que o laboratório só lançará a droga

se achar que ela é eficiente mas isto é exatamente o que é desconhecido. Podemos

associar um parâmetro θ aos estados da natureza: droga é eficiente (θ = 1), droga

não é eficiente (θ = 0) e as posśıveis ações como lança a droga (δ = 1), não lança

a droga (δ = 0). Suponha que foi posśıvel construir a seguinte tabela de perdas

levando em conta a eficiência da droga,

eficiente não eficiente lança -500 600 não lança 1500 100

Vale notar que estas perdas traduzem uma avaliação subjetiva em relação à

gravidade dos erros cometidos. Suponha agora que a incerteza sobre os estados

da natureza é descrita por P (θ = 1) = π, 0 < π < 1 avaliada na distribuição

atualizada de θ (seja a priori ou a posteriori). Note que, para δ fixo, L(δ, θ) é uma

variável aleatória discreta assumindo apenas dois valores com probabilidades π e

1− π. Assim, usando a definição de risco obtemos que

R(δ = 0) = E(L(0, θ)) = π1500 + (1− π)100 = 1400π + 100 R(δ = 1) = E(L(1, θ)) = π(−500) + (1− π)600 = −1100π + 600

Uma questão que se coloca aqui é, para que valores de π a regra de Bayes será de

lançar a droga. Não é dif́ıcil verificar que as duas ações levarão ao mesmo risco,

i.e. R(δ = 0) = R(δ = 1) se somente se π = 0, 20. Além disso, para π < 0, 20

temos que R(δ = 0) < R(δ = 1) e a regra de Bayes consiste em não lançar a

droga enquanto que π > 0, 20 implica em R(δ = 1) < R(δ = 0) e a regra de Bayes

deve ser de lançar a droga.

3.2 Estimadores de Bayes

Seja agora uma amostra aleatória X1, . . . , Xn tomada de uma distribuição com

função de (densidade) de probabilidade p(x|θ) aonde o valor do parâmetro θ é desconhecido. Em um problema de inferência como este o valor de θ deve ser

estimado a partir dos valores observados na amostra.

Se θ ∈ Θ então é razoável que os posśıveis valores de um estimador δ(X) também devam pertencer ao espaço Θ. Além disso, um bom estimador é aquele

para o qual, com alta probabilidade, o erro δ(X) − θ estará próximo de zero. Para cada possivel valor de θ e cada possivel estimativa a ∈ Θ vamos associar uma perda L(a, θ) de modo que quanto maior a distância entre a e θ maior o

3.2. ESTIMADORES DE BAYES 37

valor da perda. Neste caso, a perda esperada a posteriori é dada por

E[L(a, θ)|x] = ∫

L(a, θ)p(θ|x)dθ

e a regra de Bayes consiste em escolher a estimativa que minimiza esta perda

esperada.

Aqui vamos discutir apenas funções de perda simétricas, já que estas são mais

comumente utilizadas (para outras funções de perda ver por exemplo (Bernardo

& Smith 1994) e O’Hagan 1994). Dentre estas a mais utilizada em problemas de

estimação é certamente a função de perda quadrática, definida como L(a, θ) =

(a−θ)2. Neste caso, pode-se mostrar que o estimador de Bayes para o parâmetro θ será a média de sua distribuição atualizada.

Exemplo 3.2 : Suponha que queremos estimar a proporção θ de itens defeituosos

em um grande lote. Para isto será tomada uma amostra aleatória X1, . . . , Xn de

uma distribuição de Bernoulli com parâmetro θ. Usando uma priori conjugada

Beta(α, β) sabemos que após observar a amostra a distribuição a posteriori é

Beta(α+ t, β + n− t) onde t = ∑ni=1 xi. A média desta distribuição Beta é dada por (α + t)/(α + β + n) e portanto o estimador de Bayes de θ usando perda

quadrática é

δ(X) = α +

∑n i=1 Xi

α + β + n .

A perda quadrática é as vezes criticada por penalizar demais o erro de esti-

mação. A função de perda absoluta, definida como L(a, θ) = |a − θ|, introduz punições que crescem linearmente com o erro de estimação e pode-se mostrar que

o estimador de Bayes associado é a mediana da distribuição atualizada de θ.

Para reduzir ainda mais o efeito de erros de estimação grandes podemos con-

siderar funções que associam uma perda fixa a um erro cometido, não importando

sua magnitude. Uma tal função de perda, denominada perda 0-1, é definida como

L(a, θ) =

{

1 se |a− θ| > ǫ 0 se |a− θ| < ǫ

para todo ǫ > 0. Neste caso pode-se mostrar que o estimador de Bayes é a moda

da distribuição atualizada de θ. A moda da posteriori de θ também é chamado

de estimador de máxima verossimilhança generalizado (EMVG) e é o mais fácil

de ser obtido dentre os estimadores vistos até agora. No caso cont́ınuo devemos

obter a solução da equação ∂p(θ|x)

∂θ = 0.

38 CAPÍTULO 3. ESTIMAÇÃO

Note que isto equivale a obter a solução de

∂p(x|θ)p(θ) ∂θ

= 0

e não é necessário conhecer a expressão exata de p(θ|x).

Exemplo 3.3 : Se X1, . . . , Xn é uma amostra aleatória da N(θ, σ 2) com σ2

conhecido e usarmos a priori conjugada, i.e. θ ∼ N(µ0, τ 20 ) então a posteriori também será normal e neste caso média, mediana e moda coincidem. Portanto,

o estimador de Bayes de θ é dado por

δ(X) = τ−20 µ0 + nσ

−2X

τ−20 + nσ −2

.

Exemplo 3.4 : No exemplo 3.2 suponha que foram observados 100 itens dos

quais 10 eram defeituosos. Usando perda quadrática a estimativa de Bayes de θ

δ(x) = α + 10

α + β + 100

Assim, se a priori for Beta(1,1), ou equivalentemente U(0, 1), então δ(x) = 0, 108.

Por outro lado se especificarmos uma priori Beta(1,2), que é bem diferente da an-

terior, então δ(x) = 0, 107. Ou seja, as estimativas de Bayes são bastante próxi-

mas, e isto é uma consequência do tamanho amostral ser grande. Note também

que ambas as estimativas são próximas da proporção amostral de defeituosos 0,1,

que é a estimativa de máxima verossimilhança. Se usarmos perda 0-1 e priori

Beta(1,1) então δ(x) = 0, 1.

3.3 Estimação por Intervalos

Voltamos a enfatizar que a forma mais adequada de expressar a informação que

se tem sobre um parâmetro é através de sua distribuição a posteriori. A principal

restrição da estimação pontual é que quando estimamos um parâmetro através de

um único valor numérico toda a informação presente na distribuição a posteriori

é resumida através deste número. É importante também associar alguma infor-

mação sobre o quão precisa é a especificação deste número. Para os estimadores

vistos aqui as medidas de incerteza mais usuais são a variância ou o coeficiente de

variação para a média a posteriori, a medida de informação observada de Fisher

para a moda a posteriori, e a distância entre quartis para a mediana a posteriori.

Nesta seção vamos introduzir um compromisso entre o uso da própria dis-

tribuição a posteriori e uma estimativa pontual. Será discutido o conceito de

3.4. ESTIMAÇÃO NO MODELO NORMAL 39

intervalo de credibilidade (ou intervalo de confiança Bayesiano) baseado no dis-

tribuição a posteriori.

Definição 3.3 C é um intervalo de credibilidade de 100(1-α)%, ou ńıvel de cred-

ibilidade (ou confiança) 1− α, para θ se P (θ ∈ C) ≥ 1− α.

Note que a definição expressa de forma probabiĺıstica a pertinência ou não de

θ ao intervalo. Assim, quanto menor for o tamanho do intervalo mais concentrada

é a distribuição do parâmetro, ou seja o tamanho do intervalo informa sobre a

dispersão de θ. Além disso, a exigência de que a probabilidade acima possa ser

maior do que o ńıvel de confiança é essencialmente técnica pois queremos que o

intervalo seja o menor posśıvel, o que em geral implica em usar uma igualdade.

No entanto, a desigualdade será útil se θ tiver uma distribuição discreta onde

nem sempre é posśıvel satisfazer a igualdade.

Outro fato importante é que os intervalos de credibilidade são invariantes a

transformações 1 a 1, φ(θ). Ou seja, se C = [a, b] é um intervalo de credibilidade

100(1-α)% para θ então [φ(a), φ(b)] é um intervalo de credibilidade 100(1-α)%

para φ(θ). Note que esta propriedade também vale para intervalos de confiança

na inferência clássica.

É posśıvel construir uma infinidade de intervalos usando a definição acima mas

estamos interessados apenas naquele com o menor comprimento posśıvel. Pode-se

mostrar que intervalos de comprimento mı́nimo são obtidos tomando-se os valores

de θ com maior densidade a posteriori, e esta idéia é expressa matematicamente

na definição abaixo.

Definição 3.4 Um intervalo de credibilidade C de 100(1-α)% para θ é de máx-

ima densidade a posteriori (MDP) se C = {θ ∈ Θ : p(θ|x) ≥ k(α)} onde k(α) é a maior constante tal que P (θ ∈ C) ≥ 1− α.

Usando esta definição, todos os pontos dentro do intervalo MDP terão den-

sidade maior do que qualquer ponto fora do intervalo. Além disso, no caso de

distribuições com duas caudas, e.g. normal, t de Student, o intervalo MDP é

obtido de modo que as caudas tenham a mesma probabilidade. Um problema

com os intervalos MDP é que eles não são invariantes a transformações 1 a 1, a

não ser para transformações lineares. O mesmo problema ocorre com intervalos

de comprimento mı́nimo na inferência clássica.

3.4 Estimação no Modelo Normal

Os resultados desenvolvidos nos caṕıtulos anteriores serão aplicados ao modelo

normal para estimação da média e variância em problemas de uma ou mais

40 CAPÍTULO 3. ESTIMAÇÃO

amostras e em modelos de regressão linear. A análise será feita com priori con-

jugada e priori não informativa quando serão apontadas as semelhanças com a

análise clássica. Assim como nos caṕıtulos anteriores a abordagem aqui é in-

trodutória. Um tratamento mais completo do enfoque Bayesiano em modelos

lineares pode ser encontrado em Broemeling (1985) e Box & Tiao (1992).

Nesta seção considere uma amostra aleatória X1, · · · , Xn tomada da dis- tribuição N(θ, σ2).

3.4.1 Variância Conhecida

Se σ2 é conhecido e a priori de θ é N(µ0, τ 2 0 ) então, pelo Teorema 1.1, a posteriori

de θ é N(µ1, τ 2 1 ). Intervalos de confiança Bayesianos para θ podem então ser

constrúıdos usando o fato de que

θ − µ1 τ1

|x ∼ N(0, 1).

Assim, usando uma tabela da distribuição normal padronizada podemos obter o

valor do percentil zα/2 tal que

P

(

−zα/2 ≤ θ − µ1 τ1

≤ zα/2 )

= 1− α

e após isolar θ, obtemos que

P (

µ1 − zα/2τ1 ≤ θ ≤ µ1 + zα/2τ1 )

= 1− α.

Portanto (

µ1 − zα/2τ1;µ1 + zα/2τ1 )

é o intervalo de confiança 100(1-α)% MDP

para θ, devido à simetria da normal.

A priori não informativa pode ser obtida fazendo-se a variância da priori

tender a infinito, i.e. τ 20 → ∞. Neste caso, é fácil verificar que τ−21 → nσ−2 e µ1 → x, i.e. a média e a precisão da posteriori convergem para a média e a precisão amostrais. Média, moda e mediana a posteriori coincidem então com

a estimativa clássica de máxima verossimilhança, x. O intervalo de confiança

Bayesiano 100(1-α)% é dado por

(

x− zα/2 σ/ √ n; x+ zα/2 σ/

√ n )

e também coincide numericamente com o intervalo de confiança clássico. Aqui

entretanto a interpretação do intervalo é como uma afirmação probabiĺıstica sobre

θ.

3.4. ESTIMAÇÃO NO MODELO NORMAL 41

3.4.2 Média e Variância desconhecidas

Neste caso, usando a priori conjugada Normal-Gama vista no Caṕıtulo 2 temos

que a distribuição a posteriori marginal de θ é dada por

θ|x ∼ tn1(µ1, σ21/c1).

Portanto, média, moda e mediana a posteriori coincidem e são dadas por µ1.

Denotando por tα/2,n1 o percentil 100(1-α/2)% da distribuição tn1(0, 1) podemos

obter este percentil tal que

P

(

−tα/2,n1 ≤ √ c1 θ − µ1 σ1

≤ tα/2,n1 )

= 1− α

e após isolar θ, usando a simetria da distribuição t-Student obtemos que

(

µ1 − tα/2,n1 σ1√ c1

≤ θ ≤ µ1 + tα/2,n1 σ1√ c1

)

é o intervalo de confiança Bayesiano 100(1-α)% de MDP para θ.

No caso da variância populacional σ2 intervalos de confiança podem ser obti-

dos usando os percentis da distribuição qui-quadrado uma vez que a distribuição

a posteriori de φ é tal que n1σ 2 1φ|x ∼ χ2n1 . Denotando por

χ2 α/2,n1

e χ2α/2,n1

os percentis α/2 e 1−α/2 da distribuição qui-quadrado com n1 graus de liberdade respectivamente, podemos obter estes percentis tais que

P

(

χ2 α/2,n1

n1σ21 ≤ φ ≤

χ2α/2,n1 n1σ21

)

= 1− α.

Note que este intervalo não é de MDP já que a distribuição qui-quadrado não é

simétrica. Como σ2 = 1/φ é uma função 1 a 1 podemos usar a propriedade de

invariância e portanto (

n1σ 2 1

χ2α/2,n1 ; n1σ

2 1

χ2 α/2,n1

)

é o intervalo de confiança Bayesiano 100(1-α)% para σ2.

Um caso particular é quanto utilizamos uma priori não informativa. Vimos

na Seção 2.4 que a priori não informativa de locação e escala é p(θ, σ) ∝ 1/σ, portanto pela propriedade de invariância segue que a priori não informativa de

(θ, φ) é obtida fazendo-se p(θ, φ) ∝ φ−1 pois p(θ, σ2) ∝ σ−2. Note que este é um caso particular (degenerado) da priori conjugada natural com c0 = 0, σ

2 0 = 0 e

42 CAPÍTULO 3. ESTIMAÇÃO

n0 = −1. Neste caso a distribuição a posteriori marginal de θ fica

θ|x ∼ tn−1(x, s2/n)

sendo s2 = 1/(n − 1)∑ni=1(xi − x)2 a variância amostral. Mais uma vez média, moda e mediana a posteriori de θ coincidem com a média amostral x que é a

estimativa de máxima verossimilhança. Como √ n(θ − x)/s ∼ tn−1(0, 1) segue

que o intervalo de confiança 100(1-α)% para θ de MDP é

(

x− tα/2,n−1 s√ n ; x+ tα/2,n−1

s√ n

)

que coincide numericamente com o intervalo de confiança clássico.

Para fazer inferências sobre σ2 temos que

φ|x ∼ Gama (

n− 1 2

, (n− 1)s2

2

)

ou (n− 1)s2φ|x ∼ χ2n−1.

A estimativa pontual de σ2 utilizada é [E(φ|x)]−1 = s2 que coincide com a estimativa clássica uma vez que o estimador de máxima verossimilhança

(n − 1)S2/n é viciado e normalmente substituido por S2 (que é não viciado). Os intervalos de confiança 100(1-α)% Bayesiano e clássico também coincidem e

são dados por (

(n− 1)s2 χ2α/2,n−1

; (n− 1)s2 χ2 α/2,n−1

)

.

Mais uma vez vale enfatizar que esta coincidência com as estimativas clás-

sicas é apenas numérica uma vez que as interpretações dos intervalos diferem

radicalmente.

3.4.3 O Caso de duas Amostras

Nesta seção vamos assumir que X11, . . . , X1n1 e X21, . . . , X2n2 são amostras

aleatórias das distribuições N(θ1, σ 2 1) e N(θ2, σ

2 2) respectivamente e que as

amostras são independentes.

Para começar vamos assumir que as variâncias σ21 e σ 2 2 são conhecidas. Neste

caso, a função de verossimilhança é dada por

p(x1,x2|θ1, θ2) = p(x1|θ1)p(x2|θ2)

∝ exp {

− n1 2σ21

(θ1 − x1)2 }

exp

{

− n2 2σ22

(θ2 − x2)2 }

isto é, o produto de verossimilhanças relativas a θ1 e θ2. Assim, se assumirmos

que θ1 e θ2 são independentes a priori então eles também serão independentes a

3.4. ESTIMAÇÃO NO MODELO NORMAL 43

posteriori já que

p(θ1, θ2|x1,x2) = p(x1|θ1)p(θ1)

p(x1) × p(x2|θ2)p(θ2)

p(x2) .

Se usarmos a classe de prioris conjugadas θi ∼ N(µi, τ 2i ) então as posterioris independentes serão θi|xi ∼ N(µ∗i , τ ∗

2

i ) onde

µ∗i = τ−2i µi + niσ

−2 i xi

τ−2i + niσ −2 i

e τ ∗ 2

i = 1/(τ −2 i + niσ

−2 i ), i = 1, 2.

Em geral estaremos interessados em comparar as médias populacionais, i.e

queremos estimar β = θ1 − θ2 (por exemplo, testar se θ1 = θ2). Neste caso, a posteriori de β é facilmente obtida, devido à independência, como

β|x1,x2 ∼ N(µ∗1 − µ∗2, τ ∗ 2

1 + τ ∗2

2 )

e podemos usar µ∗1 − µ∗2 como estimativa pontual para a diferença e também construir um intervalo de credibilidade MDP para esta diferença.

(µ∗1 − µ∗2)± zα/2 √

τ ∗ 2

1 + τ ∗2 2 .

Note que se usarmos priori não informativa, i.e. fazendo τ 2i → ∞, i = 1, 2 então a posteriori fica

β|x1,x2 ∼ N (

x1 − x2, σ21 n1

+ σ22 n2

)

e o intervalo obtido coincidirá mais uma vez com o intervalo de confiança clássico.

No caso de variâncias populacionais desconhecidas porém iguais, temos que

φ = σ−21 = σ −2 2 = σ

−2. A priori conjugada pode ser constrúıda em duas etapas.

No primeiro estágio, assumimos que, dado φ, θ1 e θ2 são a priori condicionalmente

independentes, e especificamos

θi|φ ∼ N(µi, (ciφ)−1), i = 1, 2.

e no segundo estágio, especificamos a priori conjugada natural para φ, i.e.

φ ∼ Gama (

n0 2 , n0σ

2 0

2

)

.

Combinando as prioris acima não é dif́ıcil verificar que a priori conjunta de

44 CAPÍTULO 3. ESTIMAÇÃO

(θ1, θ2, φ) é

p(θ1, θ2, φ) = p(θ1|φ)p(θ2|φ)p(φ)

∝ φn0/2 exp {

−φ 2

[

n0σ 2 0 + c1(θ1 − µ1)2 + c2(θ2 − µ2)2

]}

.

Além disso, também não é dif́ıcil obter a priori condicional de β = θ1 − θ2, dado φ, como

β|φ ∼ N(µ1 − µ2, φ−1(c−11 + c−12 ))

e portanto, usando os resultados da Seção 2.3.5 segue que a distribuição a priori

marginal da diferença é

β ∼ tn0(µ1 − µ2, σ20(c−11 + c−12 )).

Podemos mais uma vez obter a posteriori conjunta em duas etapas já que θ1 e

θ2 também serão condicionalmente independentes a posteriori, dado φ. Assim, no

primeiro estágio usando os resultados obtidos anteriormente para uma amostra

segue que

θi|φ,x ∼ N(µ∗i , (c∗iφ)−1), i = 1, 2

onde

µ∗i = ciµi + nixi ci + ni

e c∗i = ci + ni.

Na segunda etapa temos que combinar a verossimilhança com a priori de

(θ1, θ2, φ). Definindo a variância amostral combinada

s2p = (n1 − 1)S21 + (n2 − 1)S22

n1 + n2 − 2

e denotando ν = n1 + n2 − 2, a função de verossimilhança pode ser escrita como

p(x1,x2|θ1, θ2, φ) = φ(n1+n2)/2 exp {

−φ 2

[

νs2 + n1(θ1 − x1)2 + n2(θ2 − x2)2 ]}

e após algum algebrismo obtemos que a posteriori é proporcional a

φ(n0+n1+n2)/2 exp

{

−φ 2

[

n0σ 2 0 + νs

2 + 2

i=1

cini c∗i

(µi − xi)2 + c∗i (θi − µ∗i )2 ]

}

.

Como esta posteriori tem o mesmo formato da priori segue por analogia que

φ|x ∼ Gama (

n∗0 2 , n∗0σ

∗2

0

2

)

3.4. ESTIMAÇÃO NO MODELO NORMAL 45

onde n∗0 = n0+n1+n2 e n ∗ 0σ

∗2

0 = n0σ 2 0 + νs

2+ ∑2

i=1 cini(µi−xi)2/c∗i . Ainda por analogia com o caso de uma amostra, a posteriori marginal da diferença é dada

por

β|x ∼ tn∗ 0 (µ∗1 − µ∗2, σ∗

2

0 (c ∗−1

1 + c ∗−1

2 )).

Assim, média, moda e mediana a posteriori de β coincidem e a estimativa

pontual é µ∗1−µ∗2. Também intervalos de credibilidade de MDP podem ser obtidos usando os percentis da distribuição t de Student. Para a variância populacional

a estimativa pontual usual é σ∗ 2

0 e intervalos podem ser constrúıdos usando os

percentis da distribuição qui-quadrado já que n∗0σ ∗2

0 φ | x ∼ χ2n∗ 0

Vejamos agora como fica a análise usando priori não informativa. Neste caso,

p(θ1, θ2, φ) ∝ φ−1 e isto equivale a um caso particular (degenerado) da priori conjugada com ci = 0, σ

2 0 = 0 e n0 = −2. Assim, temos que c∗i = ni, µ∗i = xi,

n∗0 = ν e n ∗ 0σ

∗2

0 = νs 2 e a estimativa pontual concide com a estimativa de máxima

verossimilhança β̂ = x1 − x2. O intervalo de 100(1 − α)% de MDP para β tem limites

x1 − x2 ± tα 2 ,ν sp

1

n1 +

1

n2

que coincide numericamente com o intervalo de confiança clássico.

O intervalo de 100(1 − α)% para σ2 é obtido de maneira análoga ao caso de uma amostra usando a distribuição qui-quadrado, agora com ν graus de liberdade,

i.e. (

νs2p χ2α

2 ,ν

, νs2p χ2α

2 ,ν

)

.

3.4.4 Variâncias desiguais

Até agora assumimos que as variâncias populacionais desconhecidas eram iguais

(ou pelo menos aproximadamente iguais). Na inferência clássica a violação desta

suposição leva a problemas teóricos e práticos uma vez que não é trivial encontrar

uma quantidade pivotal para β com distribuição conhecida ou tabelada. Na

verdade, se existem grandes diferenças de variabilidade entre as duas populações

pode ser mais apropriado analisar conjuntamente as consequências das diferenças

entre as médias e as variâncias. Assim, caso o pesquisador tenha interesse no

parâmetro β deve levar em conta os problemas de ordem teórica introduzidos por

uma diferença substancial entre σ21 e σ 2 2.

Do ponto de vista Bayesiano o que precisamos fazer é combinar informação a

priori com a verossimilhança e basear a estimação na distribuição a posteriori. A

função de verossimilhança agora pode ser fatorada como

p(x1,x2|θ1, θ2, σ21, σ22) = p(x1|θ1, σ21)p(x2|θ2, σ22)

46 CAPÍTULO 3. ESTIMAÇÃO

e vamos adotar prioris conjugadas normal-gama independentes com parâmetros

(µi, ci, νi, σ 2 0i) para cada uma das amostras. Fazendo as operações usuais para

cada amostra, e usando a conjugação da normal-gama, obtemos as seguintes

distribuições a posteriori independentes

θi|x ∼ tn∗ 0i (µ∗i , σ

∗2

0i /c ∗ i ) e φi|x ∼ Gama

(

n∗0i 2 , n∗0iσ

∗2

0i

2

)

, i = 1, 2.

Pode-se mostrar que β tem uma distribuição a posteriori chamada Behrens-

Fisher, que é semelhante à t de Student e é tabelada. Assim, intervalos de

credibilidade podem ser constrúıdos usando-se estes valores tabelados.

Outra situação de interesse é a comparação das duas variâncias populacionais.

Neste caso, faz mais sentido utilizar a razão de variâncias ao invés da diferença

já que elas medem a escala de uma distribuição e são sempre positivas. Neste

caso temos que obter a distribuição a posteriori de σ22/σ 2 1 = φ1/φ2. Usando a

independência a posteriori de φ1 e φ2 e após algum algebrismo pode-se mostrar

que σ∗

2

01

σ∗ 2

02

φ1 φ2

∼ F (n∗01, n∗02).

Embora sua função de distribuição não possa ser obtida analiticamente os val-

ores estão tabelados em muitos livros de estat́ıstica e também podem ser obtidos

na maioria dos pacotes computacionais. Os percentis podem então ser utilizados

na construção de intervalos de credibilidade para a razão de variâncias.

Uma propriedade bastante útil para calcular probabilidade com a distribuição

F vem do fato de que se X ∼ F (ν2, ν1) então X−1 ∼ F (ν1, ν2) por simples inver- são na razão de distribuições qui-quadrado independentes. Assim, denotando os

quantis α e 1−α da distribuição F (ν1, ν2) por Fα(ν1, ν2) e F α(ν1, ν2) respectiva- mente segue que

F α(ν1, ν2) = 1

F α(ν2, ν1) .

Note que é usual que os livros forneçam tabelas com os percentis superiores da

distribuição F para várias combinações de valores de ν1 e ν2 devido à propriedade

acima. Por exemplo, se temos os valores tabelados dos quantis 0,95 podemos obter

também um quantil 0,05. Basta procurar o quantil 0,95 inverterndo os graus de

liberdade.

Finalmente, a análise usando priori não informativa pode ser feita para

p(θ1, θ2, σ 2 1, σ

2 2) ∝ σ−21 σ−22 e será deixada como exerćıcio.

3.5. EXERCÍCIOS 47

3.5 Exerćıcios

1. Gere 2 amostras de tamanho 50 da distribuição N(0, 1). Agora construa um

intervalo MDP de 95% para a diferença entre as médias (assuma variância

conhecida igual a 1). Qual a sua conclusão?

2. Repita a análise da Seção 3.4.4 usando priori não informativa para

p(θ1, θ2, σ 2 1, σ

2 2) ∝ σ−21 σ−22 .

Caṕıtulo 4

Métodos Aproximados

4.1 Computação Bayesiana

Existem várias formas de resumir a informação descrita na distribuição a poste-

riori. Esta etapa frequentemente envolve a avaliação de probabilidades ou esper-

anças.

Neste caṕıtulo serão descritos métodos baseados em simulação, incluindo

Monte Carlo simples, Monte Carlo com função de importância, métodos de

reamostragem e Monte Carlo via cadeias de Markov (MCMC). O material apre-

sentado é introdutório e mais detalhes sobre os estes métodos podem ser obtidos

por exemplo em Gamerman (1997), Robert & Casella (1999) e Gamerman &

Lopes (2006). Outros métodos computacionalmente intensivos como técnicas de

otimização e integração numérica, bem como aproximações anaĺıticas não serão

tratados aqui e uma referência introdutória é Migon & Gamerman (1999).

Todos os algoritmos que serão vistos aqui são não determińısticos, i.e. todos

requerem a simulação de números (pseudo) aleatórios de alguma distribuição de

probabilidades. Em geral, a única limitação para o número de simulações são o

tempo de computação e a capacidade de armazenamento dos valores simulados.

Assim, se houver qualquer suspeita de que o número de simulações é insuficiente,

a abordagem mais simples consiste em simular mais valores.

4.2 Uma Palavra de Cautela

Apesar da sua grande utilidade, os métodos que serão apresentados aqui devem ser

aplicados com cautela. Devido à facilidade com que os recursos computacionais

podem ser utilizados hoje em dia, corremos o risco de apresentar uma solução para

o problema errado (o erro tipo 3) ou uma solução ruim para o problema certo.

Assim, os métodos computacionalmente intensivos não devem ser vistos como

substitutos do pensamento cŕıtico sobre o problema por parte do pesquisador.

48

4.3. O PROBLEMA GERAL DA INFERÊNCIA BAYESIANA 49

Além disso, sempre que posśıvel deve-se utilizar soluções exatas, i.e. não

aproximadas, se elas existirem. Por exemplo, em muitas situações em que pre-

cisamos calcular uma integral múltipla existe solução exata em algumas dimen-

sões, enquanto nas outras dimensões temos que usar métodos de aproximação.

4.3 O Problema Geral da Inferência Bayesiana

A distribuição a posteriori pode ser convenientemente resumida em termos de

esperanças de funções particulares do parâmetro θ, i.e.

E[g(θ)|x] = ∫

g(θ)p(θ|x)dθ

ou distribuições a posteriori marginais quando θ for multidimensional, por exem-

plo se θ = (θ1,θ2) então

p(θ1|x) = ∫

p(θ|x)dθ2.

Assim, o problema geral da inferência Bayesiana consiste em calcular tais

valores esperados segundo a distribuição a posteriori de θ. Alguns exemplos são,

1. Constante normalizadora. g(θ) = 1 e p(θ|x) = kq(θ), segue que

k =

[∫

q(θ)dθ

]−1

.

2. Se g(θ) = θ, então têm-se µ = E(θ|x), média a posteriori.

3. Quando g(θ) = (θ−µ)2, então σ2 = E((θ−µ)2|x), a variância a posteriori.

4. Se g(θ) = IA(θ), onde IA(x) = 1 se x ∈ A e zero caso contrário, então P (A | x) =

A

p(θ|x)dθ

5. Seja g(θ) = p(y|θ), onde y ⊥ x|θ. Nestas condições obtemos E[p(y|x)], a distribuição preditiva de y, uma observação futura.

Portanto, a habilidade de integrar funções, muitas vezes complexas e multi-

dimensionais, é extremamente importante em inferência Bayesiana. Inferência

exata somente será posśıvel se estas integrais puderem ser calculadas analitica-

mente, caso contrário devemos usar aproximações. Nas próximas seções iremos

apresentar métodos aproximados baseados em simulação para obtenção dessas

integrais.

50 CAPÍTULO 4. MÉTODOS APROXIMADOS

4.4 Método de Monte Carlo Simples

A idéia do método é justamente escrever a integral que se deseja calcular como

um valor esperado. Para introduzir o método considere o problema de calcular a

integral de uma função g(θ) no intervalo (a, b), i.e.

I =

∫ b

a

g(θ)dθ.

Esta integral pode ser reescrita como

I =

∫ b

a

(b− a)g(θ) 1 b− adθ = (b− a)E[g(θ)]

identificando θ como uma variável aleatória com distribuição U(a, b). Assim,

transformamos o problema de avaliar a integral no problema estat́ıstico de es-

timar uma média, E[g(θ)]. Se dispomos de uma amostra aleatória de tamanho

n, θ1, . . . , θn da distribuição uniforme no intervalo (a, b) teremos também uma

amostra de valores g(θ1), . . . , g(θn) da função g(θ) e a integral acima pode ser

estimada pela média amostral, i.e.

Î = (b− a) 1 n

n ∑

i=1

g(θi).

Não é dif́ıcil verificar que esta estimativa é não viesada já que

E(Î) = (b− a)

n

n ∑

i=1

E[g(θi)] = (b− a)E[g(θ)] = ∫ b

a

g(θ)dθ.

Podemos então usar o seguinte algoritmo

1. gere θ1, . . . , θn da distribuição U(a, b);

2. calcule g(θ1), . . . , g(θn);

3. calcule a média amostral g = ∑n

i=1 g(θi)/n

4. calcule Î = (b− a)g

Exemplo 4.1 : Suponha que queremos calcular ∫ 3

1 exp(−x)dx. A integral pode

ser reescrita como

(3− 1) ∫ 3

1

exp(−x)/(3− 1)dx

e será aproximada usando 100 valores simulados da distribuição Uniforme no

intervalo (1,3) e calculando yi = e −xi , i = 1, . . . , 100. O valor aproximado da

4.4. MÉTODO DE MONTE CARLO SIMPLES 51

integral é 2 ∑100

i=1 yi/100. Por outro lado, sabemos que exp(−x) é a função de densidade de uma v.a. X ∼ Exp(1) e portanto a integral pode ser calculada de forma exata,

∫ 3

1

exp(−x)dx = Pr(X < 3)− Pr(X < 1) = 0.3181.

Podemos escrever uma função mais geral no R cujos argumentos são o número

de simulações e os limites de integração.

> int.exp = function(n, a, b) {

+ x = runif(n, a, b)

+ y = exp(-x)

+ int.exp = (b - a) * mean(y)

+ return(int.exp)

+ }

Executando a função int.exp digamos 50 vezes com n = 10, a = 1 e b = 3

existirá uma variação considerável na estimativa da integral. Veja a Figura 4.1.

Isto se chama “erro de Monte Carlo” e decresce conforme aumentamos o número

de simulações. Repetindo o experimento com n = 1000 a variação ficará bem

menor. Na Figura 4.2 a evolução deste erro conforme se aumenta o número de

simulações fica bem evidente. Os comandos do R a seguir foram utilizados.

> n = c(20, 50, 100, 200, 500)

> y = matrix(0, ncol = length(n), nrow = 50)

> for (j in 1:length(n)) {

+ m = NULL

+ for (i in 1:50) m = c(m, int.exp(n[j], 1, 3))

+ y[, j] = m

+ }

> boxplot(data.frame(y), names = n)

A generalização é bem simples para o caso em que a integral é a esperança

matemática de uma função g(θ) onde θ tem função de densidade p(θ), i.e.

I =

∫ b

a

g(θ)p(θ)dθ = E[g(θ)]. (4.1)

Neste caso, podemos usar o mesmo algoritmo descrito acima modificando o passo

1 para gerar θ1, . . . , θn da distribuição p(θ) e calculando

Î = g = 1

n

n ∑

i=1

g(θi).

52 CAPÍTULO 4. MÉTODOS APROXIMADOS

0.20 0.25 0.30 0.35 0.40

0 2

4 6

8

Figura 4.1: Histograma de 50 estimativas de Monte Carlo da integral no Exemplo 4.1 com n = 10.

Uma vez que as gerações são independentes, pela Lei Forte dos Grandes

Números segue que Î converge quase certamente para I,

1

n

n ∑

i=1

g(θi) → E[g(θ], n → ∞.

Além disso, temos uma amostra g(θ1), . . . , g(θn) tal que

E[g(θi)] = E[g(θ)] = I e V ar[g(θi)] = σ 2 =

1

n

(g(θi)− ḡ)2

e portanto a variância do estimador pode também ser estimada como

v = 1

n2

n ∑

i=1

(g(θi)− g)2,

i.e. a aproximação pode ser tão acurada quanto se deseje bastando aumentar o

valor de n. É importante notar que n está sob nosso controle aqui, e não se trata

do tamanho da amostra de dados.

O Teorema Central do Limite também se aplica aqui de modo que para n

4.4. MÉTODO DE MONTE CARLO SIMPLES 53

20 50 100 200 500

0. 20

0. 25

0. 30

0. 35

0. 40

Figura 4.2: Boxplots para 50 estimativas da integral no Exemplo 4.1 com n=20, 50, 100, 200, e 500 simulações.

grande segue que g − E[g(θ)]√

v

tem distribuição aproximadamente N(0, 1). Podemos usar este resultado para

testar convergência e construir intervalos de confiança.

No caso multivariado a extensão também é direta. Seja θ = (θ1, . . . , θk) ′

um vetor aleatório de dimensão k com função de densidade p(θ). Neste caso os

valores gerados serão também vetores θ1, . . . ,θn e o estimador de Monte Carlo

fica

Î = 1

n

n ∑

i=1

g(θi)

Exemplo 4.2 : Suponha que queremos calcular Pr(X < 1, Y < 1) onde o ve-

tor aleatório (X, Y ) tem distribuição Normal padrão bivariada com correlação

igual a 0,5. Note que esta probabilidade é a integral de p(x, y) definida no inter-

valo acima, portanto simulando valores desta distribuição poderemos estimar esta

probabilidade como a proporção de pontos que caem neste intervalo. A Figura 4.3

apresenta um diagrama de dispersão dos valores simulados e foi obtida usando os

camandos do R abaixo.

54 CAPÍTULO 4. MÉTODOS APROXIMADOS

> Sigma = matrix(c(1, 0.5, 0.5, 1), 2, 2)

> m = c(0, 0)

> require(MASS)

> y = mvrnorm(n = 1000, mu = m, Sigma = Sigma)

> plot(y[, 1], y[, 2], xlab = "x", ylab = "y")

> abline(1, 0)

> abline(v = 1)

−3 −2 −1 0 1 2 3

− 3

− 2

− 1

0 1

2 3

x

y

Figura 4.3: Diagrama de dispersão de 1000 valores simulados da distribuição N(0,1) bivariada.

Uma grande vantagem dos métodos de simulação é que após uma amostra

de vetores aleatórios ser gerada podemos facilmente calcular caracteŕısticas das

distribuições marginais e condicionais. No Exemplo 4.2, para calcular Pr(X < 1)

basta calcular a frequência relativa de pontos (xi, yi) tais que xi < 1. Para

calcular a probabilidade condicional Pr(X < 1|Y < 1) basta selecionar somente aqueles pontos cuja segunda coordenada é menor do que 1. Depois calcula-se a

frequência relativa dos pontos restantes cuja primeira coordenada é menor do que

1.

4.4.1 Monte Carlo via Função de Importância

Em muitas situações pode ser muito custoso ou mesmo imposśıvel simular valores

da distribuição a posteriori. Neste caso, pode-se recorrer à uma função q(θ) que

seja de fácil amostragem, usualmente chamada de função de importância. O

procedimento é comumente chamado de amostragem por importância.

4.4. MÉTODO DE MONTE CARLO SIMPLES 55

−4 −2 0 2 4

0. 0

0. 1

0. 2

0. 3

0. 4

x

p( x)

−4 −2 0 2 4

0. 0

0. 1

0. 2

0. 3

0. 4

y

p( y)

−4 −2 0 2

0. 0

0. 1

0. 2

0. 3

0. 4

p( x

| y <

1)

−4 −2 0 2

0. 0

0. 1

0. 2

0. 3

0. 4

p( y

| x <

1)

Figura 4.4: Estimativas das densidades marginais e condicionais no Exemplo 4.2.

Se q(θ) for uma função de densidade definida no mesmo espaço variação de θ

então a integral (4.1) pode ser reescrita como

I =

g(θ)p(θ)

q(θ) q(θ)dx = E

[

g(θ)p(θ)

q(θ)

]

onde a esperança agora é com respeito a distribuição q. Assim, se dispomos de

uma amostra aleatória θ1, . . . , θn tomada da distribuição q o estimador de Monte

Carlo da integral acima fica

Î = 1

n

n ∑

i=1

g(θi)p(θi)

q(θi) .

e tem as mesmas propriedades do estimador de Monte Carlo simples.

Em prinćıpio não há restrições quanto à escolha da densidade de importância

q, porém na prática alguns cuidados devem ser tomados. Pode-se mostrar que

a escolha ótima no sentido de minimizar a variância do estimador consiste em

tomar q(θ) ∝ g(θ)p(θ).

56 CAPÍTULO 4. MÉTODOS APROXIMADOS

Exemplo 4.3 : Para uma única observação X suponha que

X|θ ∼ N(θ, 1) e θ ∼ Cauchy(0, 1).

Então,

p(x|θ) ∝ exp[−(x− θ)2/2] e p(θ) = 1 π(1 + θ2)

.

Portanto, a densidade a posteriori de θ é dada por

p(θ|x) = 1

1 + θ2 exp[−(x− θ)2/2]

1

1 + θ2 exp[−(x− θ)2/2]dθ

.

Suponha agora que queremos estimar θ usando função de perda quadrática. Como

vimos no Caṕıtulo 3 isto implica em tomar a média a posteriori de θ como esti-

mativa. Mas

E[θ|x] = ∫

θp(θ|x)dθ =

θ

1 + θ2 exp[−(x− θ)2/2]dθ

1

1 + θ2 exp[−(x− θ)2/2]dθ

e as integrais no numerador e denominador não têm solução anaĺıtica exata.

Uma solução aproximada via simulação de Monte Carlo pode ser obtida usando

o seguinte algoritmo,

1. gerar θ1, . . . , θn independentes da distribuição N(x, 1);

2. calcular gi = θi

1 + θ2i e g∗i =

1

1 + θ2i ;

3. calcular Ê(θ|x) = ∑n

i=1 gi ∑n

i=1 g ∗ i

.

Este exemplo ilustrou um problema que geralmente ocorre em aplicações

Bayesianas. Como a posteriori só é conhecida a menos de uma constante de

proporcionalidade as esperanças a posteriori são na verdade uma razão de inte-

grais. Neste caso, a aproximação é baseada na razão dos dois estimadores de

Monte Carlo para o numerador e denominador.

Exercicios

1. Para cada uma das distribuições N(0, 1), Gama(2,5) e Beta(2,5) gere 100,

1000 e 5000 valores independentes. Faça um gráfico com o histograma e

4.5. MÉTODOS DE REAMOSTRAGEM 57

a função de densidade superimposta. Estime a média e a variância da

distribuição. Estime a variância do estimador da média.

2. Para uma única observação X com distribuição N(θ, 1), θ desconhecido,

queremos fazer inferência sobre θ usando uma priori Cauchy(0,1). Gere um

valor de X para θ = 2, i.e. x ∼ N(2, 1).

(a) Estime θ através da sua média a posteriori usando o algoritmo do

Exemplo 4.3.

(b) Estime a variância da posteriori.

(c) Generalize o algoritmo para k observações X1, . . . , Xk da distribuição

N(θ, 1).

4.5 Métodos de Reamostragem

Existem distribuições para as quais é muito dif́ıcil ou mesmo imposśıvel simular

valores. A idéia dos métodos de reamostragem é gerar valores em duas etapas.

Na primeira etapa gera-se valores de uma distribuição auxiliar conhecida. Na

segunda etapa utiliza-se um mecanismo de correção para que os valores sejam

representativos (ao menos aproximadamente) da distribuição a posteriori. Na

prática costuma-se tomar a priori como distribuição auxiliar conforme proposto

em Smith & Gelfand (1992).

4.5.1 Método de Rejeição

Considere uma função de densidade auxiliar q(θ) da qual sabemos gerar valores.

A única restrição é que exista uma constante A finita tal que p(θ|x) < Aq(θ). O método de rejeição consiste em gerar um valor θ∗ da distribuição auxiliar q

e aceitar este valor como sendo da distribuição a posteriori com probabilidade

p(θ∗|x)/Aq(θ∗). Caso contrário, θ∗ não é aceito como um valor gerado da pos- teriori e o processo é repetido até que um valor seja aceito. O método também

funciona se ao invés da posteriori, que em geral é desconhecida, usarmos a sua

versão não normalizada, i.e p(x|θ)p(θ). Podemos então usar o seguinte algoritmo,

1. gerar um valor θ∗ da distribuição auxiliar;

2. gerar u ∼ U(0, 1);

3. se u < p(θ∗|x)/Aq(θ∗) faça θ(j) = θ∗, faça j = j + 1 e retorne ao passo 1. caso contrário retorne ao passo 1.

58 CAPÍTULO 4. MÉTODOS APROXIMADOS

Tomando a priori p(θ) como densidade auxiliar a constante A deve ser tal que

p(x|θ) < A. Esta desigualdade é satisfeita se tomarmos A como sendo o valor máximo da função de verossimilhança, i.e. A = p(x|θ̂) onde θ̂ é o estimador de máxima verossimilhança de θ. Neste caso, a probabilidade de aceitação se

simplifica para p(x|θ)/p(x|θ̂). Podemos então usar o seguinte algoritmo para gerar valores da posteriori,

1. gerar um valor θ∗ da distribuição a priori;

2. gerar u ∼ U(0, 1);

3. aceitar θ∗ como um valor da posteriori se u < p(x|θ∗)/p(x|θ̂), caso contrário rejeitar θ∗ e retornar ao passo 1.

Exemplo 4.4 : Suponha que X1, . . . , Xn ∼ N(θ, 1) e assume-se uma distribuição a priori Cauchy(0,1) para θ. A função de verossimilhança é,

p(x|θ) = n ∏

i=1

(2π)−1/2 exp

{

−(xi − θ) 2

2

}

∝ exp {

−1 2

n ∑

i=1

(xi − θ)2 }

∝ exp {

−n 2 (x− θ)2

}

e o estimador de máxima verossimilhança é θ̂ = x̄. Usando o algoritmo acima,

gera-se um valor da distribuição Cauchy(0,1) e a probabilidade de aceitação neste

caso fica simplesmente exp[−n(x̄ − θ)2/2]. A função do R a seguir obtém uma amostra de tamanho m de θ e como ilustração vamos gerar 50 observações da

distribuição N(2,1). Note que a taxa de aceitação foi extremamente baixa. Isto

ocorreu devido ao conflito entre verossimilhança e priori.

4.5. MÉTODOS DE REAMOSTRAGEM 59

> rej <- function(x, m) {

+ total = 0

+ theta = rep(0, m)

+ x.bar = mean(x)

+ n = length(x)

+ for (i in 1:m) {

+ accept = FALSE

+ while (!accept) {

+ total = total + 1

+ theta.new = rcauchy(1, 0, 1)

+ prob = exp(-0.5 * n * (theta.new - x.bar)^2)

+ u = runif(1, 0, 1)

+ if (u < prob) {

+ theta = c(theta, theta.new)

+ accept = TRUE

+ }

+ }

+ }

+ cat("\nTaxa de aceitacao", round(m/total, 4), "\n")

+ return(list(theta = theta, total = total))

+ }

> x = rnorm(n = 50, mean = 2, sd = 1)

> m = rej(x, m = 1000)

Taxa de aceitacao 0.0215

O problema é ilustrado na Figura 4.5 (gerada com os comandos abaixo) onde

se pode notar que a maioria dos valores de θ foi gerada em regiões de baixa

verossimilhança.

> x.bar = mean(x)

> x.sd = sd(x)

> curve(dnorm(x, x.bar, x.sd), xlab = expression(theta), from = -4,

+ to = 6, ylab = "", col = 1, lty = 1)

> curve(dcauchy(x, 0, 1), from = -4, to = 6, add = T, lty = 2)

> legend(-3, 0.4, legend = c("priori", "veross."), lty = c(2, 1))

> rug(m$theta)

Mudando a priori para Cauchy(2,1) obtém-se uma taxa de aceitação em torno

de 10% o que ainda constitui uma amostra pequena. Na verdade o número de

simulações deveria ser no mı́nimo 10000 neste caso.

60 CAPÍTULO 4. MÉTODOS APROXIMADOS

−4 −2 0 2 4 6

0. 0

0. 1

0. 2

0. 3

0. 4

θ

priori veross.

Figura 4.5: Verossimilhança normalizada e densidade a priori juntamente com valores simulados.

Portanto, um problema técnico associado ao método é a necessidade de se

maximizar a função de verossimilhança o que pode não ser uma tarefa simples

em modelos mais complexos. Se este for o caso então o método de rejeição

perde o seu principal atrativo que é a simplicidade. Neste caso, o método da

próxima seção passa a ser recomendado. Outro problema é que a taxa de aceitação

pode ser muito baixa. Teremos que gerar muitos valores da distribuição auxiliar

até conseguir um número suficiente de valores da distribuição a posteriori. Isto

ocorrerá se as informações da distribuição a priori e da verossimilhança forem

conflitantes já que neste caso os valores gerados terão baixa probabilidade de

serem aceitos.

4.5.2 Reamostragem Ponderada

Estes métodos usam a mesma idéia de gerar valores de uma distribuição auxiliar

porém sem a necessidade de maximização da verossimilhança. A desvantagem

é que os valores obtidos são apenas aproximadamente distribuidos segundo a

posteriori.

Suponha que temos uma amostra θ1, . . . , θn gerada da distribuição auxiliar q

4.5. MÉTODOS DE REAMOSTRAGEM 61

e a partir dela construimos os pesos

wi = p(θi|x)/q(θi)

∑n j=1 p(θj|x)/q(θj)

, i = 1, . . . , n

O método consiste em tomar uma segunda amostra (ou reamostra) de tamanho

m da distribuição discreta em θ1, . . . , θn com probabilidades w1, . . . , wn. Aqui

também não é necessário que se conheça completamente a posteriori mas apenas

o produto priori vezes verossimilhança já que neste caso os pesos não se alteram.

Tomando novamente a priori como densidade auxiliar, i.e. q(θ) = p(θ) os

pesos se simplificam para

wi = p(x|θi)

∑n j=1 p(x|θj)

, i = 1, . . . , n

e o algoritmo para geração de valores (aproximadamente) da posteriori então fica

1. gerar valores θ1, . . . , θn da distribuição a priori;

2. calcular os pesos wi, i = 1, . . . , n;

3. reamostrar valores com probabilidades w1, . . . , wn.

Este método é essencialmente um bootstrap ponderado. O mesmo problema de

informações conflitantes da priori e da verossimilhança pode ocorrer aqui. Neste

caso, apenas poucos valores gerados da priori terão alta probabilidade de apare-

cerem na reamostra.

Exemplo 4.5 : No Exemplo 4.4, utilizando reamostragem ponderada obtém-se

os gráficos da Figura 4.6.

> reamostra <- function(x, n, m) {

+ x.bar = mean(x)

+ nobs = length(x)

+ theta = rcauchy(n, 0, 1)

+ peso = exp(-0.5 * nobs * (theta - x.bar)^2)

+ aux = sum(peso)

+ peso = peso/aux

+ theta.star = sample(theta, size = m, replace = TRUE, prob = peso)

+ return(list(amostra = theta, pesos = peso, reamostra = theta.star))

+ }

62 CAPÍTULO 4. MÉTODOS APROXIMADOS

−4 −2 0 2 4 6

0. 0

0. 1

0. 2

0. 3

0. 4

θ

2.0 2.1 2.2 2.3 2.4

0. 01

0 0.

02 5

0. 04

0

θ

pe so

s

−4 −2 0 2 4 6

0. 0

0. 1

0. 2

0. 3

0. 4

θ

Figura 4.6: Verossimilhança normalizada (linha cheia), densidade a priori (linha trace- jada) e os valores amostrados (a) e reamostrados (c). Em (b) os valores de θ com pesos maiores do que 0,01.

Exerćıcios

1. Em um modelo de regressão linear simples temos que yi ∼ N(βxi, 1). Os dados observados são y = (−2, 0, 0, 0, 2) e x = (−2,−1, 0, 1, 2), e usamos uma priori vaga N(0, 4) para β. Faça inferência sobre β obtendo uma

amostra da posteriori usando reamostragem ponderada. Compare com a

estimativa de máxima verossimilhança β̂ = 0, 8.

2. Para o mesmo modelo do exerćıcio 1 e os mesmos dados suponha agora

que a variância é desconhecida, i.e. yi ∼ N(βxi, σ2). Usamos uma priori hierárquica para (β, σ2), i.e. β|σ2 ∼ N(0, σ2) e σ−2 ∼ G(0, 01, 0, 01).

(a) Obtenha uma amostra da posteriori de (β, σ2) usando reamostragem

ponderada.

4.6. MONTE CARLO VIA CADEIAS DE MARKOV 63

(b) Baseado nesta amostra, faça um histograma das distribuições

marginais de β e σ2.

(c) Estime β e σ2 usando uma aproximação para a média a posteriori.

Compare com as estimativas de máxima verossimilhança.

4.6 Monte Carlo via cadeias de Markov

Em todos os métodos de simulação vistos até agora obtém-se uma amostra da

distribuição a posteriori em um único passo. Os valores são gerados de forma

independente e não há preocupação com a convergência do algoritmo, bastando

que o tamanho da amostra seja suficientemente grande. Por isso estes métodos

são chamados não iterativos (não confundir iteração com interação). No entanto,

em muitos problemas pode ser bastante dif́ıcil, ou mesmo imposśıvel, encontrar

uma densidade de importância que seja simultaneamente uma boa aproximação

da posteriori e fácil de ser amostrada.

Os métodos de Monte Carlo via cadeias de Markov (MCMC) são uma al-

ternativa aos métodos não iterativos em problemas complexos. A idéia ainda é

obter uma amostra da distribuição a posteriori e calcular estimativas amostrais

de caracteŕısticas desta distribuição. A diferença é que aqui usaremos técnicas de

simulação iterativa, baseadas em cadeias de Markov, e assim os valores gerados

não serão mais independentes.

Nesta seção serão apresentados os métodos MCMC mais utilizados, o

amostrador de Gibbs e o algoritmo de Metropolis-Hastings. A idéia básica é

simular um passeio aleatório no espaço de θ que converge para uma distribuição

estacionária, que é a distribuição de interesse no problema. Uma discussão mais

geral sobre o tema pode ser encontrada por exemplo em Gamerman (1997) e

Gamerman & Lopes (2006).

4.6.1 Cadeias de Markov

Uma cadeia de Markov é um processo estocástico {X0, X1, . . . } tal que a dis- tribuição de Xt dados todos os valores anteriores X0, . . . , Xt−1 depende apenas

de Xt−1. Matematicamente,

P (Xt ∈ A|X0, . . . , Xt−1) = P (Xt ∈ A|Xt−1)

para qualquer subconjunto A. Os métodos MCMC requerem ainda que a cadeia

seja,

ˆ homogênea, i.e. as probabilidades de transição de um estado para outro são

invariantes;

64 CAPÍTULO 4. MÉTODOS APROXIMADOS

ˆ irredut́ıvel, i.e. cada estado pode ser atingido a partir de qualquer outro

em um número finito de iterações;

ˆ aperiódica, i.e. não haja estados absorventes.

e os algoritmos que serão vistos aqui satisfazem a estas condições.

Suponha que uma distribuição π(x), x ∈ Rd seja conhecida a menos de uma constante multiplicativa porém complexa o bastante para não ser posśıvel obter

uma amostra diretamente. Dadas as realizações {X(t), t = 0, 1, . . . } de uma cadeia de Markov que tenha π como distribuição de equilibrio então, sob as

condições acima,

X(t) t→∞−→ π(x) e 1

n

n ∑

t=1

g(X (t) i )

n→∞−→ Eπ(g(Xi)) q.c.

Ou seja, embora a cadeia seja por definição dependente a média aritmética dos

valores da cadeia é um estimador consistente da média teórica.

Uma questão importante de ordem prática é como os valores iniciais influen-

ciam o comportamento da cadeia. A idéia é que conforme o número de iterações

aumenta, a cadeia gradualmente esquece os valores iniciais e eventualmente con-

verge para uma distribuição de equiĺıbrio. Assim, em aplicações práticas é comum

que as iterações iniciais sejam descartadas, como se formassem uma amostra de

aquecimento.

4.6.2 Acurácia Numérica

Na prática teremos um número finito de iterações e tomando

ĝ = 1

n

n ∑

t=1

g(X (t) i )

como estimativa da E(g(Xi)) devemos calcular o seu erro padrão. Como a se-

quência de valores gerados é dependente pode-se mostrar que

V ar(ĝ) = s2

n

[

1 + 2 n

k=1

(

1− k n

)

ρk

]

sendo s2 a variância amostral e ρk a autocorrelação amostral de ordem k. Se

ρk > 0 ∀k então V ar(ĝ) > s2/n. Uma forma muito utilizada para o cálculo da variância do estimador é o método dos lotes aonde os valores da cadeia são

divididos em k lotes de tamanho m e cada lote tem média Bi. O erro padrão de

4.6. MONTE CARLO VIA CADEIAS DE MARKOV 65

ĝ é então estimado como

1

k(k − 1)

k ∑

i=1

(Bi −B)2

sendo m escolhido de modo que a correlação serial de ordem 1 entre as médias

dos lotes seja menor do que 0,05.

Nas próximas seções serão apresentados e discutidos os algoritmos MCMC

mais comumente utilizados.

4.6.3 Algoritmo de Metropolis-Hastings

Os algoritmos de Metropolis-Hastings usam a mesma idéia dos métodos de re-

jeição vistos no caṕıtulo anterior, i.e. um valor é gerado de uma distribuição aux-

iliar e aceito com uma dada probabilidade. Este mecanismo de correção garante

a convergência da cadeia para a distribuição de equilibrio, que neste caso é a

distribuição a posteriori.

Suponha que a cadeia esteja no estado θ e um valor θ′ é gerado de uma

distribuição proposta q(·|θ). Note que a distribuição proposta pode depender do estado atual da cadeia, por exemplo q(·|θ) poderia ser uma distribuição normal centrada em θ. O novo valor θ′ é aceito com probabilidade

α(θ, θ′) = min

(

1, π(θ′) q(θ|θ′) π(θ) q(θ′|θ)

)

. (4.2)

onde π é a distribuição de interesse.

Uma caracteŕıstica importante é que só precisamos conhecer π parcialmente,

i.e. a menos de uma constante já que neste caso a probabilidade (4.2) não se

altera. Isto é fundamental em aplicações Bayesianas aonde não conhecemos com-

pletamente a posteriori. Note também que a cadeia pode permanecer no mesmo

estado por muitas iterações e na prática costuma-se monitorar isto calculando a

porcentagem média de iterações para as quais novos valores são aceitos.

Em termos práticos, o algoritmo de Metropolis-Hastings pode ser especificado

pelos seguintes passos,

1. Inicialize o contador de iterações t = 0 e especifique um valor inicial θ(0).

2. Gere um novo valor θ′ da distribuição q(·|θ).

3. Calcule a probabilidade de aceitação α(θ, θ′) e gere u ∼ U(0, 1).

4. Se u ≤ α então aceite o novo valor e faça θ(t+1) = θ′, caso contrário rejeite e faça θ(t+1) = θ.

66 CAPÍTULO 4. MÉTODOS APROXIMADOS

5. Incremente o contador de t para t+ 1 e volte ao passo 2.

Embora a distribuição proposta possa ser escolhida arbitrariamente na prática

deve-se tomar alguns cuidados para garantir a eficiência do algoritmo. Em apli-

cações Bayesianas a distribuição de interesse é a própria posteriori, i.e. π = p(θ|x) e a probabilidade de aceitação assume uma forma particular,

α(θ, θ′) = min

{

1, p(x|θ′) p(x|θ)

p(θ′)

p(θ)

q(θ|θ′) q(θ′|θ)

}

. (4.3)

O algoritmo será ilustrado nos exemplos a seguir.

Exemplo 4.6 : Em uma certa população de animais sabe-se que cada animal

pode pertencer a uma dentre 4 linhagens genéticas com probabilidades

p1 = 1

2 +

θ

4 , p2 =

1− θ 4

, p3 = 1− θ 4

, p4 = θ

4 .

sendo 0 < θ < 1 um parâmetro desconhecido. Para qualquer θ ∈ (0, 1) é fácil verificar que pi > 0, i = 1, 2, 3, 4 e p1 + p2 + p3 + p4 = 1. Observando-se

n animais dentre os quais yi pertencem à linhagem i então o vetor aleatório

Y = (y1, y2, y3, y4) tem distribuição multinomial com parâmetros n, p1, p2, p3, p4 e portanto,

p(y|θ) = n! y1!y2!y3!y4!

py11 p y2 2 p

y3 3 p

y4 4

∝ (2 + θ)y1(1− θ)y2+y3θy4 .

Atribuindo a distribuição a priori θ ∼ U(0, 1) segue que a densidade a posteriori é proporcional à expressão acima. Então,

p(θ|y) ∝ (2 + θ)y1(1− θ)y2+y3θy4 .

Tomando a distribuição U(0, 1) como proposta então q(θ) = 1, ∀ θ e a probabil- idade (4.3) se simplifica para

α(θ, θ′) = min

{

1, p(x|θ′) p(x|θ)

}

= min

{

1,

(

2 + θ′

2 + θ

)y1 (1− θ′ 1− θ

)y2+y3 (θ′

θ

)y4 }

.

Podemos programar este algoritmo com os comandos do R a seguir.

> p <- function(x, y) {

+ (2 + x)^y[1] * (1 - x)^(y[2] + y[3]) * x^y[4]

+ }

4.6. MONTE CARLO VIA CADEIAS DE MARKOV 67

> metr0 <- function(n, y, fun, start) {

+ theta = c(start, rep(NA, n - 1))

+ taxa = 0

+ for (i in 2:n) {

+ x = runif(1)

+ A = fun(x, y)/fun(theta[i - 1], y)

+ prob = min(1, A)

+ if (runif(1) < prob) {

+ theta[i] = x

+ taxa = taxa + 1

+ }

+ else {

+ theta[i] = theta[i - 1]

+ }

+ }

+ return(list(theta = theta, taxa = taxa/n))

+ }

Suponha que foram observados 197 animais com os números de animais nas

categorias dados por y = (125, 18, 20, 34) e foi gerada uma cadeia de Markov

com 2000 valores de θ. Os valores simulados e as primeiras 30 autocorrelações

amostrais de θ estão na Figura 4.7. A cadeia parece ter convergido após algumas

iterações e podemos descartar os 100 primeiros valores (esta foi a nossa amostra

de aquecimento). Note tambem que a cadeia é altamente correlacionada ao longo

das iterações e isto é devido a alta taxa de rejeição por causa da escolha de q.

> y = c(125, 18, 20, 34)

> n = 2000

> m = metr0(n, y, fun = p, start = 0.5)

> m$taxa

[1] 0.17

Dada uma amostra com valores de θ temos também amostras de valores de

(p1, p2, p3, p4) que podem ser resumidas da seguinte forma,

> p1 = m$theta/4 + 0.5

> p2 = (1 - m$theta)/4

> p3 = p2

> p4 = m$theta/4

> z = as.mcmc(cbind(p1, p2, p3, p4))

> colnames(z) = c("p1", "p2", "p3", "p4")

> b = summary(window(z, start = 501))

> print(b, digits = 3)

68 CAPÍTULO 4. MÉTODOS APROXIMADOS

0 500 1000 1500 2000

0. 50

0. 60

0. 70

(a)

Iterations

0 5 10 15 20 25 30

− 1.

0 0.

0 0.

5 1.

0

(b)

Lag

A ut

oc or

re la

tio n

0.50 0.60 0.70

0 2

4 6

8

(c)

N = 1500 Bandwidth = 0.0106

Figura 4.7: (a) 2000 valores simulados de θ, (b) 30 primeiras autocorrelações amostrais após aquecimento, (c) Densidade a posteriori estimada.

Iterations = 501:2000

Thinning interval = 1

Number of chains = 1

Sample size per chain = 1500

1. Empirical mean and standard deviation for each variable,

plus standard error of the mean:

Mean SD Naive SE Time-series SE

p1 0.6584 0.0114 0.000294 0.000954

p2 0.0916 0.0114 0.000294 0.000954

p3 0.0916 0.0114 0.000294 0.000954

p4 0.1584 0.0114 0.000294 0.000954

4.6. MONTE CARLO VIA CADEIAS DE MARKOV 69

2. Quantiles for each variable:

2.5% 25% 50% 75% 97.5%

p1 0.6340 0.6512 0.6592 0.6656 0.678

p2 0.0721 0.0844 0.0908 0.0988 0.116

p3 0.0721 0.0844 0.0908 0.0988 0.116

p4 0.1340 0.1512 0.1592 0.1656 0.178

Exemplo 4.7 : Suponha que queremos simular valores X ∼ N(0, 1) propondo valores Y ∼ N(x, σ2). Neste caso as densidades propostas no numerador e de- nominador de (4.2) se cancelam e a probabilidade de aceitação fica

α(x, y) = min

{

1, exp

(

−1 2 (y2 − x2)

)}

.

Fixando os valores σ = 0.5 e σ = 10 foram simuladas as cadeias que aparecem na

Figura 4.8. Note que o valor de σ teve um grande impacto na taxa de aceitação

do algoritmo. Isto ocorre porque com σ = 0.5 a distribuição proposta está muito

mais próxima da distribuição de interesse do que com σ = 10.

70 CAPÍTULO 4. MÉTODOS APROXIMADOS

> metrop <- function(n, sigma) {

+ x = c(0, rep(NA, n - 1))

+ for (i in 2:n) {

+ y = rnorm(1, x[i - 1], sigma)

+ prob = min(1, exp(-0.5 * (y^2 - x[i - 1]^2)))

+ u = runif(1, 0, 1)

+ if (u < prob)

+ x[i] = y

+ else x[i] = x[i - 1]

+ }

+ return(x)

+ }

sigma=0.5

Time

0 100 200 300 400 500

− 2

0 2

sigma=10

Time

0 100 200 300 400 500

− 2

0 2

Figura 4.8: 500 valores simulados para o Exemplo 4.7 usando o algoritmo de Metropolis-Hastings com (a) σ = 0.5 e (b) σ = 10.

4.6. MONTE CARLO VIA CADEIAS DE MARKOV 71

Nos Exemplos 4.6 e 4.7 foram ilustrados casos especiais do algoritmo nos quais

a distribuição proposta não depende do estado atual ou a dependência é na forma

de um passeio aleatório. Estes casos são formalizados a seguir.

4.6.4 Casos Especiais

Um caso particular é quando a distribuição proposta não depende do estado atual

da cadeia, i.e. q(θ′|θ) = q(θ′). Em geral, q(·) deve ser uma boa aproximação de π(·), mas é mais seguro se q(·) tiver caudas mais pesadas do que π(·). A probabilidade de aceitação agora fica,

α(θ, θ′) = min

{

1, π(θ′) q(θ)

π(θ) q(θ′)

}

. (4.4)

Note que embora os valores θ′ sejam gerados de forma independente a cadeia

resultante não será i.i.d. já que a probabilidade de aceitação ainda depende de θ.

Outro caso particular é chamado algoritmo de Metropolis e considera apenas

propostas simétricas, i.e., q(θ′|θ) = q(θ|θ′) para todos os valores de θ e θ′. Neste caso a probabilidade de aceitação se reduz para

α(θ, θ′) = min

(

1, π(θ′)

π(θ)

)

.

Um algoritmo de Metropolis muito utilizado é baseado em um passeio aleatório

de modo que a probabilidade da cadeia mover-se de θ para θ′ depende apenas

da distância entre eles, i.e. q(θ′|θ) = q(|θ − θ′|). Neste caso, se usarmos uma distribuição proposta com variância σ2 duas situações extremas podem ocorrer,

1. se σ2 for muito pequena os valores gerados estarão próximos do valor atual

e quase sempre serão aceitos. Mas levará muitas iterações até o algoritmo

cobrir todo o espaço do parâmetro;

2. valores grandes de σ2 levam a uma taxa de rejeição excessivamente alta e a

cadeia se movimenta muito pouco.

Nas duas situações o algoritmo fica ineficiente e na prática temos que tentar vários

valores de σ2.

De um modo geral θ = (θ1, . . . , θd) ′ será um vetor de parâmetros de dimensão

d. Neste caso, pode ser computacionalmente mais eficiente dividir θ em k blo-

cos {θ1, . . . ,θk} e dentro de cada iteração teremos o algoritmo aplicado k vezes. Definindo o vetor θ−i = (θ1, . . . ,θi−1,θi+1, . . . ,θk) que contém todos os elemen-

tos de θ exceto θi suponha que na iteração t+1 os blocos 1, 2, . . . , i− 1 já foram atualizados, i.e.

θ−i = (θ (t+1) 1 , . . . ,θ

(t+1) i−1 ,θ

(t) i+1, . . . ,θ

(t) k ).

72 CAPÍTULO 4. MÉTODOS APROXIMADOS

Para atualizar a i-ésima componente, um valor de θi é gerado da distribuição

proposta q(·|θi,θ−i) e este valor candidato é aceito com probabilidade

α(θi,θ ′ i) = min

{

1, π(θ′i|θ−i) q(θi|θ′i,θ−i)) π(θi|θ−i) q(θ′i|θi,θ−i)

}

. (4.5)

Aqui, π(θi|θ−i) é chamada de distribuição condicional completa como será visto na próxima seção.

Exercicios

1. Assumindo que a distribuição estacionária é N(0, 1),

(a) faça 500 iterações do algoritmo de Metropolis com distribuições pro-

postas N(θ; 0, 5), N(θ; 0, 1) e N(θ, 10).

(b) faça os gráficos dos valores das cadeias ao longo das iterações. Existe

alguma indicação de convergência nos gráficos?

(c) Calcule as taxas de aceitação.

2. Suponha que a distribuição estacionária é N(0, 1).

(a) Para distribuições propostas Cauchy(0, σ), selecione experimental-

mente o valor de σ que maximiza a taxa de aceitação.

(b) Para este valor de σ faça os gráficos dos valores simulados da cadeia

ao longo das iterações e verifique se há indicação de convergência.

(c) Repita os itens anteriores com a distribuição proposta Cauchy(θ, σ).

4.6.5 Amostrador de Gibbs

No amostrador de Gibbs a cadeia irá sempre se mover para um novo valor, i.e não

existe mecanismo de aceitação-rejeição. As transições de um estado para outro

são feitas de acordo com as distribuições condicionais completas π(θi|θ−i), onde θ−i = (θ1, . . . , θi−1, θi+1, . . . , θd)

′.

Em geral, cada uma das componentes θi pode ser uni ou multidimensional.

Portanto, a distribuição condicional completa é a distribuição da i-ésima compo-

nente de θ condicionada em todas as outras componentes. Ela é obtida a partir

da distribuição conjunta como,

π(θi|θ−i) = π(θ)

π(θ)dθi

.

4.6. MONTE CARLO VIA CADEIAS DE MARKOV 73

Assim, para obter a distribuição condicional completa de xi basta pegar os termos

da distribuição conjunta que não dependem de xi.

Exemplo 4.8 : Em um modelo Bayesiano para os dados y que depende dos

parâmetros θ, λ e δ suponha que a distribuição conjunta é dada por

p(y, θ, λ, δ) ∝ p(y|θ, δ)p(θ|λ)p(λ)p(δ).

Após observar y as distribuições a posteriori de cada parâmetro dados todos os

outros são

π(θ|y, λ, δ) ∝ p(y|θ, δ)p(θ|λ) π(λ|y, θ, δ) ∝ p(θ|λ)p(λ) π(δ|y, θ, λ) ∝ p(y|θ, δ)p(δ).

Em muitas situações, a geração de uma amostra diretamente de π(θ) pode

ser custosa, complicada ou simplesmente imposśıvel. Mas se as distribuições

condicionais completas forem completamente conhecidas, então o amostrador de

Gibbs é definido pelo seguinte esquema,

1. inicialize o contador de iterações da cadeia t = 0;

2. especifique valores iniciais θ(0) = (θ (0) 1 , . . . , θ

(0) d )

′;

3. obtenha um novo valor de θ(t) a partir de θ(t−1) através da geração sucessiva

dos valores

θ (t) 1 ∼ π(θ1|θ

(t−1) 2 , θ

(t−1) 3 , . . . , θ

(t−1) d )

θ (t) 2 ∼ π(θ2|θ

(t) 1 , θ

(t−1) 3 , . . . , θ

(t−1) d )

...

θ (t) d ∼ π(θd|θ

(t) 1 , θ

(t) 2 , . . . , θ

(t) d−1)

4. Incremente o contador de t para t + 1 e retorne ao passo 2 até obter con-

vergência.

Assim, cada iteração se completa após d movimentos ao longo dos eixos coorde-

nados das componentes de θ. Após a convergência, os valores resultantes formam

uma amostra de π(θ). Vale notar que, mesmo em problema de grandes dimen-

sões todas as simulações podem ser univariadas, o que em geral é uma vantagem

computacional.

Note também que o amostrador de Gibbs é um caso especial do algoritmo de

Metropolis-Hastings, no qual os elementos de θ são atualizados um de cada vez

74 CAPÍTULO 4. MÉTODOS APROXIMADOS

(ou em blocos), tomando a distribuição condicional completa como proposta e

probabilidade de aceitação igual a 1.

Mais detalhes sobre o amostrado de Gibbs e outros algoritmos relacionados

podem ser obtidos, por exemplo, em Gamerman (1997, Cap. 5) e Robert &

Casella (1999, Cap. 7) .

Exemplo 4.9 : Suponha que Y1, . . . , Yn ∼ N(µ, σ2) com µ e σ2 desconhecidos. Definindo τ = σ−2 a função de verossimilhança é dada por

p(y|µ, τ) ∝ τn/2 exp [

−τ 2

n ∑

i=1

(yi − µ)2 ]

e especificando prioris independentes µ ∼ N(0, s2), sendo s2 a variância amostral e τ ∼ Gama(a, b), com a e b conhecidos, segue que

p(µ, τ |y) ∝ p(y|µ, τ)p(µ)p(τ)

∝ τn/2 exp [

−τ 2

n ∑

i=1

(yi − µ)2 ]

exp

[

− µ 2

2s2

]

τa−1e−bτ .

Esta distribuição conjunta não tem forma padrão mas as condicionais completas

são fáceis de obter,

p(µ|y, τ) ∝ exp [

−τ 2

n ∑

i=1

(yi − µ)2 ]

exp

[

− µ 2

2s2

]

∝ exp [

−1 2 (nτ + s−2)µ2 − 2µȳ)

]

∝ exp [

− 1 2C

(µ−m)2 ]

onde C−1 = nτ + s−2 e m = Cȳ e

p(τ |y, µ) ∝ τa+n/2−1 exp [

−τ (

b+ 1

2

n ∑

i=1

(yi − µ)2 )]

.

Segue então que

µ|y, τ ∼ N(m,C)

τ |y, µ ∼ Gama (

a+ n

2 , b+

1

2

n ∑

i=1

(yi − µ)2 )

e o amostrador de Gibbs pode ser implementado facilmente gerando valores destas

distribuições alternadamente.

Exemplo 4.10 : Em um processo de contagem no qual foram observados

4.6. MONTE CARLO VIA CADEIAS DE MARKOV 75

Y1, . . . , Yn suspeita-se que houve um ponto de mudança m tal que

Yi ∼ Poisson(λ), i = 1, . . . ,m Yi ∼ Poisson(φ), i = m+ 1, . . . , n.

O objetivo é estimar o ponto de mudança m e os parâmetros dos 2 processos de

Poisson. Assumindo-se as distribuições a priori independentes

λ ∼ Gama(a, b) φ ∼ Gama(c, d) m ∼ Uniforme{1, . . . , n}

a densidade a posteriori fica

p(λ, φ,m|y) ∝ m ∏

i=1

e−λλyi n ∏

i=m+1

e−φφyiλa−1e−bλφc−1e−dφ 1

n

∝ λa+t1−1e−(b+m)λφc+t2−1e−(d+n−m)φ 1 n

sendo t1 = ∑m

i=1 yi e t2 = ∑n

i=m+1 yi. Neste caso não é dificil verificar que as

distribuições condicionais completas ficam

p(λ|φ,m,y) ∝ λa+t1−1e−(b+m)λ ou λ|φ,m,y ∼ Gama(a+ t1, b+m) p(φ|λ,m,y) ∝ φc+t2−1e−(d+n−m)φ ou φ|λ,m,y ∼ Gama(c+ t2, d+ n−m) p(m|λ, φ,y) ∝ λt1e−mλφt2e−(n−m)φ, m = 1, . . . , n.

A função do R abaixo obtem uma amostra da posteriori conjunta simulando val-

ores destas condicionais completas.

76 CAPÍTULO 4. MÉTODOS APROXIMADOS

> Gibbs <- function(a, b, c, d, y, niter) {

+ N = length(y)

+ lambda = phi = m = matrix(0, nrow = niter)

+ lambda[1] = 1

+ phi[1] = 1

+ m[1] = 10

+ for (i in 2:niter) {

+ t1 = sum(y[1:m[i - 1]])

+ t2 = 0

+ if (m[i - 1] < N)

+ t2 = sum(y[(m[i - 1] + 1):N])

+ lambda[i] = rgamma(1, (a + t1), (b + m[i - 1]))

+ phi[i] = rgamma(1, (c + t2), (d + N - m[i - 1]))

+ prob = NULL

+ for (j in 1:N) {

+ t1 = sum(y[1:j])

+ t2 = 0

+ if (j < N) {

+ t2 = sum(y[(j + 1):N])

+ }

+ aux = (lambda[i]^t1) * exp(-j * lambda[i]) * (phi[i]^t2) *

+ exp(-(N - j) * phi[i])

+ prob = c(prob, aux)

+ }

+ soma = sum(prob)

+ probm = prob/soma

+ m[i] = sample(x = N, size = 1, prob = probm)

+ }

+ return(list(lambda = lambda, phi = phi, m = m))

+ }

Testando a função Gibbs com 40 dados simulados de processos com médias 2

e 5 e ponto de mudança 23.

> y = c(rpois(n = 22, lambda = 2), rpois(n = 18, lambda = 5))

> x = Gibbs(a = 0.1, b = 0.1, c = 0.1, d = 0.1, y = y, niter = 2000)

Podemos usar o pacote coda para analisar os valores simulados. As 1000

primeiras simulações são descartadas como amostra de aquecimento.

> library(coda)

> amostra = cbind(x$lambda, x$phi, x$m)[1001:2000, ]

4.6. MONTE CARLO VIA CADEIAS DE MARKOV 77

> theta = mcmc(amostra)

> colnames(theta) = names(x)

> summary(theta)

Iterations = 1:1000

Thinning interval = 1

Number of chains = 1

Sample size per chain = 1000

1. Empirical mean and standard deviation for each variable,

plus standard error of the mean:

Mean SD Naive SE Time-series SE

lambda 2.273 0.3247 0.01027 0.00865

phi 5.246 0.5569 0.01761 0.02049

m 21.612 1.6125 0.05099 0.06403

2. Quantiles for each variable:

2.5% 25% 50% 75% 97.5%

lambda 1.668 2.054 2.258 2.479 2.979

phi 4.213 4.843 5.230 5.610 6.398

m 18.975 21.000 22.000 22.000 24.025

A partir dos valores simulados de m podemos estimar suas probabilidades,

> tm = table(theta[, 3])/1000

> print(tm)

7 8 9 10 11 14 15 16 17 18 19 20 21

0.001 0.002 0.001 0.001 0.001 0.005 0.002 0.004 0.001 0.007 0.012 0.059 0.196

22 23 24 25 26 27

0.648 0.010 0.025 0.010 0.013 0.002

Finalmente, pode-se estimar as contagens médias condicionando nos valor de

m com maior probabilidade.

> lambda.22 = theta[, 1][theta[, 3] == 22]

> phi.22 = theta[, 2][theta[, 3] == 22]

> theta.22 = as.mcmc(cbind(lambda.22, phi.22))

78 CAPÍTULO 4. MÉTODOS APROXIMADOS

> plot(theta)

0 200 400 600 800 1000

1. 5

2. 5

Iterations

Trace of lambda

1.5 2.0 2.5 3.0 3.5

0. 0

0. 6

1. 2

N = 1000 Bandwidth = 0.08448

Density of lambda

0 200 400 600 800 1000

4. 0

5. 5

7. 0

Iterations

Trace of phi

4 5 6 7

0. 0

0. 3

0. 6

N = 1000 Bandwidth = 0.1483

Density of phi

0 200 400 600 800 1000

10 20

Iterations

Trace of m

y

D en

si ty

10 15 20 25

0. 0

0. 2

0. 4

Density of m

Figura 4.9: rtwert

4.7 Problemas de Dimensão Variável

Em muitas aplicações práticas é razoável assumir que existe incerteza também em

relação ao modelo que melhor se ajusta a um conjunto de dados. Do ponto de vista

Bayesiano esta incerteza é simplesmente incorporada ao problema de inferência

considerando-se o próprio modelo como mais um parâmetro desconhecido a ser

estimado. Assim os diferentes modelos terão uma distribuição de probabilidades.

Para isto vamos criar uma variável aleatória discreta k que funciona como

indicador de modelo e atribuir probabilidades a priori p(k) para cada modelo.

Além disso, para cada k existe um vetor de parâmetros θ(k) ∈ Rnk com

ˆ uma verossimilhança p(y|θ(k), k)

ˆ uma distribuição a priori p(θ(k)|k).

4.7. PROBLEMAS DE DIMENSÃO VARIÁVEL 79

> plot(tm)

0. 0

0. 1

0. 2

0. 3

0. 4

0. 5

0. 6

tm

7 9 11 14 16 18 20 22 24 26

Figura 4.10:

Se M é conjunto de todos os posśıveis modelos (ou modelos candidatos), então

as probabilidades a posteriori de cada posśıvel modelo são dadas por

π(k|y) = p(k) p(y|k)∑

k∈M

p(k) p(y|k) , k ∈ M

sendo p(y|k) a verossimilhança marginal obtida como

p(y|k) = ∫

p(y|θ, k)p(θ|k)dθ.

O problema aqui é que esta última integral só é analiticamente tratável em alguns

casos restritos. Além disso, se o número de modelos candidatos for muito grande

calcular (ou aproximar) p(y|k) pode ser inviável na prática.

80 CAPÍTULO 4. MÉTODOS APROXIMADOS

> plot(theta.22)

0 100 300 500

1. 5

2. 0

2. 5

3. 0

Iterations

Trace of lambda.22

1.5 2.0 2.5 3.0 3.5

0. 0

0. 4

0. 8

1. 2

N = 648 Bandwidth = 0.08688

Density of lambda.22

0 100 300 500

4. 0

5. 0

6. 0

7. 0

Iterations

Trace of phi.22

4 5 6 7

0. 0

0. 2

0. 4

0. 6

N = 648 Bandwidth = 0.1611

Density of phi.22

Figura 4.11:

Por outro lado, se for especificada a distribuição de interesse como a seguinte

posteriori conjunta,

π(θ, k|y) ∝ p(y|θ, k) p(θ|k) p(k)

e conseguirmos simular valores desta distribuição então automaticamente teremos

uma amostra aproximada de π(k|y) e π(θ|k,y). Note que neste caso estamos admitindo que a dimensão de θ pode variar ao

longo dos modelos e precisamos então construir uma cadeia com espaço de esta-

dos que muda de dimensão ao longo das iterações. Os algoritmos de Metropolis-

Hastings e o amostrador de Gibbs não podem ser utilizados já que são definidos

apenas para distribuições com dimensão fixa. Embora existam outras possibili-

dades iremos estudar os algoritmos MCMC com saltos reversiveis (Green 1995)

que são particularmente úteis no contexto de seleção Bayesiana de modelos.

4.7. PROBLEMAS DE DIMENSÃO VARIÁVEL 81

4.7.1 MCMC com Saltos Reversiveis (RJMCMC)

Este algoritmo é baseado na abordagem usual dos métodos de Metropolis-

Hastings de propor um novo valor para a cadeia e definir uma probabilidade

de aceitação. No entanto, os movimentos podem ser entre espaços de dimen-

sões diferentes como veremos a seguir. Em cada iteração o algoritmo envolve a

atualização dos parâmetros, dado o modelo, usando os métodos MCMC usuais

discutidos anteriormente e a atualização da dimensão usando o seguinte proced-

imento.

Suponha que o estado atual da cadeia é (k,θ), i.e. estamos no modelo k

com parâmetros θ e um novo modelo k′ com parâmetros θ′ é proposto com

probabilidade rk,k′ . Em geral isto significa incluir ou retirar parâmetros do modelo

atual. Vamos assumir inicialmente que o modelo proposto tem dimensão maior,

i.e. nk′ > nk e que θ ′ = g(θ,u) para uma função deterministica g e um vetor

aleatório u ∼ q(u) com dimensão nk′−nk. Então o seguinte algoritmo é utilizado,

ˆ proponha (k,θ) → (k′,θ′) com probabilidade rk,k′

ˆ gere u ∼ q(u) com dimensão nk′ − nk

ˆ faça θ′ = g(θ,u),

ˆ aceite (k′,θ′) com probabilidade min(1, A) sendo

A = π(k′,θ′)

π(k,θ) × rk′,k

rk,k′ q(u)

∂g(θ,u)

∂(θ,u)

.

Exemplo 4.11 : Sejam Y1, . . . , Yn os tempos de vida de componentes eletrônicos

sorteados ao acaso e existe incerteza em relação a distribuição dos dados. Sabe-se

que

Yi ∼ Exp(λ) (Modelo 1) ou Yi ∼ Gama(α, β) (Modelo 2), i = 1, . . . , n.

Suponha que atribuimos as probabilidades a priori p(k) = 1/2 para o indicador

de modelo e as seguintes distribuições a priori foram atribuidas aos parâmetros

dentro de cada modelo,

λ|k = 1 ∼ Gama(2, 1) α|k = 2 ∼ Gama(4, 2) e β|k = 2 ∼ Gama(4, 2).

Dado o modelo, as funções de verossimilhança ficam

p(y|λ, k = 1) = λne−λ ∑

yi

82 CAPÍTULO 4. MÉTODOS APROXIMADOS

p(y|α, β, k = 2) = β nα

Γn(α)

yα−1i e −β

∑ yi

as distribuições condicionais completas são facilmente obtidas como

λ|y, α, β, k = 1 ∼ Gama(n+ 2, 1 + ∑

yi)

β|y, α, λ, k = 2 ∼ Gama(nα + 4, 2 + ∑

yi)

p(α|y, β, λ, k = 2) ∝ β nα

Γn(α)

yα−1i α 3e−2α

A distribuição condicional completa de α não é conhecida então vamos usar o

algoritmo de Metropolis-Hastings propondo valores α′ ∼ U [α − ǫ, α + ǫ]. A função a seguir atualiza o valor de α segundo este esquema.

> mh.alpha <- function(y, n, alpha, beta, eps) {

+ z = runif(1, alpha - eps, alpha + eps)

+ if (z <= 0) {

+ acc = 0

+ }

+ else {

+ t1 = prod(y)

+ num = beta^(n * z) * t1^(z - 1)/(gamma(z)^n)

+ den = beta^(n * alpha) * t1^(alpha - 1)/(gamma(alpha)^n)

+ num = num * exp(-2 * z) * z^3

+ den = den * exp(-2 * alpha) * alpha^3

+ }

+ aceita = min(1, num/den)

+ u = runif(1)

+ newalpha = ifelse(u < aceita, z, alpha)

+ return(newalpha)

+ }

Suponha que o modelo atual é Exp(λ) e queremos propor o modelo

Gama(α, β). Um possivel esquema de atualização é o seguite,

1. gere u ∼ Gama(a, b)

2. defina (α, β) = g(λ, u) = (u, λu)

3. calcule o Jacobiano, ∣

0 1

u λ

= u

4.7. PROBLEMAS DE DIMENSÃO VARIÁVEL 83

4. aceite o novo modelo com probabilidade min(1, A) sendo

A = p(y | α, β, k = 2) p(y | λ, k = 1)

p(α)p(β)

p(λ)

u

q(u)

Note que transformação no item (2) preserva a média, ou seja E(Y ) = 1/λ sob o

modelo exponencial e E(Y ) = u/λu = 1/λ sob o modelo gama.

Se o modelo atual for Gama(α, β) e propomos o modelo Exp(λ) o esquema

reverso consiste em fazer (λ, u) = g−1(α, β) = (β/α, α). A probabilidade de

aceitação é simplesmente min(1, 1/A) substituindo u = α.

> rj.modelo <- function(y, n, lambda, alpha, beta, model, a, b) {

+ if (model == 1) {

+ u = rgamma(1, a, b)

+ alpha1 = u

+ beta1 = lambda * u

+ lambda1 = lambda

+ }

+ else {

+ lambda1 = beta/alpha

+ alpha1 = alpha

+ beta1 = beta

+ u = alpha

+ }

+ t1 = prod(y)

+ t2 = sum(y)

+ num = beta1^(n * alpha1) * t1^(alpha1 - 1) * exp(-beta1 *

+ t2)/(gamma(alpha1)^n)

+ num = num * 2^4 * alpha1^3 * exp(-2 * alpha1)/gamma(4)

+ num = num * 2^4 * beta1^3 * exp(-2 * beta1)/gamma(4) * alpha1

+ den = (lambda1^n) * exp(-lambda1 * t2)

+ den = den * lambda1 * exp(-lambda1)/gamma(2)

+ den = den * b^a * u^(a - 1) * exp(-b * u)/gamma(a)

+ u = runif(1, 0, 1)

+ if (model == 1) {

+ aceita = min(1, num/den)

+ if (u < aceita) {

+ model = 2

+ alpha = alpha1

+ beta = beta1

+ }

+ }

84 CAPÍTULO 4. MÉTODOS APROXIMADOS

+ else {

+ aceita = min(1, den/num)

+ if (u < aceita) {

+ model = 1

+ lambda = lambda1

+ }

+ }

+ if (model == 1)

+ return(list(model = model, lambda = lambda))

+ else return(list(model = model, alpha = alpha, beta = beta))

+ }

Finalmente o algoritmo pode ser implementado para atualizar tanto o modelo

quanto os parâmetros dentro do modelo.

> rjmcmc <- function(niter, nburn, y, n, a, b, eps = 0.25) {

+ x = matrix(0, nrow = niter + 1, ncol = 3)

+ x1 = matrix(0, nrow = niter - nburn, ncol = 3)

+ nv = array(0, 2)

+ nv1 = array(0, 2)

+ x[1, (1:3)] = c(1, 1, 1)

+ model = 1

+ t1 = prod(y)

+ t2 = sum(y)

+ for (i in 1:niter) {

+ if (model == 1) {

+ x[nv[1] + 1, 1] = rgamma(1, n + 2, t2 + 1)

+ }

+ else {

+ x[nv[2] + 1, 3] = rgamma(1, 4 + n * x[nv[2], 2],

+ t2 + 2)

+ x[nv[2] + 1, 2] = mh.alpha(y, n, x[nv[2], 2], x[nv[2] +

+ 1, 3], eps)

+ }

+ new = rj.modelo(y, n, x[nv[1] + 1, 1], x[nv[2] + 1, 2],

+ x[nv[2] + 1, 3], model, a, b)

+ model = new$model

+ if (model == 1) {

+ x[nv[1] + 1, 1] = new$lambda

+ nv[1] = nv[1] + 1

+ if (i > nburn) {

+ x1[nv1[1] + 1, 1] = new$lambda

4.7. PROBLEMAS DE DIMENSÃO VARIÁVEL 85

+ nv1[1] = nv1[1] + 1

+ }

+ }

+ else {

+ x[nv[2] + 1, 2] = new$alpha

+ x[nv[2] + 1, 3] = new$beta

+ nv[2] = nv[2] + 1

+ if (i > nburn) {

+ x1[nv1[2] + 1, 2] = new$alpha

+ x1[nv1[2] + 1, 3] = new$beta

+ nv1[2] = nv1[2] + 1

+ }

+ }

+ }

+ cat("Probabilidades a posteriori dos modelos", "\n")

+ print(nv1/(niter - nburn))

+ cat("Medias a posteriori dos parametros", "\n")

+ somas = apply(x1, 2, sum)

+ print(somas/c(nv1[1], nv1[2], nv1[2]))

+ return(list(x = x, nv = nv, x1 = x1, nv1 = nv1))

+ }

Vamos testar as funções acima simulando um conjunto de dados com dis-

tribuição exponencial.

> y = rexp(10, 3)

> niter = 1000

> nburn = 500

> m = rjmcmc(1000, 500, y, 10, 1, 1)

Probabilidades a posteriori dos modelos

[1] 0.8 0.2

Medias a posteriori dos parametros

[1] 3.794036 1.044988 3.439110

Assim o modelo exponencial tem probabilidade a posteriori bem maior que o

modelo gama. Podemos estar interessados em estimar os tempos médios de vida

(E(Y )) sob cada modelo.

> r1 = 1:m$nv1[1]

> r2 = 1:m$nv1[2]

> x = m$x1[, c(1, 2)]

> x[r1, 1] = 1/m$x1[r1, 1]

86 CAPÍTULO 4. MÉTODOS APROXIMADOS

> x[r2, 2] = m$x1[r2, 2]/m$x1[r2, 3]

> somas = apply(x, 2, sum)

> medias = somas/c(m$nv1[1], m$nv1[2])

> print(medias)

[1] 0.2892936 0.3186531

> prob = m$nv1/(niter - nburn)

> prob[1] * medias[1] + prob[2] * medias[2]

[1] 0.2951655

4.8 Tópicos Relacionados

4.8.1 Autocorrelação Amostral

Em uma cadeia de Markov, os valores gerados são por definição correlacionados

ao longo das iterações pois o valor de θ(t) foi gerado a partir de θ(t−1). Em

muitas situações estes valores podem ser altamente correlacionados e em geral a

autocorrelação será positiva. Ou seja, pode não haver muito ganho em termos

de informação em se armazenar todos os valores simulados da cadeia e podemos

estar desperdiçando espaço em disco, especialmente se a dimensão do problema

for muito grande.

Embora não tenha nenhuma justificativa teórica, uma abordagem prática

muito utilizada consiste em guardar os valores simulados a cada k iterações. Neste

caso, dizemos que as simulações foram feitas com thinning igual a k. Por exemplo,

se foram feitas 100 mil simulações, descartadas as 50 mil primeiras e guardados

os valores a cada 10 iterações então no final as inferências serão baseadas em uma

amostra de tamanho 5000.

Comentário

A não ser para obter esta redução de espaço ocupado em disco, descartar valores

simulados (além daqueles da amostra de aquecimento) me parece um desperdicio.

Métodos de séries temporais estão disponiveis para analisar cadeias levando em

conta as autocorrelações. Além disso pode-se tentar outros amostradores que

gerem cadeias com menor autocorrelação amostral.

4.8.2 Monitorando a Convergência

Aqui vale lembrar que a verificação de convergência (ou falta de convergência) é

responsabilidade do analista. Além disso estamos falando de convergência para

4.8. TÓPICOS RELACIONADOS 87

a distribuição alvo, que neste caso é a distribuição a posteriori, o que pode ser

extremamente dif́ıcil de se verificar na prática.

Caṕıtulo 5

Modelos Lineares

Em uma situação mais geral, a variável de interesse (variável resposta) tem sua

descrição probabiĺıstica afetada por outras variáveis (variáveis explicativas ou

covariáveis). No caso mais simples a influência sobre a resposta média é linear e

aditiva e pode ser vista como uma aproximação de primeira ordem para funções

mais complexas.

Usando uma notação matricial, o modelo linear normal pode ser escrito como

y = Xβ + ǫ,

onde y é um vetor n × 1 de observações, X é uma matriz n × p conhecida, β é um vetor p × 1 de parâmetros e ǫ é um vetor n × 1 de erros aleatórios tais que ǫi ∼ N(0, σ2) e E(ǫiǫj) = 0, para i = 1, · · · , n e j 6= i. O modelo nos diz então que, a distribuição condicional de y dados β e σ2 é normal multivariada,

i.e. y ∼ N(Xβ, σ2In) sendo In é a matriz identidade de ordem n. Definindo φ = σ−2 e usando a função de densidade da normal multivariada (ver apêndice

A) segue que

f(y|β, φ) = (2π)−n/2|φ−1In|−1/2 exp {

−1 2 (y −Xβ)′(φ−1In)−1(y −Xβ)

}

∝ φn/2 exp {

−φ 2 (y −Xβ)′(y −Xβ)

}

. (5.1)

A forma quadrática em (5.1) pode ser reescrita em termos de β̂ = (X ′X)−1X ′y

que é o estimador de máxima verossimilhança de β,

(y −Xβ)′(y −Xβ) = (y −Xβ̂ −X(β − β̂))′(y −Xβ̂ −X(β − β̂)) = (y −Xβ̂)′(y −Xβ̂) + (β −Xβ̂)′X ′X(β −Xβ̂)

−2(β −Xβ̂)X ′(y −Xβ̂) = (y −Xβ̂)′(y −Xβ̂) + (β −Xβ̂)′X ′X(β −Xβ̂)

88

89

pois X ′(y − Xβ̂) = 0. Denotando por S = (y − Xβ̂)′(y − Xβ̂) a soma de quadrados residual, podemos escrever então a função de verossimilhança como,

f(y|β, φ) ∝ φn/2 exp {

−φ 2 [(β − β̂)′X ′X(β − β̂) + S]

}

.

A distribuição a priori adotada aqui é uma generalização multivariada da

distribuição Normal-Gama vista na Seção 2.3.5. Assim, a distribuição a priori é

especificada como

β|φ ∼ Np(µ0, (C0φ)−1)

onde C0 é agora uma matriz p× p e

φ ∼ Gama (

n0 2 , n0σ

2 0

2

)

.

Com isso a densidade a priori conjunta de (β, φ) fica completamente especificada

e assim como no caso univariado a distribuição marginal de β é obtida integrando-

se p(β, φ) em relação a φ onde,

p(β, φ) ∝ φ n0+p

2 −1 exp

{

−φ 2

[

n0σ 2 0 + (β

′ − µ0)′C0(β′ − µ0) ]

}

.

É fácil verificar que

p(β) ∝ [

1 + (β − µ0)′C0(β − µ0)

n0σ20

]−(n0+p)/2

de modo que a distribuição a priori marginal de β é β ∼ tn0(µ0, σ20C−10 ). Note que, como C0 é simétrica, é necessário especificar p(p + 1)/2 de seus elementos.

Na prática, podemos simplificar esta especificação assumindo que C0 é diagonal,

i.e. que os componentes de β são não correlacionados a priori.

Combinando-se com a verossimilhança via teorema de Bayes obtem-se as

seguintes distribuições a posteriori

β|φ,y ∼ N(µ1, (C1φ)−1)

φ|y ∼ Gama (

n1 2 , n1σ

2 1

2

)

ou n1σ 2 1φ ∼ χ2n1

β|y ∼ tn1(µ1, σ21C−11 )

90 CAPÍTULO 5. MODELOS LINEARES

onde os parâmetros atualizados são

n1 = n0 + n

C1 = C0 +X ′X

µ1 = (C0 +X ′X)−1(C0µ0 +X

′Xβ̂)

n1σ 2 1 = n0σ

2 0 + (y −Xµ1)′y + (µ0 − µ1)′C0µ0

= n0σ 2 0 + (n− p)σ̂2 + (µ0 − β̂)′[C−10 +X ′X−1]−1(µ0 − β̂)

onde

σ̂2 = 1

n− p(y −Xβ̂) ′(y −Xβ̂).

Os estimadores pontuais de β e φ são dados respectivamente por µ1 e σ −2 1 .

Intervalos de confiança para βj e φ são obtidos através dos percentis das

distribuições univariadas tn1(µj, σ 2 1(C

−1 1 )jj), j = 1, · · · , p e χ2n1 . Em particular,

note que µ1 é obtida como uma ponderação matricial entre a estimativa a priori

de β e sua estimativa de máxima verossimilhança β̂. Inferência conjunta sobre

β também pode ser feita usando o fato que a forma quadrática

(β − µ1)′C1(β − µ1)/p σ21

∼ F (p, n1).

Note que o modelo visto na seção anterior é na verdade o caso mais simples

de um modelo linear quando p = 1 e X é um vetor n× 1 de 1’s. Neste caso β é um escalar podendo ser denotado por µ e o modelo se reduz a yi = µ+ ǫi.

A priori não informativa é também uma generalização multivariada da seção

anterior. Aqui o vetor β é um parâmetro de locação e φ é um parâmetro de escala,

e portanto a priori não informativa de Jeffreys é p(β, φ) ∝ φ−1. Vale notar que esta priori é um caso particular (degenerado) da priori conjugada natural

com C0 = 0 e n0 = −p. Fazendo as substituições adequadas obtém-se que as distribuições a posteriori são dadas por

β|y ∼ tn−p(β̂, s2(X ′X)−1) (n− p)s2φ|y ∼ χ2n−p

(β − β̂)′X ′X(β − β̂) s2

|y ∼ F (p, n− p)

e estimadores pontuais bem como intervalos de confiança coincidirão com os obti-

dos usando métodos clássicos.

5.1. ANÁLISE DE VARIÂNCIA COM 1 FATOR DE CLASSIFICAÇÃO 91

5.1 Análise de Variância com 1 Fator de Classi-

ficação

Considere o modelo yij = βj + ǫij , i = 1, · · · , nj e j = 1, · · · , p. Assim, todas as nj observações do grupo j têm a mesma média βj . Neste problema, o número

total de observações independentes é n = n1 + · · · + np. Em outras palavras, Y1j, · · · , Ynjj ∼ N(βj, σ2). Se os yij forem “empilhados” em um único vetor n× 1 então podemos reescrever o modelo na forma matricial y = Xβ + ǫ sendo

X =

1 0 · · · 0 ...

... ...

1 0 · · · 0 ...

... ...

0 0 · · · 1 ...

... ...

0 0 · · · 1

.

Note que X ′X = diagonal(n1, · · · , np) e a forma quadrática (β−β̂)′X ′X(β−β̂) se reduz a

p ∑

j=1

nj(βj − yj)2

e a função de verossimilhança é dada por

l(β1, · · · , βp, φ; y) ∝ φn/2 exp {

−φ 2

[

(n− p)s2 + p

j=1

nj(βj − yj)2 ]}

com

s2 = 1

n− p(y −Xβ̂) ′(y −Xβ̂).

Assumindo que βj|φ ∼ N(µj, (cjφ)−1), j = 1, · · · , p são condicionalmente independentes e que n0σ

2 0φ ∼ χ2n0 então as distribuições a posteriori são

βj|φ, y ∼ N(µ∗j , (c∗jφ)−1) n1σ

2 1φ|y ∼ χ2n1 βj|y ∼ tn1(µ∗j , σ21/c∗j)

92 CAPÍTULO 5. MODELOS LINEARES

onde

µ∗j = cjµj + njyj cj + nj

c∗j = cj + nj

n1 = n0 + n

n1σ 2 1 = n0σ

2 0 + (n− p)s2 +

p ∑

i=1

njcj cj + nj

(yj − µj)2

e os βj|φ, y permanecem independentes. A priori não informativa p(β, φ) ∝ φ−1 é obtida fazendo-se cj = 0, j = 1, · · · , p

e n0 = −p. Assim, as distribuições a posteriori marginais são dadas por

βj|y ∼ tn−p(yj, s2/nj) (n− p)s2φ ∼ χ2n−p

e as estimativas pontuais e intervalos de confiança coincidirão com os da inferência

clássica. Em particular, se estamos interessados em testar

H0 : β1 = · · · = βp = β

então pode-se mostrar que (DeGroot,1970, páginas 257 a 259) deve-se rejeitar H0 se

P

F >

p ∑

j=1

nj(yj − y)2/(p− 1)

s2

onde F ∼ F (p− 1, n− p) for pequena. Note que as hipóteses equivalentes são

H0 : α1 = · · · = αp = 0

sendo

αj = βj − β, β = 1

n

p ∑

j=1

njβj e

p ∑

j=1

njαj = 0

e αj é o efeito da j-ésima população. Neste caso, X ′X = diagonal(n1, · · · , np) e

a forma quadrática (β− β̂)′X ′X(β− β̂) fica ∑

nj(αj − yj − y)2 + n(β− yj − y)2.

Apêndice A

Lista de Distribuições

Neste apêndice são listadas as distribuições de probabilidade utilizadas no texto

para facilidade de referência. Sã apresentadas suas funções de (densidade) de

probabilidade além da média e variância. Uma revisão exaustiva de distribuições

de probabilidades pode ser encontrada em Johnson et al. (1992, 1995) e Evans

et al. (1993).

A.1 Distribuição Normal

X tem distribuição normal com parâmetros µ ∈ R e σ2 > 0, denotando-se X ∼ N(µ, σ2), se sua função de densidade é dada por

p(x|µ, σ2) = (2πσ2)−1/2 exp [

−1 2

(x− µ)2 σ2

]

, −∞ < x < ∞.

E(X) = µ e V (X) = σ2.

Quando µ = 0 e σ2 = 1 a distribuição é chamada normal padrão.

No caso vetorial, X = (X1, . . . , Xp) tem distribuição normal multivari-

ada com vetor de médias µ e matriz de variância-covariância Σ, denotando-se

X ∼ N(µ,Σ) se sua função de densidade é dada por

p(x|µ,Σ) = (2π)−p/2|Σ|−1/2 exp[−(x− µ)′Σ−1(x− µ)/2]

para µ ∈ Rp e Σ positiva-definida.

93

94 APÊNDICE A. LISTA DE DISTRIBUIÇÕES

A.2 Distribuição Log-Normal

Se X ∼ N(µ, σ2) então Y = eX tem distribuição log-normal com parâmetros µ e σ2. Portanto, sua função de densidade é dada por

p(y|µ, σ2) = 1 y (2πσ2)−1/2 exp

[

−1 2

(log(y)− µ)2 σ2

]

, −∞ < x < ∞.

E(X) = exp{µ+ σ2/2} e V (X) = exp{2µ+ σ2}(exp{σ2} − 1).

A.3 A Função Gama

Γ(α) =

∫ ∞

0

xα−1e−xdx.

Propriedades,

ˆ Usando integração por partes pode-se mostrar que,

Γ(α + 1) = αΓ(α), α > 0.

ˆ Γ(1) = 1.

ˆ Γ(1/2) = √ π.

ˆ Para n um inteiro positivo,

Γ(n+ 1) = n! e Γ

(

n+ 1

2

)

=

(

n− 1 2

)(

n− 3 2

)

. . . 3

2

1

2

√ π

A.4 Distribuição Gama

X tem distribuição Gama com parâmetros α > 0 e β > 0, denotando-se X ∼ Ga(α, β), se sua função de densidade é dada por

p(x|α, β) = β α

Γ(α) xα−1e−βx, x > 0.

E(X) = α/β e V (X) = α/β2.

Casos particulares da distribuição Gama são a distribuição de Erlang, Ga(α, 1),

a distribuição exponencial, Ga(1, β), e a distribuição qui-quadrado com ν graus

de liberdade, Ga(ν/2, 1/2).

A.5. DISTRIBUIÇÃO WISHART 95

A.5 Distribuição Wishart

Diz-se que uma matriz aleatória Ω (n× n) segue uma distribuição Wishart com parâmetro Σ e ν graus de liberdade, denotando-se Ω ∼ W (Σ, ν), se sua função de densidade é dada por,

p(Ω|Σ, ν) ∝ |Ω|(ν−n−1)/2 exp(−(1/2)tr(ΣΩ))

sendo ν ≥ n, Σ positiva-definida e tr(A) indica o traço de uma matriz A. Uma propriedade útil é que AΩA′ ∼ W (AΣA′, ν).

A.6 Distribuição Gama Inversa

X tem distribuição Gama Inversa com parâmetros α > 0 e β > 0, denotando-se

X ∼ GI(α, β), se sua função de densidade é dada por

p(x|α, β) = β α

Γ(α) x−(α+1) e−β/x, x > 0.

E(X) = β

α− 1 , para α > 1 e V (X) = β2

(α− 1)2(α− 2) , para α > 2.

Não é dif́ıcil verificar que esta é a distribuição de 1/X quando X ∼ Ga(α, β).

A.7 Distribuição Wishart Invertida

Diz-se que uma matriz aleatória Ω (n × n) segue uma distribuição Wishart- Invertida com parâmetro Σ e ν graus de liberdade, denotando-se Ω ∼ WI(Σ, ν) se sua função de densidade é dada por,

p(Ω|Σ, ν) ∝ |Ω|−(ν+n+1)/2 exp(−(1/2)tr(ΣΩ))

sendo ν ≥ n, Σ positiva-definida e tr(A) indica o traço de uma matriz A. Não é dif́ıcil verificar que Ω−1 ∼ W (Σ, ν). Outra propriedade é que AΩA′ ∼ WI(AΣA′, ν).

96 APÊNDICE A. LISTA DE DISTRIBUIÇÕES

A.8 Distribuição Beta

X tem distribuição Beta com parâmetros α > 0 e β > 0, denotando-se

X ∼ Be(α, β), se sua função de densidade é dada por

p(x|α, β) = Γ(α + β) Γ(α)Γ(β)

xα−1 (1− x)β−1, 0 < x < 1.

E(X) = α

α + β e V (X) =

αβ

(α + β)2(α + β + 1) .

A.9 Distribuição de Dirichlet

O vetor aleatórioX = (X1, . . . , Xk) tem distribuição de Dirichlet com parâmetros

α1, . . . , αk, denotada por Dk(α1, . . . , αk) se sua função de densidade conjunta é

dada por

p(x|α1, . . . , αk) = Γ(α0)

Γ(α1), . . . ,Γ(αk) xα1−11 . . . x

αk−1 k ,

k ∑

i=1

xi = 1,

para α1, . . . , αk > 0 e α0 = ∑k

i=1 αi.

E(Xi) = αi α0

, V (Xi) = (α0 − αi)αi α20(α0 + 1)

, e Cov(Xi, Xj) = − αiαj

α20(α0 + 1)

Note que a distribuição Beta é obtida como caso particular para k = 2.

A.10 Distribuição t de Student

X tem distribuição t de Student (ou simplesmente t) com parâmetros µ ∈ R, σ2 > 0 e ν > 0 (chamado graus de liberdade), denotando-se X ∼ tν(µ, σ2), se sua função de densidade é dada por

p(x|ν, µ, σ2) = Γ( ν+1 2 ) νν/2

Γ(ν 2 ) √ πσ

[

ν + (x− µ)2

σ2

]−(ν+1)/2

, x ∈ R.

E(X) = µ, para ν > 1 e V (X) = σ2 ν

ν − 2 , para ν > 2.

Um caso particular da distribuição t é a distribuição de Cauchy, denotada por

C(µ, σ2), que corresponde a ν = 1.

A.11. DISTRIBUIÇÃO F DE FISHER 97

A.11 Distribuição F de Fisher

X tem distribuição F com ν1 > 0 e ν2 > 0 graus de liberdade, denotando-se

X ∼ F (ν1, ν2), se sua função de densidade é dada por

p(x|ν1, ν2) = Γ(ν1+ν2

2 )

Γ(ν1 2 )Γ(ν2

2 ) ν ν1/2 1 ν

ν2/2 2 x

ν1/2−1(ν2 + ν1x) −(ν1+ν2)/2, x > 0.

E(X) = ν2

ν2 − 2 , para ν2 > 2 e V (X) =

2ν22(ν1 + ν2 − 2) ν1(ν2 − 4)(ν2 − 2)2

, para ν2 > 4.

A.12 Distribuição de Pareto

X tem distribuição de Pareto com parâmetros α e β denotando-se X ∼ Pareto(α, β), se sua função de densidade de probabilidade é dada por,

p(x|α, β) = α β

(

β

x

)α+1

, x > β.

E(X) = αβ

α− 1 e V (X) = αβ2

α− 2 − (

αβ

α− 1

)2

.

A.13 Distribuição Binomial

X tem distribuição binomial com parâmetros n ≥ 1 e p ∈ (0, 1), denotando-se X ∼ bin(n, p), se sua função de probabilidade é dada por

p(x|n, p) = (

n

x

)

px(1− p)n−x, x = 0, . . . , n.

E(X) = np e V (X) = np(1− p)

e um caso particular é a distribuição de Bernoulli com n = 1.

A.14 Distribuição Multinomial

O vetor aleatório X = (X1, . . . , Xk) tem distribuição multinomial com parâmet-

ros n e probabilidades θ1, . . . , θk, denotada por Mk(n, θ1, . . . , θk) se sua função de

probabilidade conjunta é dada por

p(x|θ1, . . . , θk) = n!

x1!, . . . , xk! θx11 , . . . , θ

xk k , xi = 0, . . . , n,

k ∑

i=1

xi = n,

98 APÊNDICE A. LISTA DE DISTRIBUIÇÕES

para 0 < θi < 1 e ∑k

i=1 θi = 1. Note que a distribuição binomial é um caso

particular da distribuição multinomial quando k = 2. Além disso, a distribuição

marginal de cada Xi é binomial com parâmetros n e θi e

E(Xi) = nθi, V (Xi) = nθi(1− θi), e Cov(Xi, Xj) = −nθiθj.

A.15 Distribuição de Poisson

X tem distribuição de Poisson com parâmetro θ > 0, denotando-se X ∼ Poisson(θ), se sua função de probabilidade é dada por

p(x|θ) = θ xe−θ

x! , x = 0, 1, . . .

E(X) = V (X) = θ.

A.16 Distribuição Binomial Negativa

X tem distribuição de binomial negativa com parâmetros r ≥ 1 e p ∈ (0, 1), denotando-se X ∼ BN(r, p), se sua função de probabilidade é dada por

p(x|r, p) = (

r + x− 1 x

)

pr(1− p)x, x = 0, 1, . . .

E(X) = r(1− p)

p e V (X) =

r(1− p) p2

.

Um caso particular é quando r = 1 e neste caso diz-se que X tem distribuição

geométrica com parâmetro p. Neste caso,

p(x|p) = pr(1− p)x, x = 0, 1, . . .

E(X) = 1− p p

e V (X) = 1− p p2

.

Apêndice B

Alguns Endereços Interessantes

Neste apêndice são listados alguns endereços na internet com conteúdo relativo a

abordagem Bayesiana.

ˆ Teorema de Bayes no Wikipedia: http://en.wikipedia.org/wiki/Bayes the-

orem

ˆ Bayesian Analysis - The Journal: http://ba.stat.cmu.edu/

ˆ International Society for Bayesian Analysis: http://www.bayesian.org

ˆ American Statistical Association, Section on Bayesian Statistical Science:

http://www.amstat.org/sections/SBSS

ˆ Bayes Methods Working Group of the International Biometric Society, Ger-

man Region: http://ibealt.web.med.uni-muenchen.de/bayes-ag

ˆ Encontro Brasileiro de Estat́ıstica Bayesiana:

2006 (http://www.im.ufrj.br/ebeb8),

2008 (http://www.ime.usp.br/˜ isbra/ebeb/9ebeb)

ˆ Valencia Meetings: http://www.uv.es/valenciameeting

ˆ I Workshop em Estat́ıstica Espacial e Métodos Computacionalmente Inten-

sivos: leg.ufpr.br/˜ ehlers/folder

ˆ Case Studies in Bayesian Statistics: http://lib.stat.cmu.edu/bayesworkshop/

ˆ MCMC preprints: http://www.statslab.cam.ac.uk/˜ mcmc

ˆ Projeto BUGS (Bayesian inference Using Gibbs Sampling):

http://www.mrc-bsu.cam.ac.uk/bugs

ˆ Projeto JAGS (Just Another Gibbs Sampler):

http://www-fis.iarc.fr/˜ martyn/software/jags/

99

100 APÊNDICE B. ALGUNS ENDEREÇOS INTERESSANTES

ˆ BayesX (Bayesian Inference in Structured Additive Regression Models.):

http://www.stat.uni-muenchen.de/˜ bayesx/bayesx.html

ˆ MrBayes (Bayesian estimation of phylogeny): http://mrbayes.scs.fsu.edu

ˆ Número especial do Rnews dedicado a inferencia Bayesiana e MCMC:

http://www.est.ufpr.br/R/doc/Rnews/Rnews 2006-1.pdf

ˆ CRAN Task View (Bayesian Inference):

http://cran.r-project.org/src/contrib/Views/Bayesian.html

ˆ Centro de Estudos do Risco UFSCAR:

http://www.ufscar.br/˜ des/CER/inicial.htm

Referências

Berger, J. (1985). Statistical Decision Theory and Bayesian Analysis. Springer-

Verlag: New York.

Bernardo, J. M. & Smith, A. F. M. (1994). Bayesian Theory. Wiley: New York.

Box, G. E. P. & Tiao, G. C. (1992). Bayesian Inference in Statistical Analysis.

Wiley Classics Library ed. Wiley-Interscience.

Broemeling, L. (1985). Bayesian Analysis of Linear Models. New York: Marcel

Dekker.

DeGroot, M. H. (1970). Optimal Statistical Decisions. McGraw-Hill Book Co.

Evans, M., Hastings, N. & Peacock, B. (1993). Statistical Distributions, Second

Edition (Second ed.). Wiley Interscience.

Gamerman, D. (1997). Markov chain Monte Carlo: Stochastic Simulation for

Bayesian Inference. Texts in Statistical Sciences. Chapman and Hall, Lon-

don.

Gamerman, D. & Lopes, H. (2006). Markov chain Monte Carlo: Stochastic

Simulation for Bayesian Inference. Texts in Statistical Science Series. CRC

Press.

Gelman, A., Carlin, J. B., Stern, H. S. & Rubin, D. B. (2004). Bayesian Data

Analysis (2nd ed.). Chapman and Hall: London.

Green, P. J. (1995). Reversible jump MCMC computation and Bayesian model

determination. Biometrika 82, 711–732.

Johnson, N. L., Kotz, S. & Balakrishnan, N. (1995). Continuous Univariate

Distributions (2nd ed.), Volume 2. John Wiley, New York.

Johnson, N. L., Kotz, S. & Kemp, A. W. (1992). Univariate Discrete Distri-

butions (2nd ed.). John Wiley, New York.

Migon, H. S. & Gamerman, D. (1999). Statistical Inference: An Integrated

Approach. Arnold.

O’Hagan, A. (1994). Bayesian Inference, Volume 2B. Edward Arnold, Cam-

bridge.

101

102 References.

Robert, C. P. & Casella, G. (1999). Monte Carlo Statistical Methods. Springer-

Verlag, New York.

Smith, A. F. M. & Gelfand, A. E. (1992). Bayesian statistics without tears: A

sampling-resampling perspective. The American Statistician 46, 84–88.

Até o momento nenhum comentário
Esta é apenas uma pré-visualização
3 mostrados em 106 páginas