Cristian Villegas: Doutor em Estatística pela USP (2010). Atualmente é professor doutor na ESALQ-USP e tem experiência na área de modelos simétricos e semi-paramétricos.
clobos@usp.br
Eduardo E. Ribeiro Jr: Bacharel em Estatística pela UFPR (2016). Atualmente é mestrando no PPG em Estatística e Experimentação na ESALQ-USP e tem experiência na área de modelos lineares generalizados e extensões.
jreduardo@usp.br
Roseli A. Leandro: Doutora em Estatística e Experimentação pela ESALQ-USP (1997). Atualmente é professora doutora na ESALQ-USP e tem experiência na área de inferência bayesiana.
raleandr@usp.br
Desenvolvimento e códigos-fonte (dos apps e dos slides):
https://github.com/jreduardo/minicurso-seeb
Webpage com os slides e os materiais:
http://www.leg.ufpr.br/~eduardojr/courses/seeb
Pacote R com ferramentas para análise de cadeias MCMC:
http://github.com/jreduardo/mcmctools
Aumento da utilização da inferência bayesiana para análise de dados por meio de modelos mais fidedignos a realidade;
O melhor uso dos métodos disponíveis para inferência bayesiana dependem de um bom entendimento básico da área;
Conceitos importantes de inferência bayesiana podem ser melhor ilustrados por meio de recursos interativos.
Parte I: Introdução ao Shiny;
Parte II: Conceitos de inferência bayesiana com Shiny;
Parte III: Usando Shiny para diagnóstico de cadeias MCMC.
Exemplos:
Referências originais
Motivações
Primeira publicação em periódico científico
Wikipedia: An Essay towards solving a Problem in the Doctrine of Chances Philosophical Transactions − the world's first science journal
Teorema de Bayes (teorema da probabilidade inversa)
Se \(H\) denota uma hipótese e \(D\) denota os dados o teorema estabelece que:
\[ P(H \mid D) = \dfrac{P(D \mid H) \times P(H)}{ P(D)}, \]
O teorema de Bayes mostra como probabilidades mudam à luz de evidências: dados.
Tendo especificado \(P(D \mid H)\) e \(P(H)\) o teorema fornece uma solução do problema de como aprender através dos dados.
Teorema de Bayes, no contexto de inferência estatística
\[ \pi(\boldsymbol{\theta} \mid \boldsymbol{y}) = \dfrac{\mathcal{L}(\boldsymbol{\theta} \mid \boldsymbol{y}) \pi(\boldsymbol{\theta})}{ \int \mathcal{L}(\boldsymbol{\theta} \mid \boldsymbol{y}) \pi(\boldsymbol{\theta})d\boldsymbol{\theta}} \propto \mathcal{L}(\boldsymbol{\theta} \mid \boldsymbol{y}) \pi(\boldsymbol{\theta}), \]
em que,
A distribuição a priori, \(\pi(\boldsymbol{\theta})\) que representa o que é conhecido a respeito dos parâmetros desconhecidos antes que os dados sejam observados, desempenha um importante papel na inferência bayesiana.
\(\pi(\boldsymbol{\theta})\) pode ser utilizada para representar conhecimento prévio ou ignorância relativa.
Do ponto de vista bayesiano, toda inferência a respeito do parâmetro está baseada na distribuição a posteriori, ou seja,
Dada a verossimilhança \(\mathcal{L}(\boldsymbol{\theta} \mid \boldsymbol{y})\) e a distribuição a priori \(\pi(\boldsymbol{\theta})\), o início de qualquer inferência paramétrica ou decisão sobre \(\boldsymbol{\theta}\) está na densidade a posteriori:
Sejam \(\boldsymbol{y} = \{y_{1}, y_{2}, \ldots, y_{n}\}\) uma amostra aleatória de uma distribuição binomial com parâmetro \(\theta\),
\[ f(y_i \mid \theta) = {m \choose y_i} \theta^{y_i} (1-\theta)^{m-y_i}, \quad y_i \in \{ 0, 1, 2, \ldots, m\}; \]
\[ \begin{aligned} \mathcal{L}(\theta \mid \boldsymbol{y}) &= \prod_{i=1}^{n} {m \choose y_i} \theta^{y_i} (1-\theta)^{n-y_i} \\ &= \theta^{\sum_i^n y_i} (1-\theta)^{n-\sum_i^n y_i} \prod_{i=1}^n {m \choose y_i}, \quad \theta \in [0,1]. \end{aligned} \]
\[ \pi(\theta) = \dfrac{\Gamma(\alpha+\beta)}{\Gamma(\alpha) \Gamma(\beta)} \theta^{\alpha-1}(1-\theta)^{\beta-1}, \alpha > 0, \beta > 0. \]
\[ \begin{gathered} \pi(\theta \mid \boldsymbol{y}) \propto \mathcal{L}(\theta \mid \boldsymbol{y}) \pi(\theta) = \theta^{\sum_i^{n} y_i + \alpha - 1} (1 - \theta)^{nm-\sum_i^n y_i + \beta -1} \\ \theta \mid \boldsymbol{y} \sim \text{Beta} \left( \sum_i^{n} y_i + \alpha,\; nm - \sum_{i=1}^n y_i + \beta \right) \end{gathered} \]
\[ \begin{aligned} f(\tilde{y} \mid \boldsymbol{y}) &= \int_{\mathcal{R}} f(\tilde{y} \mid \theta) \pi(\theta \mid \boldsymbol{y}) d\boldsymbol{\theta} =\\ &= \int_0^1 {\tilde{m} \choose \tilde{y}} \theta^{\tilde{y}} (1-\theta)^{ \tilde{m}-\tilde{y} } \theta^{\sum_i^{n} y_i + \alpha -1}(1 - \theta)^{nm-\sum_i^n y_i + \beta -1} d \theta,\\ &\tilde{y} \in \{ 0, 1, 2, \ldots, \tilde{m}\}. \end{aligned} \]
e portanto, \[ \tilde{y} \mid \boldsymbol{y} \sim \text{Beta-Binomial}( \tilde{m},\; \sum_{i=1}^n y_i + \alpha,\; n\tilde{m}-\sum_i^n y_i + \beta). \]
Sejam \(\boldsymbol{y}=\{y_{1}, y_{2}, \ldots, y_{n}\}\) uma amostra aleatória de uma distribuição de Poisson com parâmetro \(\theta\),
\[ f(y_{i} \mid \theta) = \dfrac{e^{-\theta} \theta^{y_i}}{y_i!}, \quad y_{i}= 0, 1, 2, \ldots; \]
\[ \mathcal{L}(\theta \mid \boldsymbol{y}) = \prod_{i=1}^{n} \dfrac{e^{-\theta} \theta^{y_i}}{y_i!} = \dfrac{e^{-n \theta} \theta^{\sum_{i=1}^{n} y_i}}{ \prod_{i=1}^{n}y_{i}!},\quad \theta >0. \]
\[ \pi(\theta) = \dfrac{\beta^\alpha}{ \Gamma(\alpha)}\theta^{\alpha-1} e^{-\beta \theta}. \]
\[ \begin{aligned} \pi(\theta \mid \boldsymbol{y}) &= \dfrac{ \mathcal{L}(\theta \mid \boldsymbol{y}) \pi(\theta)}{ \int_{0}^{\infty}\pi(\theta) \mathcal{L}(\theta \mid \boldsymbol{y}) \pi(\theta) d \theta },\\ &= \dfrac{(n+\beta)^{\alpha + n \bar{y} } \theta^{\alpha+n \bar{y}-1} e^{-(n+\beta) \theta}}{\Gamma(\alpha + n \bar{y} )},\\ \theta \mid \boldsymbol{y}&\; \sim \text{Gamma}(\alpha+n \bar{y},\; n + \beta) \end{aligned} \]
Observe que, desprezando-se a constante normalizadora, o kernel da distribuição a posteriori é dado por:
\[ \theta^{\alpha+n \bar{y}-1} ~ e^{-(n+\beta), \theta} \]
que representa o kernel de uma distribuição gama com parâmetro de forma, \(\alpha+n \bar{y}\) e parâmetro de taxa , \(n + \beta\).
Ou seja, desprezando-se a constante normalizadora, chega-se a mesma conclusão.
\[ \begin{aligned} f(\tilde{y} \mid \boldsymbol{y}) &= \int_{\mathcal{R}} f(\tilde{y} \mid \theta) \pi(\theta \mid \boldsymbol{y}) d\boldsymbol{\theta} =\\ &= \int_{0}^{\infty} \dfrac{e^{-\theta}\theta^{\tilde{y}}}{\tilde{y}!} \dfrac{\tilde{\beta}^{\tilde{\alpha}}}{\Gamma(\tilde{\alpha})} \theta^{\tilde{\alpha}-1} e^{-\tilde{\beta} \theta} d\theta\\ &= \dfrac{\tilde{\beta}^{\tilde{\alpha}}}{\Gamma(\tilde{\alpha}) \tilde{y}!} \dfrac{\Gamma(\tilde{\alpha}+\tilde{y}) \Gamma(\tilde{\beta}+1)}{\Gamma(\tilde{\alpha}+\tilde{y} + \tilde{\beta}+1)} \end{aligned} \]
considerando-se \(\tilde{\alpha}=\alpha+n \bar{y}\) e \(\tilde{\beta} = n+\beta\).
Quando a distribuição a posteriori possuir uma expressão complexa, a obtenção de resumos a posteriori pode não ser simples;
Uma forma bastante utilizada na inferência bayesiana para superar essa dificuldade é gerar uma cadeia de Markov cuja distribuição de equilíbrio seja a distribuição de interesse, e assim os resumos a posteriori serão calculados por meio da amostra obtida (sampling based)
Um algoritmo bastante utilizado é o algoritmo de Metropolis-Hastings.
Seja \(\pi(\theta \mid y)\) a distribuição alvo (target distribution) e \(\theta^{(t)}\) o valor amostral atual de \(\pi(\theta\mid y)\). O algoritmo Metropolis-Hastings desenvolve-se da seguinte maneira:
\[ r = r(\theta ^{t}, \theta^{*}) = \dfrac{\pi(\theta^{*}\mid y) \;/\; g(\theta^{*} \mid \theta^{t})}{ \pi(\theta^{t} \mid y) \;/\; g(\theta^{t} \mid \theta^{*}) } \]
Seja \(\pi(\theta \mid y)\) a distribuição alvo (target distribution) e \(\theta^{(t)}\) o valor amostral atual de \(\pi(\theta\mid y)\). O algoritmo Metropolis-Hastings desenvolve-se da seguinte maneira:
\[ r = r(\theta ^{t}, \theta^{*}) = \dfrac{\pi(\theta^{*}\mid y) \;/\; g(\theta^{*})}{ \pi(\theta^{t} \mid y) \;/\; g(\theta^{t}) } \]
Seja \(\pi(\theta \mid y) = k q(\theta \mid y)\) a distribuição alvo (target distribution) e \(\theta^{(t)}\) o valor amostral atual de \(\pi(\theta\mid y)\). O algoritmo Metropolis-Hastings desenvolve-se da seguinte maneira:
\[ r = r(\theta ^{t}, \theta^{*}) = \dfrac{q(\theta^{*}\mid y) \;/\; g(\theta^{*}\mid \theta^{t})}{ q(\theta^{t} \mid y) \;/\; g(\theta^{t}\mid \theta^{*}) } \]
Seja \(\pi(\theta \mid y) = k q(\theta \mid y)\) a distribuição alvo (target distribution) e \(\theta^{(t)}\) o valor amostral atual de \(\pi(\theta\mid y)\). O algoritmo Metropolis-Hastings desenvolve-se da seguinte maneira:
\[ r = r(\theta ^{t}, \theta^{*}) = \dfrac{q(\theta^{*}\mid y) \;/\; g(\theta^{*})}{ q(\theta^{t} \mid y) \;/\; g(\theta^{t}) } \]
Ferramentas utéis para manipulação de objetos MCMC
# install.packages("devtools")
install_github("jreduardo/mcmctools")
help(package = "mcmctools", help_type = "html")
Funções:
mcmctools::mct_select
mcmctools::mct_burnin
mcmctools::mct_thin
mcmctools::mct_merge
mcmctools::mct_shinyapp
#-------------------------------------------
# (toy) Data set
str(mtcars)
help(mtcars, help_type = "html")
splom(~mtcars, data = mtcars)
#-------------------------------------------
# Define linear model in JAGS
str_model <- "
model {
for (i in 1:n) {
mu[i] <- X[i, ] %*% beta[]
y[i] ~ dnorm(mu[i], tau)
}
for (j in 1:p) {
beta[j] ~ dnorm(0, 0.001)
}
tau ~ dgamma(0.001, 0.001)
}
"
#-------------------------------------------
# Initialize the JAGS model
X <- model.matrix(
~factor(cyl) + factor(am) + hp + drat + qsec,
data = mtcars)
jags <- jags.model(
file = textConnection(str_model),
data = list(
n = nrow(X),
p = ncol(X),
X = X,
y = mtcars$mpg),
n.chains = 3)
#-------------------------------------------
# Posterior sampling
samples <- coda.samples(
model = jags,
variable.names = c("beta", "tau", "mu"),
n.iter = 5000)
str(samples)
#-------------------------------------------
mct_shinyapp(samples)