Aplicações de Approximate Bayesian Computation a controle de

Aplicaes deApproximate Bayesian Computation

a controle de qualidade

Thiago Feitosa Campos

Tese apresentadaao

Instituto de Matemtica e Estatsticada

Universidade de So Paulopara

obteno do ttulode

Doutor em Cincias

Programa: EstatsticaOrientador: Prof. Dr. Sergio Wechsler

Coorientador: Prof. Dr. Lus Gustavo Esteves

Durante o desenvolvimento deste trabalho o autor recebeu auxlio financeiro da CAPES

So Paulo, Junho de 2015

Aplicaes de Approximate Bayesian Computation a controle dequalidade

Esta verso da tese contm as correes e alteraes sugeridaspela Comisso Julgadora durante a defesa da verso original do trabalho,realizada em 11/06/2015. Uma cpia da verso original est disponvel no

Instituto de Matemtica e Estatstica da Universidade de So Paulo.

Comisso Julgadora:

Prof. Dr. Sergio Wechsler (orientador) - IME-USP

Prof. Dr. Lus Gustavo Esteves - IME-USP

Prof. Dr. Jess Enrique Garca - UNICAMP

Profa. Dra. Laura Leticia Ramos Rifo - UNICAMP

Prof. Dr. Marcio Alves Diniz - UFSCar

Agradecimentos

A tarefa de escrever essa tese foi bastante difcil, ceratocone nos dois olhos Foi um fator com-plicador que tive que lidar durante a produo dessa tese. No teria conseguido sem o apoio deamigos e famlia, em especial minha me que sempre deu apoio incondicional durante toda minhavida. Meus amigos Alexandre Patriota e Victor Fossaluza, atualmente professores do IME-USP, eCaio Freire foram importantes, sempre ajudando a manter a tempestade mental viva.

O perodo de tempo que passei como professor substituto na UFSCar foi basante engrandecedorpara minha vida, e por isso quero agradecer a todos os professores e funcionrios do departamentode estatstica da universidade federal de So Carlos apenas por terem feito parte da minha vidanaquele momento.

Minha amiga Raquel Guimares e sua tia Catharina merecem um agradecimento especial porme ajudarem a superar momentos difceis.

Ao meus amigos Alexander Strelow, e famlia, Fernando, Felipe(Pride), e famlia, e Lo, obrigadopelos momentos de descontrao com o RPG, comidas gostosas e boa cerveja.

Alexandre Patriota e sua namorada, Fernanda, merecem um agradecimento adicional por teremme hospedado por uma semana, da ocasio da defesa da minha tese.

Desculpem aos meus amigos que esqueci de mencionar, mas agradeo a vocs por serem meusamigos.

i

Resumo

CAMPOS, T. F. Aplicaes de Approximate Bayesian Computation a controle de quali-dade. 2015. 120 f. Tese (Doutorado) - Instituto de Matemtica e Estatstica, Universidade de SoPaulo, So Paulo, 2015.Neste trabalho apresentaremos dois problemas do contexto de controle estatstico da qualidade:monitoramento on-line de qualidade e enviromental stress screening, analisados pela ptica baye-siana. Apresentaremos os problemas dos modelos bayesianos relativos a sua aplicao e, os reana-lisamos com o auxlio do ABC o que nos fornece resultados de uma maneira mais rpida, e assimpossibilita anlises diferenciadas e a previso novas observaes.

Palavras-chave: Approximate Bayesian Computarion, Controle de qualidade, Preditivismo.

iii

Abstract

CAMPOS, T. F. Applications of Approximate Bayesian Computation in quality control. 2015. 120 f. Tese (Doutorado) - Instituto de Matemtica e Estatstica, Universidade de So Paulo,So Paulo, 2015 .In this work we will present two problems in the context of statistical quality control: on line qua-lity monitoring and environmental stress screening, analyzed from the Bayesian perspective. Wewill present problems of the Bayesian models related to their application, and also we reanalyzethe problems with the assistance of ABC methods which provides results in a faster way, and soenabling differentiated analyzes and new observations forecast.

Keywords: Approximate Bayesian Computation, Quality control, Forecast.

v

Sumrio

Lista de Figuras ix

Lista de Tabelas xi

1 Introduo 11.1 Organizao do trabalho . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

2 Approximate Bayesian Computation 52.1 Introduo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.2 Funo de verossimilhana . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.3 Algoritmo ABC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2.3.1 Gnese do Approximate Bayesian Computation . . . . . . . . . . . . . . 72.3.2 Observando variveis aleatrias contnuas . . . . . . . . . . . . . . . . . . . . 82.3.3 MCMC-ABC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.3.4 Noisy ABC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.4 Toy model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.4.1 Beta-binomial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.4.2 Normal-gama . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.5 Comparando o ABC com um algoritmo MCMC . . . . . . . . . . . . . . . . . . . . 14

3 Gerando observaes de uma distribuio de Holgate 213.1 Introduo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 213.2 Distribuio de Poisson bivariada . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 213.3 Abordagem Bayesiana . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 233.4 Gerando observaes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

3.4.1 Gibbs Sampler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 243.4.2 Metropolis-Hastings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 253.4.3 ABC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 263.4.4 Comparando mtodos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

3.5 Prevendo novas observaes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 303.6 Observaes sobre o ABC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

4 ABC aplicado ao modelo Bayesiano no-paramtrico para o monitoramento on-line de qualidade de Taguchi para atributos. 454.1 Introduo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 454.2 O modelo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

4.2.1 Suposies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

vii

viii SUMRIO

4.2.2 Modelo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 464.2.3 Vantagens e desvantagens . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

4.3 ABC aplicado ao modelo no paramtrico . . . . . . . . . . . . . . . . . . . . . . . . 484.3.1 Gerando observaes de um processo de Dirichlet . . . . . . . . . . . . . . . . 484.3.2 Gerando observaes do processo de controle de qualidade. . . . . . . . . . . . 524.3.3 ABC configurado para gerar uma funo de probabilidade acumulada . . . . 524.3.4 Escolha do mtodo para gerar uma observao de um processo de Dirichlet . 54

4.4 Calculando a estimativa da proporo dos itens no-conforme ao longo do processode produo usando o amostrador de Gibbs . . . . . . . . . . . . . . . . . . . . . . . 54

4.5 Resultados das simulaes do ABC . . . . . . . . . . . . . . . . . . . . . . . . . . . . 564.5.1 Anlise dos hipermetros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

4.6 Inferncia preditiva . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 574.7 Observaes e consideraes finais . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

5 ABC aplicado a Enviromental Stress Screening 675.1 Introduo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 675.2 Abordagem Bayesiana para a metodologia ESS . . . . . . . . . . . . . . . . . . . . . 67

5.2.1 Funo de risco a priori . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 685.2.2 Anlise Bayesiana . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

5.3 Usando o ABC na metodologia ESS . . . . . . . . . . . . . . . . . . . . . . . . . . . 705.3.1 Exemplo numrico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

5.4 MCMC aplicado metodologia ESS . . . . . . . . . . . . . . . . . . . . . . . . . . 725.5 Taxa de fala crescente . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

5.5.1 Modelo com taxa de falha crescente . . . . . . . . . . . . . . . . . . . . . . . 745.5.2 Exemplo numrico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76

5.6 Previses de futuros enviromental stress screening . . . . . . . . . . . . . . . . . . . . 775.7 Consideraes finais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

6 Concluses 876.1 Sugestes para Pesquisas Futuras . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88

Referncias Bibliogrficas 89

Lista de Figuras

2.1 (a) Histograma da amostra mdia gerada no primeiro cenrio comparada com adensidade a posteriori (equao 2.12); (b) Histograma da amostra mdia geradano segundo cenrio comparada com a densidade a posteriori (equao 2.12); (c)Histograma da amostra mdia gerada no terceiro cenrio comparada com a densidadea posteriori (equao 2.12); (d) Histograma da amostra mdia gerada no quartocenrio comparada com a densidade a posteriori (equao 2.12) . . . . . . . . . . . . 16

2.2 (a) Funo de probabilidade acumulada de (, )|x; (b) Funo emprica mdia daamostra do primeiro cenrio; (c) Funo emprica mdia da amostra do segundo cenrio. 17

2.3 (a) Grfico de disperso de todas as amostras geradas no primeiro cenrio; (b) Grficode disperso de todas as amostras geradas no segundo cenrio. . . . . . . . . . . . . 18

2.4 (a) Funo de probabilidade acumulada de (, )|x; (b) Funo acumulada emp-rica mdia da amostra do primeiro cenrio; (c) Funo acumulada emprica mdiada amostra do segundo cenrio; (d) Funo acumulada emprica mdia da amostragerada pelo amostrador de Gibbs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

2.5 (a) Grfico de disperso de todas as amostras geradas no primeiro cenrio; (b) Gr-fico de disperso de todas as amostras geradas no segundo cenrio; (c) Grfico dedisperso de todas as amostras geradas pelo amostrador de Gibbs. . . . . . . . . . . . 20

3.1 Boxplot do nmero da iteraes que cada algoritmo necessitou para gerar a amostra. 333.2 Histogramas de observaes a posteriori de 1, 2 e 3 geradas pelo MCMC via

Random Walk . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 343.3 Histogramas de observaes a posteriori de 1, 2 e 3 geradas pelo MCMC via

Independence Chain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 353.4 Histogramas de observaes, geradas pelo ABC a posteriori de 1, 2 e 3 . . . . . . 363.5 Histograma em trs dimenses das amostras de 1 versus 2 geradas pelo Indepen-

dence Chain MCMC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 373.6 Histograma em trs dimenses das amostras de 1 versus 3 geradas pelo Indepen-

dence Chain MCMC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 383.7 Histograma em trs dimenses das amostras de 3 versus 3 geradas pelo Indepen-

dence Chain MCMC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 393.8 Histograma em trs dimenses das amostras de 1 versus 2 geradas pelo ABC . . 403.9 Histograma em trs dimenses das amostras de 1 versus 3 geradas pelo ABC . . 413.10 Histograma em trs dimenses das amostras de 3 versus 3 geradas pelo ABC . . . 423.11 Disgrama de disperso em trs dimenses da distribuio preditiva . . . . . . . . . . 433.12 Representao da distribuio preditiva na forma de um histograma em trs dimenses 44

ix

x LISTA DE FIGURAS

4.1 Nmero de linhas da matriz Ck . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 494.2 Nmero de linhas da matriz Ck . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 504.3 Boxplot do numero de interaes necessrias para gerar amostras utilizando o stick-

breacker e o Rapid Simulation, como amostradores da distribuio priori . . . . . . . 554.4 Estimativas de F (t)| = k geradas pelo ABC, usando o Stickbreaker e o Rapid

Simulations Of A Dirichlet Process . . . . . . . . . . . . . . . . . . . . . . . . . . . . 594.5 Estimativas de F (t)| = 1, obtidos pela equao 4.4, pelo ABC (algoritmo 14) e

Gibbs(algorimo 16) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 604.6 Estimativas de F (t)| = k para a = b = = = = 5, dado = k para k = 1, . . . , 10 614.7 Estimativas de F (t)| = 1 variando de acordo com valores de a . . . . . . . . . . . . 624.8 Estimativas de F (t)| = 1 variando de acordo com valores de b . . . . . . . . . . . . 634.9 Estimativas de F (t)| = 1 variando de acordo com valores de . . . . . . . . . . . . 644.10 Estimativas de F (t)| = 1 variando de acordo com valores de . . . . . . . . . . . . 654.11 Estimativas de F (t)| = 1 variando de acordo com valores de . . . . . . . . . . . . 66

