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

Baixe Apostila Inferência Bayesiana - Ricardo Ehlers e outras Notas de estudo em PDF para Estatística, somente na Docsity! 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 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 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) ∣ ∣ ∣ ∣ dδ dθ ∣ ∣ ∣ ∣ = (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.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. 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 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=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 . 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 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(θ(φ)) ∣ ∣ ∣ ∣ dθ dφ ∣ ∣ ∣ ∣ ∝ ∣ ∣ ∣ ∣ dθ dφ ∣ ∣ ∣ ∣ . 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.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. 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=1Xi α + β + 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 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 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.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. 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 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; 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 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 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 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 + } + } 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 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