Aula 11 - Modelos Complexos
IESP-UERJ
Deep Learning
Vimos na aula anterior que técnicas como Random Forrest e Boosting permitem lidar com formas funcionais mais complexas, superando em poder de predição os modelos paramétricos. Estas técnicas são utilizadas quando temos dados estruturados. Quando temos dados não estruturados (imagens, sons, textos) outras técnicas apresentam desemepenhos ainda melhores. Estas técnicas envolvem o uso de neural networks, isto é, várias camadas de transformação dos inputs que os aproximam cada vez mais dos outputs esperados. Essas aproximações em várias camadas são feitas até que a distância entre o resultado obtido e o resultado esperado seja mínimo. Isto é, as transformações vão sendo ajustadas iterativamente, até que se aprenda o melhor ajuste aos dados. Esse aprendizado em várias camadas é chamado de Deep Learning.
[exemplo Keras]
Zero Inflated Model
Vimos que, embora modelos não paramétricos sejam muito bons para problemas de predição, sua baixa interpretabilidade é um obstáculo para a inferência. Os modelos lineares generalizados permitem uma certa flexibilidade. Vimos isso no caso de variáveis dicotômicas, multinomiais e de contagens. Nestes modelos substituímos funções de link em determinadas distribuições para estimar os parâmetros que aproximam nosso inputs no output esperado, isto é, os parâmetros que promovem o melhor ajuste aos dados.
Quando estimamos uma função de link trabalhamos com um modelo para a distribuição de um parâmetro que descreve nossa variável de interesse. Esse parâmetro, por sua vez, é descrito por outros parâmetros e esse processo pode continuar até onde o bom senso nos leve. Isso nos dá uma grande flexibilidade para modelar relações complexas entre variáveis.
Um tipo de relação complexa ocorre quando mais de um processo concorre para gerar nossos dados. Quando temos mais de uma causa para uma mesma observação temos o que se chama Modelos de Mistura que usa mais de uma verossimilhança para a mesma variável dependente. Isso ocorre, por exemplo, em modelos de contagem onde o valor zero de nossa VD pode significar a ausência de ocorrência ou que o proceso que leva à ocorrência não era possível. Pode ser, por exemplo, que uma abstenção eleitoral se deva ao não-comparecimento ou então ao não cancelamento do título de uma pessoa já falecida. Para modelar este tipo de processo de geração de dados primeiro fazemos um modelo para a geração de zeros e depois outro modelo para a geração da contagem de abstenções e então combinamos os dois modelos em uma função de verossimilhança.
Exemplo dos monges:
\[ \text{Probabilidae de se obter um zero: Pr}(Y = 0|p,\lambda) = \text{Pr}(beber|p) + \text{Pr}(trabalhar|p) X \text{Pr}(0|\lambda) = p + (1-p)\text{exp}(-\lambda) \\ \text{Probabilidae de se obter um livro ou mais: Pr}(Y = y_{i}|p,\lambda) = (1 - \pi) \frac{\lambda^{y_i} e^{-\lambda}} {y_i!}\]
No R
library(rethinking)## Carregando pacotes exigidos: rstan## Carregando pacotes exigidos: StanHeaders## Carregando pacotes exigidos: ggplot2## rstan (Version 2.21.5, GitRev: 2e1f913d3ca3)## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)## Carregando pacotes exigidos: cmdstanr## Warning: package 'cmdstanr' was built under R version 4.3.0## This is cmdstanr version 0.5.2## - CmdStanR documentation and vignettes: mc-stan.org/cmdstanr## - CmdStan path: /home/fernando/.cmdstan/cmdstan-2.29.1## - CmdStan version: 2.29.1## Carregando pacotes exigidos: parallel## rethinking (Version 2.21)##
## Attaching package: 'rethinking'## The following object is masked from 'package:rstan':
##
## stan## The following object is masked from 'package:stats':
##
## rstudentprob_beber <- 0.2
produz <- 1
# Dias no ano
N <- 365
# Simula os dias em que bebem
bebe <- rbinom(N, 1, prob_beber)
# Simula o numero de manuscritos completos.
y <- (1-bebe)*rpois(N,produz)
# Figura
simplehist(y, xlab = "Número de manuscritos", lwd = 4)
zero_bebe <- sum(bebe)
zero_trab <- sum(y==0 & bebe ==0)
zero_tot <- sum(y==0)
lines(c(0,0), c(zero_trab, zero_tot), lwd =4, col = "blue")
# Podemos estimar um modelo deste tipo usando as duas verossimilhanças
monge.zip <- ulam(
alist(
y|y>0 ~ custom(log1m(p) + poisson_lpmf(y|lambda)),
y|y==0 ~custom(log(p + (1-p)*exp(-lambda))),
logit(p) <- ap,
log(lambda) <- al,
ap ~ dnorm(0,1),
al ~ dnorm(1,0.5)
), data=list(y=as.integer(y)), chains = 4
)## Warning in '/tmp/RtmppCs3Rw/model-ce8d345e2067.stan', line 2, column 4: Declaration
## of arrays by placing brackets after a variable name is deprecated and
## will be removed in Stan 2.32.0. Instead use the array keyword before the
## type. This can be changed automatically using the auto-format flag to
## stanc## Running MCMC with 4 sequential chains, with 1 thread(s) per chain...
##
## Chain 1 Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 1 Iteration: 100 / 1000 [ 10%] (Warmup)
## Chain 1 Iteration: 200 / 1000 [ 20%] (Warmup)
## Chain 1 Iteration: 300 / 1000 [ 30%] (Warmup)
## Chain 1 Iteration: 400 / 1000 [ 40%] (Warmup)
## Chain 1 Iteration: 500 / 1000 [ 50%] (Warmup)
## Chain 1 Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 1 Iteration: 600 / 1000 [ 60%] (Sampling)
## Chain 1 Iteration: 700 / 1000 [ 70%] (Sampling)
## Chain 1 Iteration: 800 / 1000 [ 80%] (Sampling)
## Chain 1 Iteration: 900 / 1000 [ 90%] (Sampling)
## Chain 1 Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 1 finished in 0.2 seconds.
## Chain 2 Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 2 Iteration: 100 / 1000 [ 10%] (Warmup)
## Chain 2 Iteration: 200 / 1000 [ 20%] (Warmup)
## Chain 2 Iteration: 300 / 1000 [ 30%] (Warmup)
## Chain 2 Iteration: 400 / 1000 [ 40%] (Warmup)
## Chain 2 Iteration: 500 / 1000 [ 50%] (Warmup)
## Chain 2 Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 2 Iteration: 600 / 1000 [ 60%] (Sampling)
## Chain 2 Iteration: 700 / 1000 [ 70%] (Sampling)
## Chain 2 Iteration: 800 / 1000 [ 80%] (Sampling)
## Chain 2 Iteration: 900 / 1000 [ 90%] (Sampling)
## Chain 2 Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 2 finished in 0.2 seconds.
## Chain 3 Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 3 Iteration: 100 / 1000 [ 10%] (Warmup)
## Chain 3 Iteration: 200 / 1000 [ 20%] (Warmup)
## Chain 3 Iteration: 300 / 1000 [ 30%] (Warmup)
## Chain 3 Iteration: 400 / 1000 [ 40%] (Warmup)
## Chain 3 Iteration: 500 / 1000 [ 50%] (Warmup)
## Chain 3 Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 3 Iteration: 600 / 1000 [ 60%] (Sampling)
## Chain 3 Iteration: 700 / 1000 [ 70%] (Sampling)
## Chain 3 Iteration: 800 / 1000 [ 80%] (Sampling)
## Chain 3 Iteration: 900 / 1000 [ 90%] (Sampling)
## Chain 3 Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 3 finished in 0.2 seconds.
## Chain 4 Iteration: 1 / 1000 [ 0%] (Warmup)
## Chain 4 Iteration: 100 / 1000 [ 10%] (Warmup)
## Chain 4 Iteration: 200 / 1000 [ 20%] (Warmup)
## Chain 4 Iteration: 300 / 1000 [ 30%] (Warmup)
## Chain 4 Iteration: 400 / 1000 [ 40%] (Warmup)
## Chain 4 Iteration: 500 / 1000 [ 50%] (Warmup)
## Chain 4 Iteration: 501 / 1000 [ 50%] (Sampling)
## Chain 4 Iteration: 600 / 1000 [ 60%] (Sampling)
## Chain 4 Iteration: 700 / 1000 [ 70%] (Sampling)
## Chain 4 Iteration: 800 / 1000 [ 80%] (Sampling)
## Chain 4 Iteration: 900 / 1000 [ 90%] (Sampling)
## Chain 4 Iteration: 1000 / 1000 [100%] (Sampling)
## Chain 4 finished in 0.2 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 1.1 seconds.precis(monge.zip)## mean sd 5.5% 94.5% n_eff Rhat4
## ap -1.3124613 0.33016673 -1.8642471 -0.868278 624.4423 1.005613
## al 0.0581365 0.08086644 -0.0691358 0.183062 794.5027 1.004337plot(precis(monge.zip))
prior <- extract.prior(monge.zip, n = 1e4)## Warning in '/tmp/RtmppCs3Rw/model-ce8d67c92001.stan', line 2, column 4: Declaration
## of arrays by placing brackets after a variable name is deprecated and
## will be removed in Stan 2.32.0. Instead use the array keyword before the
## type. This can be changed automatically using the auto-format flag to
## stanc## Running MCMC with 1 chain, with 1 thread(s) per chain...
##
## Chain 1 Iteration: 1 / 20000 [ 0%] (Warmup)
## Chain 1 Iteration: 100 / 20000 [ 0%] (Warmup)
## Chain 1 Iteration: 200 / 20000 [ 1%] (Warmup)
## Chain 1 Iteration: 300 / 20000 [ 1%] (Warmup)
## Chain 1 Iteration: 400 / 20000 [ 2%] (Warmup)
## Chain 1 Iteration: 500 / 20000 [ 2%] (Warmup)
## Chain 1 Iteration: 600 / 20000 [ 3%] (Warmup)
## Chain 1 Iteration: 700 / 20000 [ 3%] (Warmup)
## Chain 1 Iteration: 800 / 20000 [ 4%] (Warmup)
## Chain 1 Iteration: 900 / 20000 [ 4%] (Warmup)
## Chain 1 Iteration: 1000 / 20000 [ 5%] (Warmup)
## Chain 1 Iteration: 1100 / 20000 [ 5%] (Warmup)
## Chain 1 Iteration: 1200 / 20000 [ 6%] (Warmup)
## Chain 1 Iteration: 1300 / 20000 [ 6%] (Warmup)
## Chain 1 Iteration: 1400 / 20000 [ 7%] (Warmup)
## Chain 1 Iteration: 1500 / 20000 [ 7%] (Warmup)
## Chain 1 Iteration: 1600 / 20000 [ 8%] (Warmup)
## Chain 1 Iteration: 1700 / 20000 [ 8%] (Warmup)
## Chain 1 Iteration: 1800 / 20000 [ 9%] (Warmup)
## Chain 1 Iteration: 1900 / 20000 [ 9%] (Warmup)
## Chain 1 Iteration: 2000 / 20000 [ 10%] (Warmup)
## Chain 1 Iteration: 2100 / 20000 [ 10%] (Warmup)
## Chain 1 Iteration: 2200 / 20000 [ 11%] (Warmup)
## Chain 1 Iteration: 2300 / 20000 [ 11%] (Warmup)
## Chain 1 Iteration: 2400 / 20000 [ 12%] (Warmup)
## Chain 1 Iteration: 2500 / 20000 [ 12%] (Warmup)
## Chain 1 Iteration: 2600 / 20000 [ 13%] (Warmup)
## Chain 1 Iteration: 2700 / 20000 [ 13%] (Warmup)
## Chain 1 Iteration: 2800 / 20000 [ 14%] (Warmup)
## Chain 1 Iteration: 2900 / 20000 [ 14%] (Warmup)
## Chain 1 Iteration: 3000 / 20000 [ 15%] (Warmup)
## Chain 1 Iteration: 3100 / 20000 [ 15%] (Warmup)
## Chain 1 Iteration: 3200 / 20000 [ 16%] (Warmup)
## Chain 1 Iteration: 3300 / 20000 [ 16%] (Warmup)
## Chain 1 Iteration: 3400 / 20000 [ 17%] (Warmup)
## Chain 1 Iteration: 3500 / 20000 [ 17%] (Warmup)
## Chain 1 Iteration: 3600 / 20000 [ 18%] (Warmup)
## Chain 1 Iteration: 3700 / 20000 [ 18%] (Warmup)
## Chain 1 Iteration: 3800 / 20000 [ 19%] (Warmup)
## Chain 1 Iteration: 3900 / 20000 [ 19%] (Warmup)
## Chain 1 Iteration: 4000 / 20000 [ 20%] (Warmup)
## Chain 1 Iteration: 4100 / 20000 [ 20%] (Warmup)
## Chain 1 Iteration: 4200 / 20000 [ 21%] (Warmup)
## Chain 1 Iteration: 4300 / 20000 [ 21%] (Warmup)
## Chain 1 Iteration: 4400 / 20000 [ 22%] (Warmup)
## Chain 1 Iteration: 4500 / 20000 [ 22%] (Warmup)
## Chain 1 Iteration: 4600 / 20000 [ 23%] (Warmup)
## Chain 1 Iteration: 4700 / 20000 [ 23%] (Warmup)
## Chain 1 Iteration: 4800 / 20000 [ 24%] (Warmup)
## Chain 1 Iteration: 4900 / 20000 [ 24%] (Warmup)
## Chain 1 Iteration: 5000 / 20000 [ 25%] (Warmup)
## Chain 1 Iteration: 5100 / 20000 [ 25%] (Warmup)
## Chain 1 Iteration: 5200 / 20000 [ 26%] (Warmup)
## Chain 1 Iteration: 5300 / 20000 [ 26%] (Warmup)
## Chain 1 Iteration: 5400 / 20000 [ 27%] (Warmup)
## Chain 1 Iteration: 5500 / 20000 [ 27%] (Warmup)
## Chain 1 Iteration: 5600 / 20000 [ 28%] (Warmup)
## Chain 1 Iteration: 5700 / 20000 [ 28%] (Warmup)
## Chain 1 Iteration: 5800 / 20000 [ 29%] (Warmup)
## Chain 1 Iteration: 5900 / 20000 [ 29%] (Warmup)
## Chain 1 Iteration: 6000 / 20000 [ 30%] (Warmup)
## Chain 1 Iteration: 6100 / 20000 [ 30%] (Warmup)
## Chain 1 Iteration: 6200 / 20000 [ 31%] (Warmup)
## Chain 1 Iteration: 6300 / 20000 [ 31%] (Warmup)
## Chain 1 Iteration: 6400 / 20000 [ 32%] (Warmup)
## Chain 1 Iteration: 6500 / 20000 [ 32%] (Warmup)
## Chain 1 Iteration: 6600 / 20000 [ 33%] (Warmup)
## Chain 1 Iteration: 6700 / 20000 [ 33%] (Warmup)
## Chain 1 Iteration: 6800 / 20000 [ 34%] (Warmup)
## Chain 1 Iteration: 6900 / 20000 [ 34%] (Warmup)
## Chain 1 Iteration: 7000 / 20000 [ 35%] (Warmup)
## Chain 1 Iteration: 7100 / 20000 [ 35%] (Warmup)
## Chain 1 Iteration: 7200 / 20000 [ 36%] (Warmup)
## Chain 1 Iteration: 7300 / 20000 [ 36%] (Warmup)
## Chain 1 Iteration: 7400 / 20000 [ 37%] (Warmup)
## Chain 1 Iteration: 7500 / 20000 [ 37%] (Warmup)
## Chain 1 Iteration: 7600 / 20000 [ 38%] (Warmup)
## Chain 1 Iteration: 7700 / 20000 [ 38%] (Warmup)
## Chain 1 Iteration: 7800 / 20000 [ 39%] (Warmup)
## Chain 1 Iteration: 7900 / 20000 [ 39%] (Warmup)
## Chain 1 Iteration: 8000 / 20000 [ 40%] (Warmup)
## Chain 1 Iteration: 8100 / 20000 [ 40%] (Warmup)
## Chain 1 Iteration: 8200 / 20000 [ 41%] (Warmup)
## Chain 1 Iteration: 8300 / 20000 [ 41%] (Warmup)
## Chain 1 Iteration: 8400 / 20000 [ 42%] (Warmup)
## Chain 1 Iteration: 8500 / 20000 [ 42%] (Warmup)
## Chain 1 Iteration: 8600 / 20000 [ 43%] (Warmup)
## Chain 1 Iteration: 8700 / 20000 [ 43%] (Warmup)
## Chain 1 Iteration: 8800 / 20000 [ 44%] (Warmup)
## Chain 1 Iteration: 8900 / 20000 [ 44%] (Warmup)
## Chain 1 Iteration: 9000 / 20000 [ 45%] (Warmup)
## Chain 1 Iteration: 9100 / 20000 [ 45%] (Warmup)
## Chain 1 Iteration: 9200 / 20000 [ 46%] (Warmup)
## Chain 1 Iteration: 9300 / 20000 [ 46%] (Warmup)
## Chain 1 Iteration: 9400 / 20000 [ 47%] (Warmup)
## Chain 1 Iteration: 9500 / 20000 [ 47%] (Warmup)
## Chain 1 Iteration: 9600 / 20000 [ 48%] (Warmup)
## Chain 1 Iteration: 9700 / 20000 [ 48%] (Warmup)
## Chain 1 Iteration: 9800 / 20000 [ 49%] (Warmup)
## Chain 1 Iteration: 9900 / 20000 [ 49%] (Warmup)
## Chain 1 Iteration: 10000 / 20000 [ 50%] (Warmup)
## Chain 1 Iteration: 10001 / 20000 [ 50%] (Sampling)
## Chain 1 Iteration: 10100 / 20000 [ 50%] (Sampling)
## Chain 1 Iteration: 10200 / 20000 [ 51%] (Sampling)
## Chain 1 Iteration: 10300 / 20000 [ 51%] (Sampling)
## Chain 1 Iteration: 10400 / 20000 [ 52%] (Sampling)
## Chain 1 Iteration: 10500 / 20000 [ 52%] (Sampling)
## Chain 1 Iteration: 10600 / 20000 [ 53%] (Sampling)
## Chain 1 Iteration: 10700 / 20000 [ 53%] (Sampling)
## Chain 1 Iteration: 10800 / 20000 [ 54%] (Sampling)
## Chain 1 Iteration: 10900 / 20000 [ 54%] (Sampling)
## Chain 1 Iteration: 11000 / 20000 [ 55%] (Sampling)
## Chain 1 Iteration: 11100 / 20000 [ 55%] (Sampling)
## Chain 1 Iteration: 11200 / 20000 [ 56%] (Sampling)
## Chain 1 Iteration: 11300 / 20000 [ 56%] (Sampling)
## Chain 1 Iteration: 11400 / 20000 [ 57%] (Sampling)
## Chain 1 Iteration: 11500 / 20000 [ 57%] (Sampling)
## Chain 1 Iteration: 11600 / 20000 [ 58%] (Sampling)
## Chain 1 Iteration: 11700 / 20000 [ 58%] (Sampling)
## Chain 1 Iteration: 11800 / 20000 [ 59%] (Sampling)
## Chain 1 Iteration: 11900 / 20000 [ 59%] (Sampling)
## Chain 1 Iteration: 12000 / 20000 [ 60%] (Sampling)
## Chain 1 Iteration: 12100 / 20000 [ 60%] (Sampling)
## Chain 1 Iteration: 12200 / 20000 [ 61%] (Sampling)
## Chain 1 Iteration: 12300 / 20000 [ 61%] (Sampling)
## Chain 1 Iteration: 12400 / 20000 [ 62%] (Sampling)
## Chain 1 Iteration: 12500 / 20000 [ 62%] (Sampling)
## Chain 1 Iteration: 12600 / 20000 [ 63%] (Sampling)
## Chain 1 Iteration: 12700 / 20000 [ 63%] (Sampling)
## Chain 1 Iteration: 12800 / 20000 [ 64%] (Sampling)
## Chain 1 Iteration: 12900 / 20000 [ 64%] (Sampling)
## Chain 1 Iteration: 13000 / 20000 [ 65%] (Sampling)
## Chain 1 Iteration: 13100 / 20000 [ 65%] (Sampling)
## Chain 1 Iteration: 13200 / 20000 [ 66%] (Sampling)
## Chain 1 Iteration: 13300 / 20000 [ 66%] (Sampling)
## Chain 1 Iteration: 13400 / 20000 [ 67%] (Sampling)
## Chain 1 Iteration: 13500 / 20000 [ 67%] (Sampling)
## Chain 1 Iteration: 13600 / 20000 [ 68%] (Sampling)
## Chain 1 Iteration: 13700 / 20000 [ 68%] (Sampling)
## Chain 1 Iteration: 13800 / 20000 [ 69%] (Sampling)
## Chain 1 Iteration: 13900 / 20000 [ 69%] (Sampling)
## Chain 1 Iteration: 14000 / 20000 [ 70%] (Sampling)
## Chain 1 Iteration: 14100 / 20000 [ 70%] (Sampling)
## Chain 1 Iteration: 14200 / 20000 [ 71%] (Sampling)
## Chain 1 Iteration: 14300 / 20000 [ 71%] (Sampling)
## Chain 1 Iteration: 14400 / 20000 [ 72%] (Sampling)
## Chain 1 Iteration: 14500 / 20000 [ 72%] (Sampling)
## Chain 1 Iteration: 14600 / 20000 [ 73%] (Sampling)
## Chain 1 Iteration: 14700 / 20000 [ 73%] (Sampling)
## Chain 1 Iteration: 14800 / 20000 [ 74%] (Sampling)
## Chain 1 Iteration: 14900 / 20000 [ 74%] (Sampling)
## Chain 1 Iteration: 15000 / 20000 [ 75%] (Sampling)
## Chain 1 Iteration: 15100 / 20000 [ 75%] (Sampling)
## Chain 1 Iteration: 15200 / 20000 [ 76%] (Sampling)
## Chain 1 Iteration: 15300 / 20000 [ 76%] (Sampling)
## Chain 1 Iteration: 15400 / 20000 [ 77%] (Sampling)
## Chain 1 Iteration: 15500 / 20000 [ 77%] (Sampling)
## Chain 1 Iteration: 15600 / 20000 [ 78%] (Sampling)
## Chain 1 Iteration: 15700 / 20000 [ 78%] (Sampling)
## Chain 1 Iteration: 15800 / 20000 [ 79%] (Sampling)
## Chain 1 Iteration: 15900 / 20000 [ 79%] (Sampling)
## Chain 1 Iteration: 16000 / 20000 [ 80%] (Sampling)
## Chain 1 Iteration: 16100 / 20000 [ 80%] (Sampling)
## Chain 1 Iteration: 16200 / 20000 [ 81%] (Sampling)
## Chain 1 Iteration: 16300 / 20000 [ 81%] (Sampling)
## Chain 1 Iteration: 16400 / 20000 [ 82%] (Sampling)
## Chain 1 Iteration: 16500 / 20000 [ 82%] (Sampling)
## Chain 1 Iteration: 16600 / 20000 [ 83%] (Sampling)
## Chain 1 Iteration: 16700 / 20000 [ 83%] (Sampling)
## Chain 1 Iteration: 16800 / 20000 [ 84%] (Sampling)
## Chain 1 Iteration: 16900 / 20000 [ 84%] (Sampling)
## Chain 1 Iteration: 17000 / 20000 [ 85%] (Sampling)
## Chain 1 Iteration: 17100 / 20000 [ 85%] (Sampling)
## Chain 1 Iteration: 17200 / 20000 [ 86%] (Sampling)
## Chain 1 Iteration: 17300 / 20000 [ 86%] (Sampling)
## Chain 1 Iteration: 17400 / 20000 [ 87%] (Sampling)
## Chain 1 Iteration: 17500 / 20000 [ 87%] (Sampling)
## Chain 1 Iteration: 17600 / 20000 [ 88%] (Sampling)
## Chain 1 Iteration: 17700 / 20000 [ 88%] (Sampling)
## Chain 1 Iteration: 17800 / 20000 [ 89%] (Sampling)
## Chain 1 Iteration: 17900 / 20000 [ 89%] (Sampling)
## Chain 1 Iteration: 18000 / 20000 [ 90%] (Sampling)
## Chain 1 Iteration: 18100 / 20000 [ 90%] (Sampling)
## Chain 1 Iteration: 18200 / 20000 [ 91%] (Sampling)
## Chain 1 Iteration: 18300 / 20000 [ 91%] (Sampling)
## Chain 1 Iteration: 18400 / 20000 [ 92%] (Sampling)
## Chain 1 Iteration: 18500 / 20000 [ 92%] (Sampling)
## Chain 1 Iteration: 18600 / 20000 [ 93%] (Sampling)
## Chain 1 Iteration: 18700 / 20000 [ 93%] (Sampling)
## Chain 1 Iteration: 18800 / 20000 [ 94%] (Sampling)
## Chain 1 Iteration: 18900 / 20000 [ 94%] (Sampling)
## Chain 1 Iteration: 19000 / 20000 [ 95%] (Sampling)
## Chain 1 Iteration: 19100 / 20000 [ 95%] (Sampling)
## Chain 1 Iteration: 19200 / 20000 [ 96%] (Sampling)
## Chain 1 Iteration: 19300 / 20000 [ 96%] (Sampling)
## Chain 1 Iteration: 19400 / 20000 [ 97%] (Sampling)
## Chain 1 Iteration: 19500 / 20000 [ 97%] (Sampling)
## Chain 1 Iteration: 19600 / 20000 [ 98%] (Sampling)
## Chain 1 Iteration: 19700 / 20000 [ 98%] (Sampling)
## Chain 1 Iteration: 19800 / 20000 [ 99%] (Sampling)
## Chain 1 Iteration: 19900 / 20000 [ 99%] (Sampling)
## Chain 1 Iteration: 20000 / 20000 [100%] (Sampling)
## Chain 1 finished in 0.2 seconds.dens(inv_logit(prior$ap))
dens(exp(prior$al))
post <- extract.samples(monge.zip)
mean(inv_logit(post$ap))## [1] 0.217132mean(exp(post$al))## [1] 1.063321Markov Chain Monte Carlo - MCMC
Para estimar os parâmetros do modelo de Poisson inflado em zeros usamos a função ulam do pacote rethinking. Na verdade esta função é só uma forma didática que McElreath usou para apresentar o Stan, um algoritmo desenvolvido por Gelman et al. para análise bayesiana que usa processos estocásticos para estimar parâmetros. No lugar de maximizar algo (como a verossimilhança) trata-se de explorar o espaço paramétrico de forma aleatória a partir de um ponto inicial. Processo deste tipo são conhecidos como Markov Chain Monte Carlo (MCMC).
Um algoritmo que implementa MCMC é o Metropolis. Esse algoritmo parte de um valor \(x\) de um parâmetro e sorteia-se um novo valor \(x1\) a partir de uma distribuição arbitrária. Calcula-se uma taxa de aceitação \(a = f(x)/f(x1)\) onde \(f(.)\) é a verossimilhança. Gera-se um número aleatório de uma distribuição uniforme e se esse número for menor que \(a\) atualizamos o valor do parâmetro para \(x = x1\), do contrário mantemos \(x\).Esse processo faz com que se explore todo o espaço paramétrico às vezes aceitando a atualização e às vezes permanecendo no mesmo lugar.
Para ilustrar esse processo, McElreath (2020) faz um analogia com a história de um Rei que tem que visitar as ilhas de seu reino permanecendo em cada uma por um tempo proporcional à população e permancendo o menos possível no mar. O Rei segue uma regra simples: se mover para uma ilha ao lado se sua população for maior do que a da ilha atual ou, sendo menor, mover para essa ilha se a proporção da população destino sobre a população atual for maior que um valor sorteado de uma distribuição uniforme entre 0 e 1. Deste modo o rei viajaria por todas as ilhas e ficaria um tempo proporcional à população em cada uma. As ilhas seriam o equivalente aos parâmetros, a população equivalente às probabilidades de cada parâmetro e o tempo de estadia equivalente às amostras de cada parâmetros.
McElreath sugere o seguinte algoritmo para simular a jornada do Rei
num_weeks <- 1e5
positions <- rep(0, num_weeks)
current <- 10
for(i in 1:num_weeks){
positions[i] <- current
proposal <- current + sample(c(-1,1), size =1)
if(proposal < 1) proposal <- 10
if(proposal > 10) proposal <- 1
prob_move <- proposal/current
current <- ifelse(runif(1) < prob_move, proposal, current)
}
plot(1:100, positions[1:100]) #impossível dizer com antecedência onde o Rei vai estar
plot(table(positions)) # mas ele terá explorado todas as ilhas ficando em cada uma um tmpo proporcional ao tamanho.
Esse algoritmo é o ancestral direto de outros algoritmos da famíla MCMC. Um dos mais conhecidos e utilizados é o Gibbs Sampling onde a distribuição do valor de parâmetro proposto é ela mesma atualizada. O Gibbs sampling estima cada parâmetro por vez como uma ‘retirada’ de uma distribuição condicionando nos dados e nos valores dos outros parâmetros. Um exemlo está no cap. 18 de Gelman e Hill
[Exemplo Gibbs Sampler]
Softwares como o BUGS (Bayesian inference UsingGibbs Sampling) ou JAGS (Just Another Gibbs Sampler) implementam o Gibbs Sampler e é o utilizado por Gelman e Hill para ajustar os modelos multinível do livro de 2007.
Recentemente Gelman e outros parceiros desenvolveram o Stan (implementado no R pelo pacote rstan) que usa outra estratégia de estimação conhecida como Hamiltonian Monte Carlo. Neste método imagina-se que o vetor de parâmetros nos dá a posição de uma partícula no log da distribuição posterior (que como vimos acima pode parecer uma colina). Damos então um impulso nesta partícula, como um peteleco em uma bolinha de gude, e assim que a partícula perde momento e para, registramos o valor do parâmetro e recomeçamos (ver pg. 275 do McElreath). Esse método permite superar um problema com o Gibbs Sampler ou o Metrópolis. Como esses algoritmos exploram aos poucos a vizinhança eles podem ficar presos em certos platôs nos dados, zonas de alta correlação entre os parâmetros. Ao dar saltos o HMC evita esse problema.