5.1 Histograma das amostras da distribuio de probabilidade a posteriori de p, e e ggeradas pelo algoritmo 18 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

5.2 Funo de risco a posteriori . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 805.3 Histograma das amostras de (, p, g, e) gerados pelo MCMC, algoritmo 19 . . . . . 815.4 Funo de risco a priori quando a taxa de falha crescente . . . . . . . . . . . . . . . 825.5 Histograma das amostras da distribuio a posteriori de (, p, g, e, kg, ke) . . . . . 835.6 Funo de risco a posteriori . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 845.7 Grfico de barras das distribuies preditivas marginais de x, y, z e m suponto taxa

de falha constante. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 855.8 Grfico de barras das distribuies preditivas marginais de x, y, z e m suponto taxa

de falha crescente. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86

Lista de Tabelas

2.1 Sumrio do nmero de interaes necessrias para gerar uma observao usando oamostrador de Gibbs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

3.1 Sumrio do nmero de iteraes necessrias para gerar uma observao da distribui-o de Holgate usando os algoritmos 9, 10 e 11 . . . . . . . . . . . . . . . . . . . . . 29

3.2 Distribuio preditiva . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

4.1 Estatsticas descritivas do nmero de interaes que o ABC levou para gerar asestimativas de F (t)| = k, usando tanto o Stickbreaker como o Rapid Simulation . . 54

4.2 Distribuio preditiva do nmero de inspees necessrias para detectar um itemno-conforme, num prximo ciclo de produo . . . . . . . . . . . . . . . . . . . . . . 58

xi

xii LISTA DE TABELAS

Lista de Algoritmos

1 ABC -Rejeio 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 ABC -Rejeio 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 MCMC-ABC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 ABC - Rejeio 2 (Beta-Binomial) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 ABC - Rejeio 2 (Normal-gama) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 Amostrador de Gibbs (Normal-gama) . . . . . . . . . . . . . . . . . . . . . . . 147 Amostrador de Gibbs para gerar observaes da distribuio Holgate . . . . . . . . . 258 Metropolis-Hastings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 259 Random Walk Chain Metropolis-Hastings . . . . . . . . . . . . . . . . . . . . . . . . 2710 Independence Chain Metropolis-Hastings . . . . . . . . . . . . . . . . . . . . . . . . . 2711 ABC - Rejeio 2 (Holgate) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2812 Gerador de processo de Dirichlet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5213 Gerador de processo de controle de qualidade de Taguchi . . . . . . . . . . . . . . . . 5214 Processo de estimao de F (t)| = k . . . . . . . . . . . . . . . . . . . . . . . . . . . 5315 Amostrador de Gibbs para obter o estimador bayesiano no-paramtrico da funo

de sobrevivncia para o caso de censura intervalar quando a priori processo deDirichlet. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

16 amostrador de Gibbs para obter observaes da distribuio 4.4 . . . . . . . . . . . . 5617 ABC para gerar amostras de (, p, g, e|x, y, z,m) . . . . . . . . . . . . . . . . . . 7018 ABC para gerar amostras de (, p, g, e|x, y, z,m) . . . . . . . . . . . . . . . . . . 7219 MCMC via independence chain para gerar amostras do modelo 5.8 . . . . . . . . . . 7320 ABC para gerar amostras de (, p, g, e, kg, ke|x, y, z,m) . . . . . . . . . . . . . . 76

xiii

xiv LISTA DE ALGORITMOS

Captulo 1

Introduo

Os fenmenos da natureza sempre chamaram ateno do homem. Na pr-historia o conhecimentosobre os fenmenos naturais representava a diferena entre viver ou morrer. A capacidade de poderexplicar, e com isso poder prever tais fenmenos se tornou importante para a evoluo humana.Tal cenrio demonstra o importante papel que um esforo de explicao e previso teve em termosprticos tanto quanto em termos explicativos

Com a necessidade de obter novas explicaes sobre o que era contemplado pelo ser humano,sugiram novas formas de responder aos questionamentos levantados perante tais fenmenos. Opensamento matemtico tornou-se uma alternativa para explicar tais contemplaes. Por considerara aleatoriedade, caos dos fenmenos observados, a abordagem estatstica a forma mais adaptadatal tarefa.

Modelos estatsticos vem evoluindo de maneira similar s explicaes dos fenmenos observados.Essa evoluo do modelo deve-se, em grande parte, a evoluo tecnolgica. Atualmente, modeloscom suposies cada vez mais sofisticadas podem ser criados, avaliados e usados no mundo prtico.Dos modelos, mais especificamente modelos estatsticos, possuidores de sofisticaes tericas, queaumentam o nvel de complexidade do mesmo, os que mais se beneficiaram do avano tecnolgicoso os modelos sob a ptica Bayesiana. Um exemplo so modelos que no tinham um teorema queestipulava uma famlia conjugada e, com o avano tecnolgico, passaram a ser frequentes.

Com a constante sofisticao dos modelos, desejo latente de alguns estatsticos, um fenmenoaflorou, a extrema complexidade computacional dos mesmos. At que ponto a contnua sofisticaode modelos estatsticos possvel, levando em considerao a complexidade gerada por ela? Quaiscustos a utilizao de um modelo super sofisticado pode trazer ao seu usurio? Quais os ganhos eperdas tem-se por utilizar um modelo menos sofisticado, mas mais acessvel?

O controle estatstico do processo de produo de um determinado item uma rea de estudoem que a complexidade de modelos estatstica uma caracterstica indesejvel. Quanto maior acomplexidade do modelo utilizado maior o tempo necessrio para diagnosticar uma possvel quebrana qualidade do processo de produo.

Assegurar o controle de qualidade de uma linha de produo essencial para evitar a comerci-alizao de produtos defeituosos. Tsunemi et al. (2012) e Barlow et al. (1994) propem modelosBayesianos para dois procedimentos de monitoramento da qualidade do processo de produo. Oprimeiro trata o monitoramento de um processo de produo em linha segundo a abordagem deTaguchi (ver Taguchi et al. (1989)). Tsunemi et al. (2012) apresentam uma abordagem Bayesianano-paramtrica para modelar a verificao de um item no-conforme durante uma inspeo dociclo de produo.

J em Barlow et al. (1994) o contexto abordado o da ps-produo. assumido que ositens produzidos podem atender ou no a qualidade estipulada. Um lote destes itens submetidoa um ensaio de estresse simulando, de forma acelerada, o tempo de uso. O trabalho em questoapresenta uma abordagem bayesiana para estimar a proporo de itens de baixa qualidade queforam produzidos, bem como o tempo mdio da vida til dos itens de boa e baixa qualidade.

Em ambos os modelos h um fato em comum, uma complexidade computacional intensa, aoponto de inviabilizar o seu uso. Tal complexidade devida a suposies tericas extremamentesofisticadas, feitas nos seus respectivos trabalhos.

1

2 INTRODUO 1.1

Quando um modelo em estudo complicado de analisar, pelo fato de no ter uma forma analticaamigvel, a opo para a utilizao do mesmo, a gerao de observaes deste modelo em estudo.

Existem diversos mtodos para gerar observaes de um modelo estatstico, todos baseados nadistribuio de probabilidade associada ao modelo. O mais simples utilizar a funo inversa dadistribuio acumulada (F1).

Quando a complexidade envolvida na anlise grande, so necessrios mtodos computacionaismais potentes. Para os modelos Bayesianos mais sofisticados os mtodos mais usuais so baseadosem Markov Chain Monte Carlo, MCMC. O amostrador de Gibbs e o algoritmo de Metropolis-Hastings so os algoritmos mais utilizados.

O lado negativo dessa estratgia de anlise o fato de que esses algoritmos dependem direta,ou indiretamente, da forma da funo de verossimilhana. Dependendo da complexidade da fun-o de verossimilhana, como o caso do modelo apresentado por Tsunemi et al. (2012), gerarobservaes do modelo continua sendo uma tarefa extremamente difcil. Os modelos apresentadosem Tsunemi et al. (2012) e Barlow et al. (1994) so dois modelos que apresentam uma comple-xidade computacional to grande que impossibilita at mesmo a utilizao de mtodos MCMC,como ser visto e discutido mais a frente nesse trabalho.

Em casos como esse que o Approximate Bayesian Computation (ABC) se destaca. Em suaessncia est o fato de no utilizar, de forma numrica, a funo de verossimilhana. Essa premissa o grande atrativo do ABC e, tambm, um ponto de grande discusso. Pela sua sua ideia motriz,o ABC o candidato ideal para gerar observaes dos modelos supracitados.

O ABC bastante usado em reas biolgicas: Beaumont et al. (2002) implementam o mtodoem gentica de populaes e discutem o problema de inmeros parmetros indesejados nos modelosestatsticos utilizados nessa rea de pesquisa. Em Cornuet et al. (2008) proposto um programacomputacional, Do It Yourself Approximate Bayesian Computation (DIY ABC ) para analisar dadosde populao gentica e aferir informaes sobre a histria evolutiva da mesma. Beaumont (2010)descreve o ABC, bem como tcnicas de MCMC nos contextos de populaes genticas, demografiade colonizao, histria demogrfica humana, filogeografia, entre outros.

No presente trabalho introduzido o uso do ABC no contexto de controle de qualidade. Comisso esperado que as anlises dos modelos de Tsunemi et al. (2012) e Barlow et al. (1994) fiquemmais palpveis e, assim, fazer com eles sejam utilizados com mais frequncia . Outra ambio apossibilidade estender os resultados, j conhecidos, para outros casos igualmente interessantes. Porexemplo, relaxar a suposio de taxa de falha constante, presente em Barlow et al. (1994), usandouma taxa de falha crescente.

Com o objetivo de manter a complexidade computacional no nvel mais baixo possvel, decidiu-sepor apenas implementar a verso mais simples do ABC.

1.1 Organizao do trabalhoEsse trabalho se organiza da seguinte forma. O captulo 2 apresenta o algoritmo Approximate

Bayesian Computation em seus detalhes. Discute-se o fato do ABC usar, ou no, a funode verossimilhana de um modelo. O funcionamento ilustrado atravs de toy models. Suas pe-culiaridades so discutidas, o que muda nas amostras geradas quando mudamos os argumentoscomputacionais do ABC e, por fim, o desempenho do ABC comparado com o do MCMC .

No captulo 3 apresentado um breve resumo sobre a distribuio de Holgate, que possui umafuno de probabilidade, e por consequncia uma funo de verossimilhana, analiticamente com-plicada. O ABC, e o MCMC , so configurados para gerar observaes de uma distribuio deprobabilidade a posteriori, distribuio dos dados observados a distribuio Holgate, e seus de-sempenhos comparados. Acredita-se que esse pode ser um exemplo esclarecedor de que o ABCpode gerar amostras de um modelo sem necessitar de sua funo de verossimilhana.

Nos captulos 4 e 5 seguido o seguinte roteiro, para ambos os captulos. Apresentar o pro-blema em estudo, revisando os artigos de origem do problema, justificar o uso do ABC, aplicar omesmo, compara-lo com o MCMC e, por fim usar os dados gerados para obter a distribuio deprobabilidade preditiva de um nova observao do experimento de interesse.

1.1 ORGANIZAO DO TRABALHO 3

Finalmente, no captulo 6, conclui-se o trabalho, e so apresentadas diretrizes para futurostrabalhos.

