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

Baixe Apostila Inferência Bayesiana e outras Notas de estudo em PDF para Estatística, somente na Docsity! INTRODUÇÃO À INFERÊNCIA BAYESIANA RICARDO S. EHLERS Laboratório de Estat́ıstica e Geoinformação Universidade Federal do Paraná Primeira publicação 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 c© RICARDO SANDES EHLERS 2003-2007 Prefácio O objetivo principal deste texto é oferecer um material didático básico para um curso introdutório de Inferência Bayesiana a ńıvel de graduação. Ele pode ser adotado em cursos de Bacharelado em Estat́ıstica bem como em outros cursos de graduação e de pós-graduação aonde os alunos tenham conhecimentos básicos de probabilidade e cálculo. Algum conhecimento sobre estimação de máxima verossimilhança também é útil porém não essencial. O texto se originou de notas de aulas da disciplina de Inferência Bayesiana minis- trada no programa de Bacharelado em Estat́ıstica da Universidade Federal do Paraná. A idéia é apresentar o enfoque Bayesiano como alternativa à abordagem clássica estabelecendo algumas comparações inevitáveis. O texto não se propõe a ser exaustivo nem deve ser visto como um livro de receitas com soluções Bayesianas para problemas de análise de dados. O manuscrito foi preparado usando o LATEX e todas as ilustrações foram produzi- das no pacote estat́ıstico R (gratuito e de código aberto) que pode ser obtido em http://www.r-project.org/ Em vários exemplos são fornecidos também os comandos do R que foram utilizados e mostradas as sáıdas resultantes de modo que o leitor é encorajado a reproduzi-los. Este texto certamente não está livre de erros, e comentários e sugestões dos leito- res são bem vindos. Citar este texto como: Ehlers, R.S. (2007) Introdução à Inferência Bayesiana. Dispońıvel em http://leg.ufpr.br/~ ehlers/bayes. Acesso em: ... Ricardo S. Ehlers Curitiba, novembro de 2007. i 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 assu- mir diferentes graus. Do ponto de vista Bayesiano, estes diferentes graus de incerteza são representados através de modelos probabiĺısticos para θ. Neste contexto, é na- tural 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 consi- derados quantidades aleatórias. 1.1 Teorema de Bayes Considere uma quantidade de interesse desconhecida θ (tipicamente não observá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 normalizadora 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 verossimilhanç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) 1 2 CAPÍTULO 1. INTRODUÇÃO 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 substitúıda por uma proporcionalidade. Esta forma simplificada do teorema de Bayes será útil em pro- blemas que envolvam estimação de parâmetros já que o denominador é apenas uma constante normalizadora. Em outras situações, como seleçã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 constante normalizadora da posteriori pode ser facilmente recuperada pois p(θ|x) = kp(x|θ)p(θ) onde k−1 = ∫ p(x|θ)p(θ)dθ = Eθ[p(X|θ)] = p(x) 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 geoestat́ıstica) o maior interesse é 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 quanti- dade 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 estat́ısticos a hipótese de independência condicional entre 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 anaĺıtica e precisará ser obtida por algum método de aproximação. 1.1. TEOREMA DE BAYES 3 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 pessoa, “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 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 6 CAPÍTULO 1. INTRODUÇÃO Na Figura 1.1(b) esta função de densidade está representada para alguns valores de µ e σ2. Novamente note como informações a priori bastante diferentes podem ser representadas. Finalmente, embora existam outras possibilidades, vamos atribuir uma distri- buição a priori θ ∼ Beta(a, b) i.e. (ver Apêndice A), p(θ) ∝ θa−1(1 − θ)b−1, a, 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 famı́lia de distribuições a priori para θ. Algumas possibilidades estão representadas na Figura 1.1(c). 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 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 ) onde µ1 = τ−20 µ0 + σ −2x τ−20 + σ −2 e τ−21 = τ −2 0 + σ −2. 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 podemos reescrever as informações na forma de equações com erros não correlacionados. Assim, X = θ + ǫ, ǫ ∼ N(0, σ2) θ = µ0 + w, w ∼ N(0, τ20 ) tal que Cov(θ, ǫ) = Cov(µ,w) = 0. Portanto a distribuição (incondicional) de X é normal pois ele resulta de uma soma de variáveis aleatórias com distribuição normal. Além disso, E(X) = E(θ) + E(ǫ) = µ0 V ar(X) = V ar(θ) + V ar(ǫ) = τ20 + σ 2 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) (a) 0.0 0.2 0.4 0.6 0.8 1.0 0 1 2 3 4 θ p(θ ) N(−1,0.5) N(1,1) N(0,4) (b) 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) (c) Figura 1.1: Densidades a priori para o parâmetro θ no Exemplo 1.2. (a) Normal truncada, (b) transformação loǵıstica e (c) Beta(a, b). 8 CAPÍTULO 1. INTRODUÇÃO Conclusão, X ∼ N(µ0, τ20 + σ2). Exemplo 1.3 : (Box & Tiao, 1992) Os f́ısicos A e B desejam determinar uma cons- tante 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 é dif́ıcil verificar que para o f́ısico A: P (860 < θ < 940) ≈ 0, 95 para o f́ısico 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 τ−21 = 0, 00312 (au- mento de 25%). • para o f́ısico B: precisão(θ) passou de τ−20 = 0, 000156 para τ−21 = 0, 000781 (aumento de 400%). A situação está representada graficamente na Figura 1.2 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 f́ısico A enquanto que a incerteza do f́ısico B foi bastante reduzida. Os comandos do R abaixo podem ser usados nos cálculos. 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)) } 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 classificado 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: 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 re- presentar (probabilisticamente) o conhecimento que se tem sobre θ antes da realização do experimento. Neste caṕıtulo 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 famı́lia 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 pri- oris 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 indiscri- minada de prioris conjugadas. Essencialmente, o problema é que a priori 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. 11 12 CAPÍTULO 2. DISTRIBUIÇÕES A PRIORI (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 onde 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 distribuiçõ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 Estat́ıstica, tanto cont́ınuas quanto discretas. Uma caracteŕıstica essencial desta famı́lia é que existe uma estat́ıstica suficiente com dimensão fixa. Ve- remos adiante que a classe conjugada de distribuições é muito fácil de caracterizar. Definição 2.2 A famı́lia de distribuições com função de (densidade) de probabilidade p(x|θ) pertence à famı́lia 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 estat́ıstica suficiente para θ. Neste caso, a classe conjugada é facilmente identificada como, p(θ) = k(α, β) exp{αφ(θ) + βb(θ)}. 2.2. CONJUGAÇÃO NA FAMÍLIA EXPONENCIAL 13 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 estensã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 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 prioris conju- gadas Beta(1,1), Beta(2,2) e Beta(1,3). As funções de densidade desta distribuições juntamente com a função de verossimilhança normalizada e as respectivas 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. No Exemplo 2.2 suponha novamente que n = 12, X = 9 e usamos as prioris conju- gadas 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. 16 CAPÍTULO 2. DISTRIBUIÇÕES A PRIORI Finalmente, a definição de famı́lia exponencial pode ser extendida ao caso multi- paramétrico, i.e. p(x|θ) = [ n ∏ i=1 a(xi) ] exp    r ∑ j=1 [ n ∑ i=1 uj(xi) ] φj(θ) + nb(θ)    onde θ = (θ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 verssimilhanç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 proporcio- nalidade. 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 o Teorema 1.1 com as devidas substituições, i.e. a distribuição a posteriori de θ dado x é N(µ1, τ 2 1 ) onde µ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 opcio- nalmente fazer os gráficos das densidades. 2.3. PRINCIPAIS FAMÍLIAS CONJUGADAS 17 norm.norm = function(x,sigma,mu0,tau0,plot=F){ n = length(x) xbar = mean(x) ep = sigma/sqrt(n) sigma2 = sigma**2 precisao = n*(1/sigma2)+(1/tau0) mu1 = (n*(1/sigma2)*xbar+(1/tau0)*mu0)/precisao 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(precisao)),add=T,col=3) } } 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 distri- buições Gama que é fechada por amostragem. Assim, a priori conjugada natural de θ é Gama com parâmetros positivos α e β, i.e. p(θ) ∝ θα−1e−βθ, α, β > 0 θ > 0. A densidade a posteriori fica p(θ|x) ∝ θα+t−1 exp {−(β + n)θ} que corresponde à densidade Gama(α+ t, β + n). 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 . 18 CAPÍTULO 2. DISTRIBUIÇÕES A PRIORI 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 catego- rias 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 k − 1 parâmetros já que temos a seguinte restrição ∑pi=1 θi = 1. Além disso, a restrição ∑p i=1Xi = 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 . Note que esta é uma generalização da distribuição binomial que apenas duas catego- rias. Não é dif́ıcil mostrar que esta distribuição também pertence à famı́lia exponen- cial. 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ância desco- nhecida 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}. 2.4. PRIORI NÃO INFORMATIVA 21 (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 “ignorâ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 pensada como um ponto de partida quando não se consegue fazer uma elicitação detalhada do “verda- deiro” conhecimento a priori. Neste sentido, serão apresentadas aqui algumas formas de “como” fazer enquanto discussões mais detalhadas são encontradas em Berger (1985), Box e Tiao (1992), Bernardo e 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 distribuiçã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θ = ∞. (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(θ(φ)) ∣ ∣ ∣ ∣ dθ dφ ∣ ∣ ∣ ∣ ∝ ∣ ∣ ∣ ∣ dθ dφ ∣ ∣ ∣ ∣ . Na prática, como estaremos interessados na distribuição a posteriori não daremos muita importância à impropriedade da distribuição a priori. No entanto devemos 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. 22 CAPÍTULO 2. DISTRIBUIÇÕES A PRIORI 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 cur- vatura mais precisa é a informação contida na verossimilhança, ou equivalentemente 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 probabilidade p(x|θ). A priori não informativa de Jeffreys tem função de densidade dada por p(θ) ∝ [I(θ)]1/2. 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. 2.4. PRIORI NÃO INFORMATIVA 23 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 conveniente- mente. 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 quanti- dade θ 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 exemplos importantes são a distribuição normal com variância conhecida, e a distribuição nor- mal 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 quan- tidade σ 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. 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 informa- tiva. 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. 26 CAPÍTULO 2. DISTRIBUIÇÕES A PRIORI 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. 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) 2.6 Problemas 1. Mostre que a famı́lia de distribuições Beta é conjugada em relação às distri- buições amostrais binomial, geométrica e binomial negativa. 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 θ conhe- cido. 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 (θ, φ). 2.6. PROBLEMAS 27 (a) Determine os parâmetros (µ0, c0, n0, σ 2 0) utilizando as seguintes informações a priori: E(θ) = 0, P (|θ| < 1, 412) = 0, 5, E(φ) = 2 e E(φ2) = 5. (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 esboce os gráficos das distribuições a priori, a posteriori e da função de verossimi- lhança, com φ fixo. (c) Calcule P (|Y | > 1|x) onde Y é uma observação tomada da mesma po- pulação. 5. Suponha que o tempo, em minutos, para atendimento a clientes segue uma dis- tribuição exponencial com parâmetro θ desconhecido. Com base na experiê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. SejaX1, . . . , 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 posteriori 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 é desconhe- cido 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. SejaX1, . . . , 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. 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. 28 CAPÍTULO 2. DISTRIBUIÇÕES A PRIORI Mostre que esta famı́lia é conjugada ao modelo normal com média µ conhecida e variância θ desconhecida. 11. Suponha que X = (X1, X2, X3) tenha distribuição trinomial com parâmetros 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âmetros 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 obtidas naquele exemplo obtenha a distribuição a posteriori de θ (incondicional). Esboce e comente os gráficos das densidades a priori e posteriori. 17. Se X ∼ Binomial Negativa(v, θ) obtenha a priori de Jeffreys para θ. 18. Se X ∼ Geometrica(θ) obtenha a priori de Jeffreys para θ. 3.2. ESTIMADORES DE BAYES 31 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 = ∑n i=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=1Xi α+ β + n . A perda quadrática é as vezes criticada por penalizar demais o erro de estimaçã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 conside- rar 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. 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 . 32 CAPÍTULO 3. ESTIMAÇÃO 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 anterior, então δ(x) = 0, 107. Ou seja, as estimativas de Bayes são bastante próximas, 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 informaçã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 distribuição a posteriori e uma estimativa pontual. Será discutido o conceito de intervalo de cre- dibilidade (ou intervalo de confiança Bayesiano) baseado no distribuição a posteriori. Definição 3.3 C é um intervalo de credibilidade de 100(1-α)%, ou ńıvel de credibi- lidade (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 trans- formaçõ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 3.4. ESTIMAÇÃO NO MODELO NORMAL 33 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áxima 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 densidade 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 trans- formaçõ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 nor- mal para estimação da média e variância em problemas de uma ou mais amostras e em modelos de regressão linear. A análise será feita com priori conjugada e priori não informativa quando serão apontadas as semelhanças com a análise clássica. As- sim como nos caṕıtulos anteriores a abordagem aqui é introdutória. Um tratamento mais completo do enfoque Bayesiano em modelos lineares pode ser encontrado em Broemeling (1985) e Box e Tiao (1992). Nesta seção considere uma amostra aleatória X1, · · · , Xn tomada da distribuiçã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. 36 CAPÍTULO 3. ESTIMAÇÃO 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 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 quere- mos estimar β = θ1−θ2 (por exemplo, testar se θ1 = theta2). 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 inde- pendentes, 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 (θ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 )) 3.4. ESTIMAÇÃO NO MODELO NORMAL 37 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 ) 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 estima- tiva 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 conju- gada com ci = 0, σ 2 0 = 0 e n0 = −2. Assim, temos que c∗i = ni, µ∗i = xi, n∗0 = ν e 38 CAPÍTULO 3. ESTIMAÇÃO n∗0σ ∗2 0 = νs 2 e a estimativa pontual concide com a estimativa de máxima verossimi- lhanç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 su- posiçã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) 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). 4.2. O PROBLEMA GERAL DA INFERÊNCIA BAYESIANA 41 4.2 O Problema Geral da Inferência Bayesiana A distribuição a posteriori pode ser convenientemente resumida em termos de espe- ranç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, i.e. p(θ1|x) = ∫ p(θ|x)dθ2 onde θ = (θ1,θ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 = var(θ) = 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 multidimensi- onais, é extremamente importante em inferência Bayesiana. Inferência exata somente será posśıvel se estas integrais puderem ser calculadas analiticamente, caso contrário devemos usar aproximações. Nas próximas seções iremos apresentar métodos aproxi- mados baseados em simulação para obtenção dessas integrais. 4.3 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θ. 42 CAPÍTULO 4. COMPUTAÇÃO BAYESIANA 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, trans- formamos o problema de avaliar a integral no problema estat́ıstico de estimar 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 valo- res 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 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){ # Calcula a integral de exp(-x) no intervalo (a,b) x = runif(n,a,b) y = exp(-x) int.exp = (b-a)*mean(y) return(int.exp) } 4.3. MÉTODO DE MONTE CARLO SIMPLES 43 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. Isto se chama “erro de Monte Carlo” e decresce conforme aumentamos o número de simulações. Repetindo o expe- rimento com n = 1000 a variação ficará bem menor. Na Figura 4.1 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) 20 50 100 200 500 0. 25 0. 30 0. 35 0. 40 Figura 4.1: Boxplots para 50 estimativas da integral no Exemplo 4.1 com n=20, 50, 100, 200, e 500 simulações. A generalização é bem simples para o caso em que a integral é a esperança ma- temá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) 46 CAPÍTULO 4. COMPUTAÇÃO BAYESIANA 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(θ). Exemplo 4.3 : Para uma única observação X com distribuição N(θ, 1), θ desconhe- cido, e priori Cauchy(0,1) segue que 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 estimativa. 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 Bayesia- nas. 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 integrais. Neste caso, a apro- ximação é baseada na razão dos dois estimadores de Monte Carlo para o numerador e denominador. Exerćıcios 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 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. 4.4. MÉTODOS DE REAMOSTRAGEM 47 2. Para uma única observação X com distribuição N(θ, 1), θ desconhecido, quere- mos 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çõesX1, . . . , Xk da distribuiçãoN(θ, 1). 4.4 Métodos de Reamostragem Existem distribuições para as quais é muito dif́ıcil ou mesmo imposśıvel simular va- lores. 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 to- mar a priori como distribuição auxiliar conforme proposto em Smith e Gelfand (1992) . 4.4.1 Método de Rejeição Considere uma densidade auxiliar q(θ) da qual sabemos gerar valores. A única res- triçã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 va- lor como sendo da distribuição a posteriori com probabilidade p(θ|x)/Aq(θ). Caso contrário, θ∗ não é aceito como uma valor gerado da posteriori 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(θ). 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 item 1. 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. 48 CAPÍTULO 4. COMPUTAÇÃO BAYESIANA Outro problema é que a taxa de aceitação pode ser muito baixa, i.e. teremos que gerar muitos valores da distribuição auxiliar até conseguir um número suficiente de valores da posteriori. Isto ocorrerá se as informações da priori e da verossimilhança forem conflitantes já que neste caso os valores gerados terão baixa probabilidade de serem aceitos. Exemplo 4.4 : Suponha que X1, . . . , Xn ∼ N(θ, 1) e assume-se uma priori Cau- chy(0,1) para θ. A função de verossimilhança é p(x|θ) ∝ exp[−n(x̄− θ)2/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 abaixo obtém uma amostra de tamanho m de θ e como ilustração vamos gerar 50 observações da N(2,1). rej = function(x,m){ x.bar = mean(x) n = length(x) # theta = rcauchy(m,2,1) theta = rcauchy(m,0,1) u = runif(m,0,1) peso = exp(-0.5*n*(theta-x.bar)**2) acc.theta = theta[u<peso] acc=mean(u<peso) cat(’\nTaxa de aceitacao’,acc,’\n’) return(list(acc=acc,theta=theta,acc.theta=acc.theta)) } x=rnorm(50,2,1) m=rej(x,1000) Taxa de aceitacao 0.018 Note que a taxa de aceitação é extremamente baixa, somente 18 dentre 1000 valores de θ foram aceitos o que constitui uma amostra muito pequena. Isto ocorreu devido ao conflito entre verossimilhança e priori. O problema é ilustrado na Figura 4.3 onde se pode notar que a maioria dos valores de θ foi gerada em regiões de baixa verossimilhança. 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. 4.4.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. 4.5. MONTE CARLO VIA CADEIAS DE MARKOV 51 tematicamente, 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; • 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 influenciam 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 converge 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.5.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 sequê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 ĝ é então estimado como √ √ √ √ 1 k(k − 1) k ∑ i=1 (Bi −B)2 52 CAPÍTULO 4. COMPUTAÇÃO BAYESIANA 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.5.3 Algoritmo de Metropolis-Hastings Os algoritmos de Metropolis-Hastings usam a mesma idéia dos métodos de rejeição vistos no caṕıtulo anterior, i.e. um valor é gerado de uma distribuição auxiliar e aceito com uma dada probabilidade. Este mecanismo de correção garante que 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 completamente 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) = θ. 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 aplicaçõ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. 4.5. MONTE CARLO VIA CADEIAS DE MARKOV 53 Exemplo 4.5 : 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 + θ 2 , 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 uma priori θ ∼ U(0, 1) segue que a posteriori é proporcional à expressão acima. Tomando a distribuição U(0, 1) como proposta então q(θ) = 1, ∀ θ e a probabilidade (4.3) se simplifica para α(θ, θ′) = min { 1, p(x|θ′) p(x|θ) } = min { 1, ( 2 + θ′ 2 + θ )y1 (1 − θ′ 1 − θ )y2+y3 (θ′ θ )y4 } . 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 1000 valores de θ. Os valores simulados e as primeiras 30 autocorrelações amostrais de θ estão na Figura 4.4. 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 também 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. Os resultados foram obtidos usando os comandos do R a seguir. p = function(x,y) (2+x)^y[1] * (1-x)^(y[2]+y[3]) * x^y[4] metr = function(n,y,p,start){ theta = matrix(NA, nrow=n) theta[1] = start for (i in 2:n) { x = runif(1) A = p(x,y)/p(theta[i-1],y) prob = min(1,A) u = runif(1) taxa= 0 if (u < prob) { theta[i] = x taxa = taxa + 1 } else theta[i] = theta[i-1] 56 CAPÍTULO 4. COMPUTAÇÃO BAYESIANA 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 blocos {θ1, . . . ,θk} e dentro de cada iteração teremos o algoritmo aplicado k vezes. Definindo o vetor x−i = (x1, . . . ,xi−1,xi+1, . . . ,xk) que contém todos os elementos de x exceto xi suponha que na iteração t+ 1 os blocos 1, 2, . . . , i− 1 já foram atualizados, i.e. x−i = (x (t+1) 1 , . . . ,x (t+1) i−1 ,x (t) i+1, . . . ,x (t) k ). Para atualizar a i-ésima componente, um valor de xi é gerado da distribuição proposta q(·|xi,x−i) e este valor candidato é aceito com probabilidade α(xi,x ′ i) = min { 1, π(x′i|x−i) q(xi|x′i,x−i)) π(xi|x−i) q(x′i|xi,x−i) } . (4.5) Aqui, π(xi|x−i) é chamada de distribuição condicional completa como será visto na próxima seção. Exerćıcios 1. Assumindo que a distribuição estacionária é N(0, 1), (a) faça 500 iterações do algoritmo de Metropolis com distribuições propostas 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 experimentalmente 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.5. MONTE CARLO VIA CADEIAS DE MARKOV 57 0 100 200 300 400 500 −2 −1 0 1 2 3 (a) 0 100 200 300 400 500 −2 −1 0 1 2 (b) Figura 4.5: 1000 valores simulados para o Exemplo 4.6 usando o algoritmo de Metropolis- Hastings com (a) σ = 0.5 e (b) σ = 10. 58 CAPÍTULO 4. COMPUTAÇÃO BAYESIANA 4.5.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. Por- tanto, a distribuição condicional completa é a distribuição da i-ésima componente de θ condicionada em todas as outras componentes. Ela é obtida a partir da distribuição conjunta como, π(θi|θ−i) = π(θ) ∫ π(θ)dθi . 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.7 : 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.6. PROBLEMAS DE DIMENSÃO VARIÁVEL 61 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) } print(round(table(m[nburn+1:n])/(n-nburn),6)) return(list(lambda=lambda, phi=phi, m=m)) } 4.6 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 esti- mado. 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 indica- dor 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). 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. 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) 62 CAPÍTULO 4. COMPUTAÇÃO BAYESIANA 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 estados 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 dis- tribuições com dimensão fixa. Embora existam outras possibilidades iremos estudar os algoritmos MCMC com saltos reverśıveis (Green 1995) que são particularmente úteis no contexto de seleção Bayesiana de modelos. 4.6.1 MCMC com Saltos Reverśıveis (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 en- tanto, os movimentos podem ser entre espaços de dimensõ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 procedimento. Suponha que o estado atual é (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.10 : 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 a priori p(k) = 1/2 para o indicador de modelo e as seguintes prioris 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). 4.6. PROBLEMAS DE DIMENSÃO VARIÁVEL 63 Dado o modelo, as funções de verossimilhança ficam p(y|λ, k = 1) = λne−λ P yi p(y|α, β, k = 2) = β nα Γn(α) ∏ yα−1i e −β P yi as 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 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 } # rejeita o novo valor 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) # prob. aceitacao u = runif(1) if (u < aceita) newalpha = z else newalpha = alpha 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. 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) 66 CAPÍTULO 4. COMPUTAÇÃO BAYESIANA 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 com o conjunto de dados que está no arquivo gama1.dat. y = scan(’http://leg.est.ufpr.br/ce718/gama1.dat’) curve(dgamma(x,10,10),from=0,to=5) curve(dgamma(x,1,1),from=0,to=5,add=T) curve(dgamma(x,3,3),from=0,to=5,add=T) m = rjmcmc(1000,500,y,10,1,1) Probabilidades a posteriori dos modelos [1] 0.718 0.282 Medias a posteriori dos parametros [1] 1.247682 1.197084 1.516506 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] # medias a posteriori dos tempos medios de vida x = m$x1[,c(1,2)] x[r1,1] = 1/m$x1[r1,1] 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]) medias [1] 0.8672224 0.8244537 # Calculando uma estimativa ponderada do tempo medio de vida prob = m$nv1/(niter-nburn) prob[1]*medias[1] + prob[2]*medias[2] [1] 0.8551617 4.7 Tópicos Relacionados 4.7.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 4.7. TÓPICOS RELACIONADOS 67 situações estes valores podem ser altamente correlacionados e em geral a autocor- relaçã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 desperd́ıcio. Métodos de séries temporais estão dispońıveis 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.7.2 Monitorando a Convergência Aqui vale lembrar que a verificação de convergência (ou falta de convergência) é res- ponsabilidade do analista. Além disso estamos falando de convergência para a distri- buição alvo, que neste caso é a distribuição a posteriori, o que pode ser extremamente dif́ıcil de se verificar na prática. 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ão apresentadas suas funções de (densidade) de probabili- dade além da média e variância. Uma revisão exaustiva de distribuições de probabili- dades pode ser encontrada em Johnson et al. (1994), Johnson et al. (1995) e Johnson et al. (1992). A.1 Distribuição Normal X tem distribuição normal com parâmetros µ e σ2, denotando-se X ∼ N(µ, σ2), se sua função de densidade é dada por p(x|µ, σ2) = (2πσ2)−1/2 exp[−(x− µ)2/2σ2], −∞ < x <∞, para −∞ < µ < ∞ e σ2 > 0. Quando µ = 0 e σ2 = 1 a distribuição é chamada normal padrão. A distribuição log-normal é definida como a distribuição de eX . No caso vetorial, X = (X1, . . . , Xp) tem distribuição normal multivariada 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. A.2 Distribuição Gama X tem distribuição Gama com parâmetros α e β, denotando-se X ∼ Ga(α, β), se sua função de densidade é dada por p(x|α, β) = β α Γ(α) xα−1e−βx, x > 0, para α, β > 0. E(X) = α/β e V (X) = α/β2. 68 A.10. DISTRIBUIÇÃO BINOMIAL 71 A.10 Distribuição Binomial X tem distribuição binomial com parâmetros n e p, 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 para n ≥ 1 e 0 < p < 1. E(X) = np e V (X) = np(1 − p) e um caso particular é a distribuição de Bernoulli com n = 1. A.11 Distribuição Multinomial O vetor aleatório X = (X1, . . . , Xk) tem distribuição multinomial com parâmetros n e probabilidades θ1, . . . , θk, denotada por Mk(n, θ1, . . . , θk) se sua função de probabi- lidade conjunta é dada por p(x|θ1, . . . , θk) = n! x1!, . . . , xk! θx11 , . . . , θ xk k , xi = 0, . . . , n, k ∑ i=1 xi = n, para 0 < θi < 1 e ∑k i=1 θi = 1. Note que a distribuição binomial é um caso especial da 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.12 Distribuição de Poisson X tem distribuição de Poisson com parâmetro θ, denotando-se X ∼ Poisson(θ), se sua função de probabilidade é dada por p(x|θ) = θ xe−θ x! , x = 0, 1, . . . para θ > 0. E(X) = V (X) = θ. A.13 Distribuição Binomial Negativa X tem distribuição de binomial negativa com parâmetros r e p, 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, . . . 72 APÊNDICE A. LISTA DE DISTRIBUIÇÕES para r ≥ 1 e 0 < p < 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. 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 theorem • 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, German 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 Intensivos: 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/ • BayesX (Bayesian Inference in Structured Additive Regression Models.): http://www.stat.uni-muenchen.de/˜ bayesx/bayesx.html 73