4 INTRODUO 1.1

Captulo 2

Approximate Bayesian Computation

2.1 IntroduoNa inferncia Bayesiana h uma necessidade recorrente de gerar observaes de dados de inte-

resse, e os motivos para tal necessidade so variados. Em grande parte das anlises bayesianas essanecessidade oriunda ao fato da distribuio de probabilidade a posteriori no possuir uma formaanaltica fechada, o que bastante frequente quando no possvel trabalhar dentro da famlias dedistribuies conjugadas usuais, que so modelos mais simples e amplamente conhecidos.

Quando no se tem uma forma analtica fechada, a obteno de resumos (como esperana,varincia, mximo, mnimo, entre outros) a posteriori fica bastante dificultada. A gerao de obser-vaes de uma distribuio de probabilidade de grande ajuda neste caso, j que existem mtodosnumricos que possibilitam calcular probabilidades, estimar densidades, fazer qualquer manipulao(inferncia) que se deseje, e para isso basta ter um nmero de observaes suficientemente grande.

Alguns dos algoritmos que geram dados observveis (tambm chamados de amostradores) maisutilizados so baseados em MCMC , sigla para Markov Chain Monte Carlo, um algoritmo quepropem gerar observaes de uma deseja da distribuio de probabilidade atravs da evoluode cadeias de Markov. Pode-se dizer que o amostrador de Gibbs (Casella e George, 1992) umadas variantes do MCMC mais utilizadas em anlises devido a facilidade de implement-lo. Outravariante basante utilizada o algoritmo de Metropolis-Hasting (Gamerman e Lopes, 2006), maisrobusto que o amostrador de Gibbs. Em esmagadora maioria os mtodos de simulao necessitam,de alguma forma, da funo de verossimilhana, geralmente atravs de manipulaes algbricas.

O ABC, tambm conhecido como likelihood-free method, um amostrador que gera dadosobservveis sem necessitar calcular, ou efetuar clculos, em uma funo de verossimilhana. Suagnese, funcionamento e peculiaridades sero apresentado a seguir, bem como uma breve discussosobre o que verossimilhana, e se o ABC a utiliza ou no.

2.2 Funo de verossimilhana usualmente dito que a funo de verossimilhana nada mais que a funo de probabilidade

(funo densidade de probabilidade) vista em funo do parmetro, usualmente retratado por ,com uma observao amostral x, podendo ser uma nica observao ou uma amostra com vriasobservaes, mas fixa, j observada (DeGroot, 1989, captulo 6). Esse fato faz com que a funo doparmetro no seja mais uma funo de probabilidade, mas sim uma funo, de qualquer, exceto deprobabilidade ou densidade de probabilidade. Em Fisher (1922, pag. 310) dito The Likelihoodthat any parameter (or set of parameters) should have any assigned value (or set of values) isproportional that if this were so, the totality of observations should be that observed. Portanto afuno pode ser escrita como:

L() f(x|). (2.1)

A verossimilhana de algum parmetro (), dado a observao de uma amostra (x) de tamanho1(um), ou seja, uma nica observao oriunda de uma populao (ou varivel) discreta, pode serescrita da forma:

L(|x) = P (x|). (2.2)

Considerando agora uma amostra de tamanho n, x = (x1, x2, . . . , xn) (supondo observaes inde-

5

6 APPROXIMATE BAYESIAN COMPUTATION 2.2

pendentes e identicamente distribudas), pode-se escrever:

L(|x) =ni=1

P (xi|). (2.3)

Analogamente, quando observada uma amostra, de tamanho um e de tamanho n, respectiva-mente, de uma populao (ou varivel) continua, tem-se:

L(|x) = f(x|);

L(|x) =ni=1

f(xi|).

Em Bayarri e DeGroot (1992), alm da definio supracitada, dito que a funo de verossi-milhana, usualmente, tambm pode ser considerada como sendo proporcional a funo de proba-bilidade, ou funo de densidade de probabilidade, de quantidade observadas em um determinadoexperimento, condicionadas a quantidades no observadas neste mesmo experimento,

L() f(o que foi observado|o que no foi observado) (2.4)

De maneira mais simples, esta forma de pensar abre um precedente para interpretar a funo deverossimilhana como sendo o condicionamento da amostra observada nas amostras que poderiamser observadas, mas no foram. Assim, no h um interesse direto no parmetro , caso haja umno problema em estudo.

Note que a nica diferena entre as duas interpretaes exatamente o que considerado,ou no, parmetro. Na primeira abordagem, equao 2.1, o parmetro, geralmente, consideradouma quantidade fixa, mas desconhecida. J na segunda abordagem, equao 2.4, isso pode nonecessariamente acontecer. Esse fato fica mais claro com o exemplo contido em Bayarri e DeGroot(1992), exemplo 2.2.1.

Exemplo 2.2.1. Considere um problema em que, para cada valor do parmetro , uma varivelaleatria Y tem funo densidade de probabilidade f(y|). Considera-se que Y no pode ser, por sis, observada, mas sim pode-se observar uma pertubao aleatria X de Y . Isto , assume-se queo que observado uma realizao da varivel aleatria X com funo densidade de probabilidadecondicional f(x|y, ) = f(x|y).

Qual a funo de verossimilhana para esse problema? Existem diversas possibilidades. Pri-meiro deve-se decidir onde o valor desconhecido de y de Y deve ser includa na funo de veros-similhana. Pode-se argumentar que, j que y no pode ser, nem nunca ser, observado, maisapropriado que no seja includa ao se definir a funo de verossimilhana, e assim definir:

L0() f(x|). (2.5)

Por outro lado, j que a formulao bsica do problema em funo de y, ento pode-se pensarque y deve ser levada em considerao na funo de verossimilhana (este caso seria bvio sehouvesse interesse tambm em y). Neste caso precisa-se decidir em como considerar y, como umarealizao de uma varivel aleatria, e portanto ficando a esquerda de |, da seguinte forma,

L1() f(x, y|), (2.6)

ou considerar y como uma quantidade no observada, e desconhecida, e assim fazendo-o ficar adireita de |, resultando uma funo de verossimilhana da seguinte forma:

L2() f(x|y, ) = f(x|y). (2.7)

Qual funo de verossimilhana usar no , a principio, fcil de escolher, e necessrio realizar

2.3 ALGORITMO ABC 7

um julgamento de carter subjetivo para poder decidir qual usar.

A motivao desse exemplo mostrar que a funo de verossimilhana pode ser ambgua e dif-cil de ser definida mesmo em problemas com um nvel baixo de complexidade. Assim, a afirmaofeita por Bayarri e DeGroot (1992) (...there is not such a thing as a likelihood function that canbe unambigouosly defined in all statistical problems.) ganha bastante fora. No h uma definionica da funo de verossimilhana para todos os problemas estatsticos. Com isso em mente umquestionamento pode ser levantado: um conceito, to difundido na estatstica, pode ser delimitadosomente pela definio de uma nica funo, que se mostrou ambgua? Pode-se considerar que omodelo probabilstico, do experimento de interesse, um sinnimo da verossimilhana, no sentidode que ambos podem ser facilmente confundidos? O exemplo 2.2.1 deixa claro que funo de veros-similhana e modelo probabilstico, funo (densidade) de probabilidade, so distintas, apesar deuma depender da outra.

2.3 Algoritmo ABCEm qualquer experimento sempre so feitas suposies sobre os dados observveis. Essas supo-

sies sempre so projetadas na funo de verossimilhana, e algumas dessas suposies tornam afuno de verossimilhana intratveis, deixando qualquer clculo impossvel, se no computacional-mente bastante custoso. Nos casos em que a intratabilidade da funo de verossimilhana no toextrema, algoritmos MCMC so bastante utilizados, mas o que pode ser feito quando a intrata-bilidade to grande que um algoritmo MCMC demoraria um tempo maior do que o disponvelpara gerar uma amostra desejvel? Nesta seo est apresentado uma alternativa para casos comoesse, o algoritmo ABC (Approximate Bayesian Computation), que no utiliza a funo deverossimilhana. Por esse motivo, esse mtodo chamado de likelihood-free.

2.3.1 Gnese do Approximate Bayesian Computation

A primeira apario do ABC, ainda como um conceito apareceu em Rubin (1984). nesse tra-balho dito que: Suppose the model for data set X is given by f(X|)p() where is the parameterwhose posterior distribution is to be calculated. Of course, the posterior distribuition of given ob-served X is calculated via Bayesians theorem, but how is conceptual content of this theorem easilyconveyed? Consider the simple description.

Suppose we first draw equally likely values of from p(), and label this 1, . . . , s. The j , j =1, . . . , s can be thought of as representing the possible populations that might have generated theobserved X. For each j, we now draw an X, that might have been observerd under the full modelf(X|)p(). Now some of the X will look just like the observerd X and many will not; of course,subject to the degree of rounding and the number of possible values of X, s might have to be verylarge in order to find generated Xj that agree with observerd X, but this creates no problem for ourconceptual experiment. Suppose we collect together all Xj match the observerd X, and then all jthat correspond to these Xj. This collection of j represents the values of that could have generatedthe observerd X; formally, this collection of values represents the posterior distribution of .

De maneira resumida, a ideia motriz do ABC a suposio de que, um valor oriundo dadistribuio de |X = x capaz de gerar, simular, um conjunto de dados igual , ou suficientementeprximos como mostrado a seguir, do conjunto de dados originalmente observado.

Essa ideia simples, desenvolvida por Rubin (1984), acabou se tornando uma alternativa com-putacional para o processo inferencial. Os casos que apresentam uma funo de verossimilhanaintratvel, beneficiam-se bastante do algoritmo ABC. Isso pode ser visto mais detalhadamente noscaptulos 4 e 5.

A transcrio, para pseudo-cdigo, da ideia presente em Rubin (1984) est presente no algo-ritmo 1. Este algortimo tambm pode ser encontrado em outros trabalhos, que contam a historia doABC, como por exemplo Marin et al. (2011) que, alm disso apresentam outras verses do ABC.

Para melhor entendimento do algoritmo 1, suponha que se deseja gerar amostras da distribui-o a posteriori de uma quantidade no observvel de um determinado experimento. Seja essaquantidade desconhecida e () sua distribuio de probabilidade a priori. Seja X uma varivelobservvel e x o valor observado de X, podendo ser uma nica ou vmrias observaes. Suponha


um modelo probabilstico para a varivel X dado o valor da quantidade desconhecida . Denota-sepor f(x|) tal modelo probabilstico. A verossimilhana gerada pela observao do experimento, jdefinida anteriormente, L(|x). Assim a distribuio de probabilidade a posteriori de |X = x dada por (|X = x) L(|x)().

Considerando um caso em que se tem uma funo de verossimilhana intratvel, mas em quese sabe como gerar uma observao de f(|), um exemplo disso est bem detalhada no captulo 3.O algoritmo 1, gera um da distribuio priori, denominado (passo 3). Ento, com geradouma amostra y oriunda de f(|) (passo 4). Se a amostra simulada igual a amostra observada,ento uma observao de (|X = x) (Passo 5).

Perceba que o ABC necessita, de maneira vital, que haja uma maneira de replicar o experimentoobservado: se no se sabe como fazer isso, ento no h como implementar qualquer verso do ABC.

Note que em nenhum momento foi feita nenhuma manipulao, algbrica ou numrica, na funode verossimilhana. O processo inferencial, a maneira de como a informao de obtida atravsda amostra observada, x, baseia-se somente no ato de replicar o experimento observado e anotarseus resultados. Por esse motivo o ABC tambm chamado de likelihood-free method.

Algoritmo 1 ABC -Rejeio 11: for i = 1 to N do2: repeat3: Gerar de ()4: Gerar y de f(|)5: until x=y6: i =

7: end for

A prova matemtica de que uma amostra gerada pelo ABC pode ser considerada como umaamostra da distribuio a posteriori de interesse dada a seguir. Seja f(i) a densidade de 1 geradapelo algoritmo 1. Assim,

f(i) yX

(i)f(y|i)Ix(y) = (i)f(x|i)

(i|x),

onde X o conjunto, enumervel, de valores possveis de x.Uma demonstrao deveras simples, mas o suficiente para dar credibilidade matemtica ao

ABC.

2.3.2 Observando variveis aleatrias contnuas

Do ponto de vista computacional, a condio de parada do algoritmo 1, passo 5, inalcanvel,quando os valores x e y pertencem ao conjunto dos nmeros reais.

Sofrendo desse problema, observando valores contnuos, Pritchard et al. (1999) utilizam oABC de forma modificada. Com o intuito de relaxar a condio de igualdade do algoritmo 1. Issofoi possvel utilizando mecanismos,que hoje so chamados de tunning parameters.

Devido ao fato de observar uma amostra de tamanho grande, Pritchard et al. sumarizam aamostra, condensando asim a informao contida nela e assume que a condio de igualdade en-tre as amostras, geradas e observadas, pode ser substituda por: a amostra gerada, devidamentesumarizada, est suficientemente perto da amostra observada, tambm devidamente sumarizada.

Com esse relaxamento, a ideia original de Rubin (1984) se modifica. Agora, se um valor

gerado de () pode gerar uma amostra y na vizinhana da amostra observada x, ento esse valor especfico considerado como uma realizao da distribuio a posteriori (|x).

A transcrio, para pseudo-cdigo, do algoritmo apresentado em Pritchard et al. (1999) encontra-se logo a baixo, algoritmo 2. Nele so presentes os seguintes tunning parameters:

2.3 ALGORITMO ABC 9

() a funo sumrio, que utilizada para resumir a informao da amostra em uma ordem degrandeza menor que a da amostra(e.g. estatsticas suficientes);

() > 0 a funo distncia, que vai informar quo longe, ou perto, as duas amostras, observadae a simulada, esto;

> 0 a distncia mxima tolerada entre as duas amostras.

Algoritmo 2 ABC -Rejeio 21: for i = 1 to N do2: repeat3: Gerar de ()4: Gerar y de f(|)5: until ((x) (y)) 6: i =

7: end for

A demonstrao terica de que a amostra gerada pelo algoritmo 2 segue, aproximadamente, adistribuio de probabilidade a posteriori de interesse segue abaixo.

Considere o conjunto A,x = {y X |((x) (y)) }, todos os valores y que satisfizeram acondio de aceitao do algoritmo 2, para um determinado. Note que no contexto que envolve oalgoritmo 2, quando observamos variveis aleatrias continuas, X no enumervel.

A distribuio a posteriori conjunta de (i, y) da forma:

(, y|x) =()f(y|)IA,x(y)A,x ()f(y|)d

.

Assim,

(|x) =(, y|x)dy (|x). (2.8)

A qualidade da aproximao, que pode no estar muito clara na demonstrao acima, dependeda escolha conjunta de e . A escolha de um sumrio de amostra representativa o suficientealiado com uma tolerncia pequena o suficiente, deve produzir uma aproximao, suficientementeboa.

Quando igual a zero, mesmo caso do ABC-Rejeio quando se observam variveis discre-tas, tem-se a amostra exata da distribuio alvo, mas, quando so observadas variveis contnuas,escolher um igual a zero revela-se uma impossibilidade computacional. A escolha do valor de uma negociao entre os custos computacionais e a qualidade da aproximao. Uma discusso maisaprofundada sobre as diferenas entre (|x) e (|x), atravs de resultados empricos, pode serencontrada em Sisson et al. (2007), e resultados tericos sobre limites superiores, em casos bemespecficos, para podem ser vistos em Dean et al. (2014).

Quando se tem um conjunto de dados com vrias dimenses, sumrios de amostra de dimen-ses menores, condensando a maior quantidade de informao possvel da maneira mais simplesque se possa, so essenciais. Mas nem sempre possvel condensar a informao dos dados damaneira que mais til. A m escolha de funes sumrios, alm de aumentar a diferena entre(|x) e (|x), podem levar a uma perda de informao o que implica numa inferncia de mqualidade (Csillry et al., 2010). Uma reviso pode mtodos para escolher pode ser encontradaem Blum et al. (2013).

importante frisar que em cada caso, a escolha de e deve ser feita depois de investigaes daqualidade da amostra gerada, em contra partida com custos computacionais envolvidos na execuodo algoritmo.


2.3.3 MCMC-ABC

Outras variaes do algoritmo ABC, desenvolvidas para superar dificuldades encontradas emdiversos problemas, podem ser encontradas em Marin et al. (2011). Uma dessas variaes oMCMC-ABC que mescla o ABC com a classe de algoritmos MCMC . A motivao para taljuno foi a tendncia em gerar pontos com baixa probabilidade a posteriori que o ABC-Rejeiotem em alguns casos. O MCMC -ABC fez sua primeira apario em Marjoram et al. (2003),com a ideia de adicionar uma nova condio de aceitao para , condio essa tambm presenteno algoritmo de Metropolis-Hastings, alm da adio de um kernel de transio para gerar novoss dependendo do ultimo aceito. A transcrio do MCMC-ABC pode ser vista a seguir,algoritmo 3.

Algoritmo 3 MCMC-ABC

1: Use o algoritmo 2 para obter uma realizao ((0), y(0))2: for t = 1 to N do3: Gerar do kernel q(|(t1))4: Gerar y de f(|)5: Gerar u de Uniforme(0, 1)6: if u (

)q((t1)|))((t1))q(|(t1)) e ((x) (y)) then

7: ((t), y(t)) = (, y)8: else9: ((t), y(t)) = ((t1), y(t1))

10: end if11: end for

Ao custo de aumentar a complexidade computacional do ABC-Rejeio, que pode ser imprati-cvel em algumas situaes, o MCMCABC promete diminuir, talvez at eliminar, a tendnciado ABC-Rejeio gerar pontos em regies de baixa probabilidade a posteriori.

2.3.4 Noisy ABC

Outra variao do ABC, tambm comentada em Marin et al. (2011), mas apresentada pelaprimeira vez em Wilkinson (2008) chama-se Noisy ABC. A particularidade desse mtodo aceitarque o ABC gera amostras perfeitas, mas da distribuio errada.

Semelhante a um modelo linear, ou modelo linear generalizado, onde os dados so observaesde um modelo adicionado de um erro aleatrio, Wilkinson (2008) interpreta a distribuio geradapelo ABC como se fosse o melhor modelo de predio possvel e assume que h uma discrepnciaentre esse melhor modelo de predio e os dados observados. Essa discrepncia pode representartanto um erro de medida nos dados como um erro de modelo. Fazendo esta discrepncia ser comoo erro aleatrio de um modelo linear, possvel escrever os dados observados da seguinte forma:

x = P (|) + ,

uma realizao do modelo com o melhor argumento , adicionado de um termo aleatrio . Por serum termo aleatrio assumido uma distribuio para , denominada por ().

A condio de aceitao do ABC, na metodologia presente em Wilkinson (2008), baseia-se nadistribuio () avaliada no ponto x y, o que torna a regra de aceitao varivel a cada iteraodo algoritmo. Detalhes em como construir () encontram-se em Wilkinson (2008).

Todos as simulaes desse captulo, bem como de todo o trabalho, feitas usando o ABC foramgeradas pelo ABC-Rejeio. Nenhuma outra variao foi utilizada, pois tendo em vista a altacomplexidade ds modelos discutidos neste trabalho, manter a simplicidade computacional do ABCtornou-se prioridade.

2.4 TOY MODEL 11

2.4 Toy modelPara ilustrar o funcionamento do ABC, ele foi implementado para gerar amostras de duas distri-

buies a posteriori de modelos Bayesianos j conhecidos, toy models. Em um dos modelos a varivelobservada discreta e no outro uma varivel contnua. Isso permite observar os algoritmos 1 e 2em ao. Os modelos em estudo so:

1. Considerando X uma sequncia de n ensaios de Bernoullis

Beta(a, b),X| Bernoulli(),

|X = x Beta(a+ni=1

xi, b+ 1ni=1

xi); (2.9)

2. Considerando uma amostra de tamanho n

Gama(0, 0),| Normal(0, ()1),

X|, Normal(, 1), |X = x Gama(n, n), (2.10)

|, X = x Normal(n, [(n+ ) ]1); (2.11)

onde

n = 0 + n/2,

n = 0 +ni=1(xix)2

2 +n(x0)2

2(n+) ,

n =nx+0n+ .

Como sabido, em outros mtodos computacionais que, tambm utilizam parmetros de tun-ning, Random Walk Chain MCMC por exemplo (Gamerman e Lopes, 2006, pag. 198), no existeuma regra geral para determinar o melhor valor para o parmetro de tunning. Por isso foram es-colhidos alguns valores para esses parmetros, gerando assim alguns cenrios diferentes em que oABC ser implementado.

Na escolha de (), tem-se a crena que usar um vetor cujas componentes so estatsticas su-ficientes, ou utilizar um vetor cujas as componentes so funes da amostra que esto parte deestatsticas suficientes, como o boxplot (um vetor cujas componentes so os trs quartis), que parte de uma estatstica suficiente(estatsticas de ordem), devem gera amostras de boa qualidade.Dependendo das caractersticas amostrais (tendncia central, disperso, assimetria, curtose, etc.)de x que deseja-se replicar na amostra gerada y ser necessrio uma funo () especifica.

Para escolher a funo distncia, (), a distncia euclideana parece natural para isso, mastambm foram testadas outras composies de distncias.

O nvel de tolerncia funciona como em qualquer outro algoritmo: quanto menor, mais prximoesto as duas amostras, e mais tempo computacional gasto para gerar uma amostra y que satisfaaa condio de aceitao desejada.

2.4.1 Beta-binomial

O primeiro dos toy models a ser estudado o modelo beta-binomial. O ABC configurado paraeste caso est apresentado em seguida, algoritmo 4.

Considere a realizao de 10 ensaios de Bernoulli, x = (1 0 1 0 0 0 0 0 0 1) e a distribuio a priori beta(1/2, 1/2). Pela equao 2.9 tem-se que,

|X = x Beta(7/2, 15/2). (2.12)


Algoritmo 4 ABC - Rejeio 2 (Beta-Binomial)1: for i = 1 to N do2: repeat3: Gerar de Beta(a, b)4: Gerar y de Bernoulli()5: until ((x) (y)) 6: i =

7: end for

O algoritmo ABC-Rejeio 2 (Beta-Binomial) foi configurado em quatro cenrios diferentes

1. (x) = x, ((x) (y)) = (x y)2 e = 102;

2. (x) = x, ((x) (y)) a distncia euclidiana e = 102

3. (x) = x, ((x) (y)) = (x y)2 e = 106;

4. (x) um vetor cujas componentes so os quantis de ordem 25%, 50%, 75%, ((x) (y)) a distncia euclidiana = 106.

Note que, quando o algoritmo 2 implementado com () descrito como no segundo cenrio,independente das escolhas de () e de , o algoritmo executado acaba sendo o algoritmo 1.

Para cada cenrio foram gerados 100 amostras de tamanho 1000 cada uma, e com essas amostrasobtiveram-se as seguintes informaes:

1. Para gerar uma amostra de tamanho 1000 foi necessrio gerar, em mdia, 7153 , em outraspalavras necessrio gerar 8 para aceitar 1;

2. Para gerar uma amostra de tamanho 1000 foi necessrio gerar, em mdia, 1840000 , emoutras palavras necessrio gerar 1841 para aceitar 1;

3. Para gerar uma amostra de tamanho 1000 foi necessrio gerar, em mdia, 153200 , em outraspalavras necessrio gerar 154 para aceitar 1;

4. Para gerar uma amostra de tamanho 1000 foi necessrio gerar, em mdia, 153300 , em outraspalavras necessrio gerar 154 para aceitar 1.

Ainda utilizando as rplicas, os grficos da figura 2.1 foram feitos com a inteno de saber quoaproximado est a amostra da posteriori gerada pelo ABC da densidade da distribuio a posterioridesejada (equao 2.12).

Em cada cenrio, cada uma das rplicas foi ordenada, gerando assim as estatsticas de ordem.Calculou-se a mdia de cada uma das estatstica de ordem nas 100 rplicas, gerando assim umaamostra composta pelas mdias de cada estatstica de ordem. Esta amostra composta pelas mdiasdas estatsticas de ordem ser denominada amostra mdia. Na figura 2.1, esto representados cadauma das amostras mdias para cada cenrio.

Percebe-se que a configurao do cenrio 1 no to boa quanto as dos outros trs cenrios, poisfoi a que mais se afastou da distribuio alvo, equao 2.9, embora tenha sido a que necessitou demenos interaes para a amostra. A amostra gerada com as configuraes do cenrio 2, teoricamente,tem a melhor qualidade, pelo fato de usar toda a amostra x como funo sumrio. Mas foi o casoem que o ABC mais foi dispendioso computacionalmente, j que necessitou mais de 10 vezesmais interaes que o ABC nos outros dois cenrios, para gerar a amostra. As configuraes doscenrios 3 e 4 geraram amostras muito prximas da densidade alvo, tanto quanto as do cenrio 2,e com uma taxa de aceitao bem maior, mostrando-se as melhores alternativas, dentre os quatrocenrios, aliando qualidade da amostra com consumo de recursos computacionais intermedirio.Como comentado anteriormente, a escolha de uma negociao entre recursos computacionais equalidade da amostra.

2.4 TOY MODEL 13

2.4.2 Normal-gama

O segundo toy model estudado neste captulo um modelo conjugado bastante conhecido, mo-delo normal-gama. O algoritmo 2 foi configurado para gerar amostras da distribuio de (, )|X =x apresentado em seguida.

Algoritmo 5 ABC - Rejeio 2 (Normal-gama)1: for i = 1 to N do2: repeat3: Gerar de Gama(0, 0)4: Gerar de Normal(0, ( )1)5: Gerar y de Normal(, 1)6: until ((x) (y)) 7: i =

8: i =

9: end for

Os hiperparmetros foram escolhidos da seguinte forma, = 1, 0 = 0.45, 0 = 1.53 e 0 = 0.81.Novamente foram escolhidos cenrios diferentes para estudar o desempenho do ABC. Estes estoapresentados a seguir:

1. (x) =

(ni=1

xi,

ni=1

x2i

), ((x) (y)) a distncia euclidiana e = 102;

2. (x) um vetor cujas as componentes so os quantis de ordem 25%, 50%, 75%, ((x)(y)) a distancia euclidiana e = 102.

Novamente, para cada cenrio, foram gerados 100 amostras de tamanho 1000, com uma mesmaamostra observada para cada um dos cenrios. A amostra observada nesse exemplo, x, tem tamanho

100 e possui

(ni=1

xi,ni=1

x2i

)= (40.46047, 45.42041) eQuantis(25%, 50%, 75%) = (0.01946887, 0.46127636, 0.83666056).

A distribuio a posteriori de (, )|x pode ser escrita da seguinte forma:

|X = x Gama(51.53, 15.33598) (2.13)|,X = x Normal(0.4050541, (101 )1). (2.14)

Note que (x) um vetor em qualquer um dos dois cenrios.Antes de apresentar os resultados das simulaes interessante ressaltar que no modelo beta-

binomial possui um parmetro escalar. J o modelo normal-gama tem um vetor de parmetros(, ). No h como usar estatsticas de ordem, como foram usadas no toy model com o modeloBeta-Binomial, para os vetores observados. A alternativa utilizada foi calcular a funo emprica deprobabilidade acumulada para cada uma das rplicas geradas, e ento calcular a mdia das funesempricas ponto a ponto. Assim h a possibilidade de comparar a funo emprica mdia de cadacenrio com funo de probabilidade acumulada da distribuio alvo.

Os seguintes resultados foram obtidos, para cada cenrio:

1. Para gerar uma amostra da distribuio a posteriori de tamanho 1000 foi necessrio gerar, emmdia, 14915067 . Significa que necessrio gerar 14916 para aceitar 1 desses 14916;

2. Para gerar uma amostra da distribuio a posteriori de tamanho 1000 foi necessrio gerar, emmdia, 298800.1 . Significa que necessrio gerar 299 para aceitar 1 desses 299.

Na figura 2.2 so apresentados as curvas de nveis das funes de probabilidade acumulada dadistribuio alvo, figura (a), e das funes empricas mdia de cada um dos cenrios, figuras (b),(c). Percebe-se que dos dois cenrios o que mais se aproxima da distribuio alvo o segundo.


Este cenrio no trabalha com uma estatstica suficiente como funo sumrio, mas sim uma partedela. Como falado anteriormente, a suficincia, utilizada como sumrio, ou parte da suficincia,aliada com a escolha da funo de distncia, deve gerar boas amostras. imprescindvel investigarcada caso a fim de determinar a melhor combinao possvel. Quando as estatsticas suficientesforam usadas como funo sumrio, primeiro cenrio, amostras y suficientemente prximas dasamostras x (passo 6 do algoritmo 5) so geradas em grande quantidade, o que diminui a taxa derejeio de . Tal fato no acontece no segundo cenrio, que utiliza uma outra estatstica, nosuficincia e, dificulta, restringe, a gerao de amostras y prximas de x aumentando assim a taxade rejeio de . Esse fato pode explicar a razo da grande gerao de pontos em uma regio debaixa probabilidade a posteriori.

De modo geral, as suspeitas sobre a influncia dos parmetros de tunning foram comprovadas.Autilizao de estatsticas suficientes, ou parte delas, como resumo das amostras, associada com adistncia euclidiana fazendo o papel de (), tende a uma amostra mais plausvel com a distribuioalvo. A tolerncia influencia a velocidade computacional e a qualidade da amostra como prevamos.

O fato do ABC gerar pontos amostrais regies de baixa probabilidade a posteriori comentadoem Marin et al. (2011). O MCMC-ABC , j exposto anteriormente foi criado com o intuito deresolver esse problema.

Como os modelos estatsticos apresentados nos captulos seguintes tem uma complexidade com-putacional elevada, decidiu-se implementar somente o ABC-Rejeio para no elevar inda mais acomplexidade dos clculos. Por esse motivo, os desempenho do MCMC-ABC no foi apresentadoaqui. mas os dados sobre ele esto disponveis em Campos e Wechsler (2012).

2.5 Comparando o ABC com um algoritmo MCMCApenas com as informaes apresentadas na seo anterior, no se pode afirmar se o ABC

um bom amostrador, se gera amostras confiveis e se de rpida execuo computacional. Por isso necessrio compar-lo com a classe de algoritmos MCMC .

Foi decidido comparar as amostras geradas, pelo ABC, para o segundo toy-model, modelonormal-gama, com as amostras geradas por um algoritmo MCMC. O escolhido para este modelofoi o amostrador de Gibbs (veja Casella e George, 1992).

Como mostrado anteriormente o modelo normal-gama apresenta as seguintes distribuies aposteriori,

|X = x Gama(n, n),|, X = x Normal(n, [(n+ ) ]1);

deste modo,

(, |X = x) =

(n+ )

2e

12

(n)2(n+) n(n)

n1en . (2.15)

No difcil ver que as distribuies condicionais completas, (|,X = x) e ( |,X = x),que so necessrias para implementar o amostrador de Gibbs, coincidem com as distribuies aposteriori de |X = x e |, X = x. Utilizando a mesma amostra observada que foi utilizada naseo 2.4.2, pode-se obter osseguintes valores n = 51.53, n = 15.33598 e n = 0.4050541. Oprocesso para obter uma amostra da distribuio de , |X = x usando o amostrador de Gibbsencontra-se descrito abaixo, algoritmo 6.

Algoritmo 6 Amostrador de Gibbs (Normal-gama)

1: Escolher valores iniciais (0) e ((0)) e inicializar o contador j = 12: repeat3: Gerar (j) de Gama(51.53, 15.33598)4: Gerar (j) de Normal(0.4050541, (101 (j))1)5: j = j + 16: until Condio de convergncia

2.5 COMPARANDO O ABC COM UM ALGORITMO MCMC 15

Foi escolhido utilizar um critrio grfico para diagnosticar a convergncia da cadeia, monitorandoa mdia ergdica das duas cadeias, uma para e outra para . Foi estipulado que se a mdia ergdicana iterao j for igual a mdia ergdica da interao j 1, com pelo menos 4 casas decimais deigualdade, e se essa igualdade se mantiver h por pelo menos 500 iteraes, as cadeias tero atingidoa convergncia.

Por um motivo de comparao, decidiu-se iniciar sempre um novo par de cadeias deMarkov paracada nova observao da amostra gerada. Deste modo possvel comparar o nmero de interaesque o ABC necessitou para gerar uma observao com o nmero de iteraes que o amostrador deGibbs necessitou para gerar uma observao, descontando as 500 iteraes obrigatrias para a cons-tatao da convergncia, alm de no haver preocupao com a dependncia entre as observaes.

Tabela 2.1: Sumrio do nmero de interaes necessrias para gerar uma observao usando o amostradorde Gibbs

Mnimo 1o Quartil Mediana Mdia 3o Quartil Mximo Varincia858 1249 1367 1386 1504 2334 34778,59

Ao comparar a informao apresentada na tabela 2.1 com as informaes contidas na seo 2.4.2,nota-se que, desconsiderando as 500 interaes obrigatrias para a constatao de convergncia,o amostrador de Gibbs necessitou, em mdia, de 1386 interaes para gerar uma observao daamostra da distribuio de interesse, enquanto o ABC, no melhor dos casos necessitou de apenas299 interaes, em mdia, para gerar uma observao.

Comparando os dois algoritmos apenas por esse critrio, tem-se a impresso de que o ABC um algoritmo melhor, mas essa deciso pode estar equivocada. Tambm necessrio comparara qualidade da amostra gerada por cada um dos amostradores. As amostras geradas pelo ABC,primeiro e terceiro cenrio, e Gibbs foram comparadas, nas figuras 2.4 e 2.5. Estes cenrios, apre-sentados na seo 2.4.2, foram escolhidos por terem gerado as amostras que mais se aproximaramcom a distribuio alvo (ver figura 2.2).

Com o auxlio das figuras 2.4 e 2.5 pode-se dizer que, para o modelo normal-gama, o amostradorde Gibbs conseguiu gerar melhores amostras, embora necessitasse de um pouco mais de 6 vezesmais interaes que o ABC no cenrio que gerou a melhor amostra, mas ainda assim o desempenhocomputacional do amostrador de Gibbs considerado bom. Os grficos na figura 2.4 confirmam oque j foi dito, o algoritmo ABC-rejeio tende a gerar pontos em baixa rea de probabilidade aposteriori. O amostrador de Gibbs conseguiu gerar de maneira satisfatria a cauda da distribuioalvo.

Note que, mesmo com as escalas dos grficos na figura 2.5 sendo diferentes, possvel ver queas amostras gerada pelo ABC mais dispersas que a amostra gerada pelo amostrador de Gibbs, umaexplicao para isso pode ser a tendncia do ABC-Rejeio de gerar pontos em regies de baixaprobabilidade a posteriori.

Se fosse necessrio escolher um algoritmo como o melhor seria o amostrador de Gibbs poismostrou-se mais confivel quando foi utilizado. Deve-se lembrar que este modelo em especifico simples, e j era esperada essa derrota do ABC. Mas para casos com uma funo de verossimi-lhana de forma analtica intratvel, o ABC se mostrar um algoritmo mais eficiente do que vistoat aqui, para gerar observaes de modelos em estudo, como ser mostrado nos prximos captulos.


(a)

Den

sity

0.0 0.2 0.4 0.6 0.8 1.0

0.0

0.5

1.0

1.5

2.0

2.5

(b)

Den

sity

0.0 0.2 0.4 0.6 0.8 1.0

0.0

0.5

1.0

1.5

2.0

2.5

(c)

Den

sity

0.0 0.2 0.4 0.6 0.8 1.0

0.0

0.5

1.0

1.5

2.0

2.5

(d)

Den

sity

0.0 0.2 0.4 0.6 0.8 1.0

0.0

0.5

1.0

1.5

2.0

2.5

Figura 2.1: (a) Histograma da amostra mdia gerada no primeiro cenrio comparada com a densidade aposteriori (equao 2.12); (b) Histograma da amostra mdia gerada no segundo cenrio comparada com adensidade a posteriori (equao 2.12); (c) Histograma da amostra mdia gerada no terceiro cenrio com-parada com a densidade a posteriori (equao 2.12); (d) Histograma da amostra mdia gerada no quartocenrio comparada com a densidade a posteriori (equao 2.12)


(a)

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

0.9

1

0 1 2 3 4 5

05

1015

20

(b)

0.1 0.2 0.3 0.4 0.5 0.6

0.7 0.8

0.9

1

0 1 2 3 4 5

05

1015

20

(c)

0.1 0.2 0.3 0.4

0.5 0.6 0.7

0.8 0.9

1

0 1 2 3 4 5

05

1015

20

Figura 2.2: (a) Funo de probabilidade acumulada de (, )|x; (b) Funo emprica mdia da amostra doprimeiro cenrio; (c) Funo emprica mdia da amostra do segundo cenrio.


1.5 1.0 0.5 0.0 0.5

05

1015

20

(a)

2 1 0 1 2

05

1015

(b)

Figura 2.3: (a) Grfico de disperso de todas as amostras geradas no primeiro cenrio; (b) Grfico dedisperso de todas as amostras geradas no segundo cenrio.


(a)

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

0.9

1

0 1 2 3 4 5

05

1015

20

(b)

0.1 0.2 0.3 0.4 0.5 0.6

0.7 0.8

0.9

1

0 1 2 3 4 5

05

1015

20

(c)

0.1 0.2 0.3 0.4

0.5 0.6 0.7

0.8 0.9

1

0 1 2 3 4 5

05

1015

20

(d)

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

0.9

1

0 1 2 3 4 5

05

1015

20

Figura 2.4: (a) Funo de probabilidade acumulada de (, )|x; (b) Funo acumulada emprica mdia daamostra do primeiro cenrio; (c) Funo acumulada emprica mdia da amostra do segundo cenrio; (d)Funo acumulada emprica mdia da amostra gerada pelo amostrador de Gibbs.


1.5 1.0 0.5 0.0 0.5

05

1015

20

(a)

2 1 0 1 2

05

1015

(b)

0.2 0.3 0.4 0.5 0.6

23

45

(c)

Figura 2.5: (a) Grfico de disperso de todas as amostras geradas no primeiro cenrio; (b) Grfico dedisperso de todas as amostras geradas no segundo cenrio; (c) Grfico de disperso de todas as amostrasgeradas pelo amostrador de Gibbs.

Captulo 3

Gerando observaes de uma distribuio de Holgate

3.1 IntroduoNo presente captulo apresentado um outro Toy Model, mais complexo que os apresentados no

captulo anterior. Esse novo exemplo, considera os dados observados obedecendo a distribuio deHolgate. Essa complexidade consequncia de uma maior sofisticao matemtica, que exemplificade forma mais clara o fato de que o ABC no necessita da funo de verossimilhana para geraruma amostra de uma determinada distribuio a posteriori.

Esse grau complexidade, apresentada e analisada no decorrer desse captulo, dificulta a imple-mentao de algum algoritmo para gerar observaes desta distribuio. Ser mostrado como oalgoritmo de Metropolis-Hastings (Gamerman e Lopes, 2006) pode ser usado para tal tarefa. Comonos exemplos passados, o ABC ser posto em ao, e consequentemente, as amostras geradas porcada um dos dois algoritmos sero comparadas.

3.2 Distribuio de Poisson bivariadaExistem diferentes formas de representar a distribuio Poisson bivaraida. A que ser tratada

aqui, distribuio de Holgate, como dito em Chire (2013), construdo pelo mtodo da reduotrivariada. Basicamente, esse mtodo consiste em reduzir a dimenso de um vetor, de trs dimensespara duas, atravs de uma operao nas componentes do vetor original, aqui essa operao ser asoma.

A construo de uma distribuio Holgate feita da seguinte forma:

Considere U1, U2, U3 variveis aleatrias (independentes) Poisson com parmetros 1, 2, 3,respectivamente;

Seja X1 = U1 +U3 e X2 = U2 +U3. O par (X1, X2) tem distribuio Holgate com parmetros(1, 2, 3) (ou (1 + 3, 2 + 3, 3) como visto em Stern e Zacks (2002));

A funo de probabilidade do par (X1, X2) da forma:

P (X1 = x, X2 = y|1, 2, 3) = e(1+2+3)min(x,y)i=0

xi1 yi2

i3

(x i)!(y i)!i!(3.1)

Notem que, quando U3 igual a zero, as componentes do par (X1, X2) so independentes.As distribuies de probabilidade marginais das variveis X1 e X2 so da forma

P (X1 = x|1, 3) =e(1+3)(1 + 3)

x

x!, (3.2)

P (X2 = x|2, 3) =e(2+3)(2 + 3)

x

x!, (3.3)

i. e., palavras X1 Poisson(1 + 3) e X2 Poisson(2 + 3). Uma consequncia desse fato acovarincia de X1 e X2 que igual a 3, COV (X1, X2) = 3 (Chire, 2013, pag. 11).

COV (X1, X2) = COV (U1 + U3, U2 + U3)

21

22 GERANDO OBSERVAES DE UMA DISTRIBUIO DE HOLGATE 3.2

= COV (U1 + U3, U2) + COV (U1 + U3, U3)

= COV (U1 + U2) + COV (U1 + U3) + COV (U3, U2) + COV (U3, U3)

= COV (U3, U3)

= V ar(U3)

= 3.

Portanto 3 uma mdida de dependncia entre X1 e X2 e,

X1,X2 =3

(1 + 3)(2 + 3). (3.4)

Note que 0 X1,X2 1, pois 1, 2, 3 > 0 e, alm disso, a funo de probabilidade do par(X1, X2), quando X1,X2 = 0, da seguinte forma,

P (X1 = x,X2 = y) =x1

y2e(1+2)

x!y!, x, y = 0, 1, . . . (3.5)

Considere agora um conjunto de observaes independentes x = [(x1, x2)1, (x1, x2)2, . . . , (x1, x2)n].A funo de verossimilhana dada da seguinte forma:

L(, x) = en(1+2+3)nj=1

min(xj ,yj)i=0

xji1

yji2

i3

(xj i)!(yj i)!i!, (3.6)

onde = (1, 2, 3). Em Chire (2013) a funo de verossimilhana apresentada com umadiferente parametrizao, em funo das esperanas de X1 e X2. Sejam 1 = 1 + 3 = E[X1] e2 = 2 +3 = E[X2]. A funo de verossimilhana, parametrizada em funo de = (1, 2, 3) da forma:

L(, x) = en(1+2+3)nj=1

min(xj ,yj)i=0

(1 3)xji(2 3)yjii3(xj i)!(yj i)!i!

. (3.7)

Essa reparametrizao, originalmente apresentada em Holgate (1964), de acordo com o prprio,torna o processo de estimao via mxima verossimilhana matematicamente mais simples, masainda assim continua bastante trabalhoso.

Todo o processo de estimao pontual, via mtodo dos momentos e mxima verossimilhana, apresentado em Chire (2013), bem como toda a metodologia para teste de hipteses, via teste darazo de verossimilhana, referente a hiptese de independncia entre X1 e X2, isto , 3 = 0.

As estimativas pontuais de , como mostradas em Chire (2013), via mtodo dos momentos soas seguintes:

1 = X1 m1,1, (3.8)2 = X2 m1,1, (3.9)3 = m1,1, (3.10)

onde m1,1 = 1n

(i,j)D

(iX1)(jX2)nij , nij =n

k=1 I{(i,j)}((X1, X2)k) a frequncia de cada par

(i, j), para i, j = 0, 1, . . ., e D = {(i, j) : nij > 0} representa o conjunto de dados observados. Porconsequncia,

n =

(i,j)D

nij .

J os estimadores via mxima verossimilhana, mais trabalhosos e complicados de serem obtidos,

3.3 ABORDAGEM BAYESIANA 23

so:

1 = X1 3, (3.11)2 = X2 3. (3.12)

O estimador de mxima verossimilhana para 3 obtido atravs da raiz da seguinte equao, comrelao a 3,

1

n

ni=1

p(xi 1, yi 1)p(xi, yi)

= 1, (3.13)

utilizando 1 = X1 e 2 = X2, no intervalo [0,min(X1, X2)]. A funo p(xi, yi) a funo deprobabilidade P (X1 = xi, X2 = yi) utilizando a reparametrizao presente na equao 3.7;

3.3 Abordagem BayesianaEm Arruda (2000) h a necessidade de modelar, probabilisticamente, resultados de partidas de

futebol, com a finalidade de poder prever os resultados das partidas seguintes, vitria, empate ouderrota do time com o mando de campo.

O placar do jogo pode ser encarado como uma varivel aleatria de duas dimenses, com cadadimenso, marginalmente, seguindo uma distribuio univariada de Poisson. Portanto, pode-seconsiderar que esse vetor, aleatrio, segue uma distribuio bivariada de Poisson. A distribuiode Holgate uma possvel distribuio de Poisson bivariada. Chire (2013) apresenta uma outraopo, a distribuio de Poisson generalizada bivariada.

O modelo Bayesiano abordado aqui, para estudar dados com essas caractersticas menciona-das acima, utiliza a distribuio de Holgate como distribuio dos dados observados, ou a seremobservados.

Considere as seguintes distribuies de probabilidade, a priori, para (1, 2, 3).

1 Gama(1, 1), (3.14)2 Gama(2, 2), (3.15)3 Gama(3, 3). (3.16)

Considere agora que 1, 2 e 3 so, a priori, independentes, e, tambm, a seguinte parametrizaoda distribuio gama. Seja x Gama(a, b), ento:

fX(x) =ba

(a)xa1ebx, a > 0, b > 0, x > 0.

Portanto, a distribuio a priori do trio ordenado (1, 2, 3)

(1, 2, 3) 111 e11212 e

22313 e33 . (3.17)

Utilizando o teorema de Bayes com a priori, equao 3.17, e a funo de verossimilhana,euqo 3.6, tem-se a seguinte distribuio de probabilidade a posteriori:

(|x) 111 212

313 e

(1(n+1)+2(n+2)+3(n+3))nj=1

min(xj ,yj)i=0

xj11

yji2

i3

(xj 1)!(yj i)!i!. (3.18)

Notem que o fatornj=1

min(xj ,yj)i=0

xj11

yji2

i3

(xj 1)!(yj i)!i! uma quantidade que, dificulta a obteno

da constante de normalizao da distribuio a posteriori, e por consequncia, deixa complexaqualquer tentativa de inferncia da distribuio a posteriori.

exatamente por esse motivo Chire (2013) apresenta uma abordagem alternativa, introduzindo


variveis latentes no modelo. Considere as variveis U1, U2 e U3 usadas para construir a distribuiode Holgate na seo 3.2 e o vetor aumentado (X1, X2, U3). A funo de probabilidade deste vetor dada por:

P (X1 = x,X2 = y, U3 = u|) = P (U1 = x u, U2 = y u, U3 = u|)

=xu1

yu2

u3e(1+2+3)

(x u)!(y u)!u!.

Note que x, y = 0, 1, . . . e u = 0, 1, . . . ,min(x, y).Suponha uma amostra aleatriaD = {(x1, y1, u1), . . . , (xn, yn, un)}. A funo de verossimilhana

baseada em D da seguinte forma:

L(D,) =ni=1

P (X1 = xi, X2 = yi, U3 = ui|)

=ni=1

xiui1 yiui2

ui3 e(1+2+3)

(xi ui)!(yi ui)!ui!

ni=1(xiui)

1 en1

ni=1(yiui)

2 en2

ni=1 ui

3 en3 (3.19)

Note que a funo de verossimilhana 3.19 tem uma forma mais simples, e amigvel, que afuno de verossimilhana 3.6, pois no possui o um produtrio de somatrios na sua expresso.

Utilizando o teorema de Bayes com a distribuio a priori 3.17 e a funo de verossimilhana 3.19tem-se a seguinte funo densidade de probabilidade a posteriori:

(|D) ni=1(xiui)+11

1 e(n+1)1

ni=1(yiui)+21

2 e(n+2)2

ni=1 ui+31

3 e(n+3)3 .

(3.20)A densidade da distribuio a posteriori 3.20 tem uma forma analtica mais amigvel que a

densidade a posteriori 3.18, mas ainda assim um problema para obter estimativas e realizar todoo processo inferencial de maneira direta, apenas utilizando a forma analtica densidade. Por issoh a necessidade de mtodos que geram observaes da distribuio de Holgate e assim realizar oprocesso inferencial to desejados.

3.4 Gerando observaesGeralmente algoritmos MCMC so os primeiros que vm a mente de um estatstico, e destes,

o mais simples, de um certo ponto de vista, o algoritmo de Gibbs (veja Casella e George, 1992).

3.4.1 Gibbs Sampler

Para implementar o algoritmo de Gibbs necessrio obter as formas analticas das distribuiescondicionais completas, e o mais importante, saber como gerar observaes delas. A partir dadistribuio a posteriori 3.18 as seguintes condicionais completas so obtidas:

(1|2, 3,x) 111 e1(n+1)

nj=1

min(xj ,yj)i=0

xj11

yji2

i3

(xj 1)!(yj i)!i!; (3.21)

(2|1, 3,x) 212 e2(n+2)

nj=1

min(xj ,yj)i=0

xj11

yji2

i3

(xj 1)!(yj i)!i!; (3.22)

(3|1, 2,x) 313 e3(n+3)

nj=1

min(xj ,yj)i=0

xj11

yji2

i3

(xj 1)!(yj i)!i!. (3.23)

Novamente tem-se a complicao gerada pelo produto de somas. Devido a essa complicao no possvel obter a constante de normalizao das distribuies condicionais completas.

3.4 GERANDO OBSERVAES 25

Os mtodos numricos de integrao, utilizados para calcular a constante de normalizao dacondicionais completas, no obtiveram sucesso e, consequentemente, no foi possvel obter-las, im-possibilitando a gerao observaes de 3.21, 3.22 e 3.23.

Chire (2013) apresenta uma configurao do amostrador de Gibbs utilizando a densidade 3.20,com dados aumentados. As distribuies condicionais completas so da forma:

1|2, 3, D Gama

(ni=1

(xi ui) + 1, n+ 1

), (3.24)

2|1, 3, D Gama

(ni=1

(yi ui) + 2, n+ 2

), (3.25)

3|2, 1, D Gama

(ni=1

ui + 3, n+ 3

). (3.26)

Com

P (u|x, y,) ni=1

1

(xi ui)!(yi ui)!ui!

(321

)ui. (3.27)

O algoritmo de Gibbs para gerar amostras da distribuio a posteriori da densidade 3.20, apre-sentado em Chire (2013), encontra-se a seguir.

Algoritmo 7 Amostrador de Gibbs para gerar observaes da distribuio Holgate

1: Escolher (0) = ((0)1 , (0)2 ,

(0)3 ) e inicializar o contador j = 0;

2: Faa t = t+ 1;3: Gerar u(t) a partir de 3.27.4: Gerar (t)1 ,

(t)2 e

(t)3 , a partir de 3.24, 3.25 e 3.26 respectivamente.

5: Repetir os passos 2, 3 e 4 at a convergncia ser diagnosticada.

3.4.2 Metropolis-Hastings

O algoritmo de Metropolis-Hastings (Veja Gamerman e Lopes, 2006, cap. 6), no caso abordadoneste captulo, tem uma vantagem sobre o algoritmo de Gibbs. A despreocupao com o clculoda constante de normalizao das equaes 3.21 3.22 3.23. Essa vantagem oriunda da prpriaestrutura do algoritmo.

Algoritmo 8 Metropolis-Hastings

1: Escolher (0) e inicializar o contador j = 12: Mover a cadeia para um novo valor gerado a partir da densidade q((j1), ).3: Calcular a probabilidade de aceitao da mudana da cadeia, ((j1), ). Se aceito, (j) = ,

caso contrario, j = (j1) e a cadeia no se move.4: Atualizar o contador j para j + 1 e retornar at o passo 2 at a convergncia ser atingida.

A quantidade ((j1), ) do algoritmo 8 representa a probabilidade de que seja aceito comoum novo elemento da trajetria da cadeia de Markov, sendo calculado pela seguinte expresso:

(, ) = min

{1,()q(, )

()q(, )

}. (3.28)

A funo q(, ) chamada de kernel de transio, uma verso da matriz de transio(espao deestados discreto) para quando o espao de estados contnuo. Representa, como o prprio nome jdiz, a transio da cadeia de Markov do estado para o estado . A funo () representa a densi-dade alvo, da qual h interesse em gerar uma amostra. Como h uma razo de funes, ()/(),tudo o que no depender de e ser simplificado, incluindo a constante de normalizao.


Dentre os casos especiais do algoritmo deMetropolis-Hastings apresentados em Gamerman e Lopes(2006) foram geradas amostras de teste usando os caso especiais Random Walk Chain e Indepen-dence Chain.

Para implementar o algoritmo de Metropolis-Hastings, atravs da tcnica Random Walk Chainos resultados das equaes 3.21 a 3.23 foram aproveitados, decidiu-se evoluir uma cadeia para cadaum dos i. Uma cadeia influenciar no na outra atravs do calculo de (, ).

Como comentado em Gamerman e Lopes (2006), no algoritmo de Metropolis-Hastings via ran-dom walk chain, a evoluo da cadeia de Markov ocorre da seguinte forma:

(j) = (j1) + wj , (3.29)

onde wj uma varivel aleatria com distribuio independente da cadeia. Geralmente as variveiswj so independentes e identicamente distribudas com densidade fw. Deste modo o kernel detransio para este caso especial q(, ) = fw( ).

Decidiu-se criar uma cadeia de Markov para cada um dos s e, portanto, foram definidas trsdistribuies para os rudos, w(1)j Normal

(0,

1+ni=1 xi

(1+n)2

), w(2)j Normal

(0,

2+ni=1 yi

(2+n)2

)e w(3)j Normal

(0, 3+cov(x,y)

(3+n)2

), onde um parmetro de tunning que ajusta a velocidade

de convergncia em oposio a taxa de aceitao do algoritmo. As varincias das distribuiesdos rudos foram escolhidas da maneira apresentada acima, pelo motivo de restringir o espao devariao das cadeias, baseando-se na amostra observada. Espera-se que isso melhore o desempenhocomputacional do algoritmo.

Pelo fato das distribuies de w(1)j , w(2)j e w

(1)j serem simtricas em torno do valor zero, h

uma simetria nas cadeias de Markov, portanto, qi(, ) = qi(, ) e i(, ) = min{

1,i ()

i ()

},

onde as distribuies i() referem-se as equaes 3.21 a 3.23.A outra tcnica explorada foi oMetropolis-Hastings via Independence Chain, onde o estado atual

da cadeia foi gerado de forma independente do estado anterior. Como Gamerman e Lopes (2006)explica, apesar da cadeia evoluir com independncia entre os estados, o algoritmo de Metropolis-Hastings ainda apresenta a propriedade markoviana, pois a probabilidade, (, ), de aceitar ou noum valor como parte da cadeia no atual instante, ainda depende do estado da cadeia no instanteanterior.

Para implementar o algoritmo de Metropolis-Hastings via Independence decidiu-se evoluir umacadeia tridimensional para o trio ordenado = (1, 2, 3), da seguinte maneira:

f(l1, l2, l3) =11

(1)l111 e

1l1 22

(2)l212 e

2l2 33

(3)l313 e

3l3 . (3.30)

Deste modo, (, ) = min{

1,()

3i=1

i1ei

()3i=1

i1ei

}, onde () refere-se a distribuio 3.18.

O critrio de convergncia adotado ser especificado mais adiante, mas pode-se adiantar quebaseia-se na mdia ergdica da cadeia. Se a mdia ergdica manter-se, com um certo nvel detolerncia, estvel por um tempo razoavelmente longo, considera-se que a cadeia convergiu.

Os pseudo-cdigos usados para implementar o Random Walk Chain Metropolis-Hastings e In-dependence Chain Metropolis-Hastings so apresentados a seguir, algoritmos 9 e 10

3.4.3 ABC

Os detalhes sobre o funcionamento do ABC j foram expostos no captulo 2, discutindo-seaqui apenas adaptaes especificas do ABC para gerar observaes para o modelo da distribuioHolgate.

Apesar de observar uma varivel de contagem-placar de um jogo-foi implementada a verso doABC para variveis contnuas, algoritmo 2. Tal fato motivado por critrios computacionais. Emsimulaes preliminares foi constatado que o algoritmo necessitava de um nmero bastante elevadode iteraes para poder aceitar um elemento gerado como parte de uma amostra de interesse. Ao


Algoritmo 9 Random Walk Chain Metropolis-Hastings1: Iniciar contador j = 12: Iniciar as cadeias (0)1 ,

(0)2 ,

(0)1

3: Gerar 1 = (j1)1 + w

1j , onde w

1j Normal

(0,

1+ni=1 xi

(1+n)2

)4: Gerar u de uma uniforme(0, 1)5: if u < 1(

(j1)1 ,

1) then

6: (j)1 1

7: else8:

(j)1

(j1)1

9: end if10: Gerar 2 =

(j1)2 + w

2j , onde w

2j Normal

(0,

2+ni=1 xi

(2+n)2


(j1)2 ,

2) then

13: (j)2 2

14: else15:

(j)2

(j1)2

16: end if17: Gerar 3 =

(j1)3 + w

3j , onde w

3j Normal

(0,

3+ni=1 xi

(3+n)2


(j1)3 ,

3) then

20: (j)3 3

21: else22:

(j)3

(j1)3

23: end if24: Repetir os passos 3 a 23 at diagnosticar a convergncia

Algoritmo 10 Independence Chain Metropolis-Hastings1: Iniciar contador j = 12: Iniciar a cadeia (0)

3: Gerar a partir de 3.304: Gerar u de uma uniforme(0, 1)5: if u < ((j1),) then6: (j) 7: else8: (j) (j1)9: end if

10: Repetir os passos 3 a 9 at diagnosticar a convergncia


usar a verso do ABC para variveis contnuas, a taxa de rejeio diminui, poupando assim recursoscomputacionais, em detrimento da qualidade da amostra gerada. Espera-se que alguns pontos aceitospelo ABC-Rejeio, para variveis contnuas, no sejam aceitos pela verso do ABC-Rejeio paravariveis discretas, o que implica numa aproximao no to boa quanto a desejada.

O vetor de sumrio de amostra ser um vetor com a primeira componente igual mdia dasprimeiras componentes das amostras. A segunda componente do vetor de sumrio a mdia dassegundas componentes da amostra, e por fim, a ltima componente do vetor de sumrio a co-varincia da primeira e segunda componentes da amostra((x) = (x1,,x2,, cov(x1,,x2,))). Talescolha foi motivada pela facilidade em obter tais quantidades, sendo mais simples de calcular queos estimadores de mxima verossimilhana e mtodo dos momentos.

Note que nesse toy model (x) no uma estatstica suficiente para (1.2, 3). A tarefa deencontrar um estimador suficiente para no nada simples, tendo em vista a complexidade dafuno de probabilidade,e por consequncia da funo de verossimilhana.

Definiu-se a funo distncia () como sendo a distancia euclidiana entre os dois vetores desumrio, o da amostra observada (x) e o da amostra gerada (y).

O algoritmo ABC, configurado para gerar amostras de uma distribuio de Holgate est apre-sentado a seguir, algoritmo 11

Algoritmo 11 ABC - Rejeio 2 (Holgate)1: for i = 1 to N do2: repeat3: Gerar 1 de Gama(1, 1)4: Gerar 2 de Gama(2, 2)5: Gerar 3 de Gama(3, 3)6: Gerar u1 de Poisson(1)7: Gerar u2 de Poisson(2)8: Gerar u3 de Poisson(3)9: y = (u1 + u3, u2 + u3)

10: until ((x) (y)) 11: i =

12: end for

3.4.4 Comparando mtodos

Com a finalidade de comparar o desempenho dos trs algoritmos, 9, 10 e 11, foram geradasobservaes da mesma distribuio a posteriori.

Para tal distribuio, foram escolhidos os hiperparmetros da seguinte forma: (1, 2, 3) =(1, 2, 3) = (1, 1, 1). Como amostra observada, x considere a sequncia de pares encontradosem Arruda (2000, pag. 12),

x = ((1, 0), (1, 2), (2, 3), (3, 2), (3, 0), (2, 1)),

considere tambm que tal amostra consiste de observaes permutveis.O critrio usado para diagnosticar a convergncia das cadeias de Markov, algoritmos 9 e 10, foi

usar a mdia ergdica da cadeia. No algoritmo que utiliza a metodologia de passeio aleatrio as, trscadeias, uma para cada i, as trs mdias ergdicas devem se manter iguais, com uma tolernciade 105, por 500 iteraes. J o algoritmo que usa a metodologia de cadeias independentes, amdia ergdica da cadeia, referente a deve se manter igual com uma tolerncia de 105, por 500iteraes.

O ABC implementado para o modelo Holgate no necessita de um critrio de convergncia,mas necessita de uma tolerncia, o quanto a amostra gerada est prxima da amostra observada.Sabemos que essa tolerncia influencia na velocidade do algoritmo e pensando nisso foi impostauma tolerncia de 105 para manter uma certa justia entre os algoritmos trabalhados, mas em


testes preliminares usar um to pequeno revelou-se impraticvel. Este fato expe uma fraquezado ABC. Por fim houve a necessidade da escolher um maior. Neste caso, a tolerncia escolhidafoi 102.

O amostrador de Gibbs no foi implementado por dois motivos. O primeiro que ele j implementado e analisado em Chire (2013). O segundo motivo que ainda no claro como gerarvalores da distribuio de u|x, y,, equao 3.27. Foi considerado que a complexidade de gerarobservaes de u|x, y, equivalente a efetuar clculos na funo de verossimilhana 3.6.

Para os trs algoritmos foram geradas 10 amostras de tamanho 1000. Para os algoritmos baseadosem MCMC, a cada observao gerada, era iniciada uma nova cadeia, com estados iniciais geradosaleatoriamente, para assegurar a independncia das observaes da amostra final.

Algumas estatsticas descritivas do nmero de iteraes que cada algoritmo levou para geraruma observao do modelo Holgate esto dispostas na tabela 3.1. possvel ver que a metodologiaMCMC, baseada em passeio aleatrio, aparenta ser a mais eficiente na gerao de amostras, pois a metodologia mais rpida, necessitando de menos iteraes, seguida pela metodologia do ABC edo MCMC, com em cadeias independentes.

Tabela 3.1: Sumrio do nmero de iteraes necessrias para gerar uma observao da distribuio deHolgate usando os algoritmos 9, 10 e 11

Mnimo 1o Quartil Mediana Mdia 3o Quartil Mximo VarinciaRandom Walk MCMC 508 547 560 564,2 575 371 634,592Independence Chain MCMC 9588 15489 19648 21689 26214 55365 60560197ABC 1 1500 3688 5409 7456 57340 30582628

Os valores presente na coluna varincia da tabela 3.1 podem parecer confusos. Como a varinciado nmero de iteraes do algoritmo ABC pode ser maior que a varincia nmero de iteraes doalgoritmo Independence Chain MCMC, a amplitude do nmero de iteraes do ABC maior? Aresposta para essa pergunta encontra-se nas informaes que podem ser extradas no grfico dafigura 3.1.

possvel ver que 75% das iteraes do ABC foram menores que o valor mnimo das iteraesdo Independence Chain MCMC. Deste modo, a ordem de grandeza do nmero de iteraes doIndependence Chain MCMC bem maior que a ordem de grandeza do nmero de iteraes doABC, o que explica o fato da varincia do nmero da iteraes do ABC quase a metade davarincia do nmero da iteraes do Independence Chain MCMC, mesmo tendo uma amplitudemaior.

Na figura 3.2 esto apresentados os histogramas de 1, 2 e 3, cujas amostras foram geradas pelametodologia de MCMC via passeio aleatrio. As 10 amostras de tamanho 1000 foram consideradascomo uma nica amostra de tamanho 10000, como feito nos toy models do captulo 2. Percebeu-seque h uma quantidade demasiadamente grande de observaes na cauda direita dos histogramas,indicando uma bimodalidade das distribuies.

Este comportamento bimodal dos histogramas pode ser indcio de uma dependncia entre 1,2 e 3, em que 2 e 3 dependem de 1, o que contraria as crenas a priori sobre . No existemmotivos para acreditar que 1, 2 e 3 apresentem uma dependncia a posteriori.

Uma outra explicao para este fato advm da taxa de aceitao do algoritmo. Durante o tempoem que o algoritmo 9 estava em execuo percebeu-se que as cadeias passeavam algumas regiesespecificas, restringindo sua variao a essas regies, ocasionando poucas mudanas na trajetriada cadeia.

Por fim, o estimador de Bayes, considerando essa amostra, de = (1, 2, 3), sob perdaquadrtica, = (0.50390, 1.66900, 2.13800).

J na figura 3.3 observamv-se os histogramas de 1, 2 e 3, considerando a amostra gerada peloalgoritmo 10. possvel ver que neste caso j no constatada a bimodalidade nos histogramas de2 e 3 que so apresentados na figura 3.2, o que pode indicar que esse fenmeno pode ter acontecidodevido ao fato de que no algoritmo 9 atualiza primeiro a cadeia referente a 1 e depois as cadeias


subsequentes, j com o algoritmo 10 h apenas uma cadeia sendo evoluda, e assim 1, 2 e 3 soatualizados conjuntamente, no havendo a possibilidade de uma componente de (1, 2, 3) ficarestagnada em um estado especfico.

Os histogramas apresentam um comportamento semelhante ao de distribuio Gama, o queest dentro do esperado, j que a distribuio a priori de uma distribuio trivariada cujasdistribuies marginais so Gamas independentes. Finalizando os comentrios sobre a amostragerada pelo algoritmo 10 apresentada a estimativa de Bayes, com relao a perda quadrtica, de como sendo o vetor (0.6417, 0.638200, 0.7361).

Por fim, na figura 3.4 esto os histogramas da amostra, de tamanho 10000, a posteriori de 1, 2e 3 gerada pelo algoritmo 11. Nota-se que h comportamento que lembra uma distribuio Gamapara os trs histogramas, como na figura 3.3, o que j era esperado. O vetor de estimativas de Bayes,a posteriori, com relao a perda quadrtica, para (1, 2, 3) (1.318548, 0.7509752, 0.6750414),respectivamente. Pode-se considerar que tanto este vetor de estimativas, quanto as estimativas daamostra apresentada na figura 3.3, so estimativas plausveis pra .

Perceba que as amostras gerada pelo ABC so mais dispersas que as amostras geradas peloMetropolis-Hastings via Random Walk Chain, mais um fato que pode ser atribudo tendnciaque o ABC-Rejeio tem de gerar observaes em regies com baixa probabilidade posteriori.Comparando as amostras geradas por esses dois mtodos, com a amostra gerada pelo algoritmo deGibbs, apresentado em Chire (2013, pag. 53), pode-se constatar que as dua amostras so coerentes,salvo as mudanas nos hiperparmetros.

Como as amostras geradas pelo ABC e pelo Independence Chain MCMC se mostraram bem se-melhantes, marginalmente, decidiu-se por construir os histogramas em trs dimenses dos s, dois adois, gerados por cada um dos mtodos. Nas figuras 3.5, 3.6 e 3.6 esto representados os histogramas,em trs dimenses, dos s gerados pelo Independence Chain MCMC. E nas figuras 3.8, 3.9 e 3.9esto representados os histogramas, em trs dimenses, dos s gerados pelo ABC. Ao compararos histogramas, de duas e trs dimenses, das amostras geradas pelos dois mtodos, possvel ver,no

Documents

Aplicações de Approximate Bayesian Computation a controle de