⏱️ 0:00
Docente: Fabio Cop (fcferreira@unifesp.br) | Curso: Bacharelado Interdisciplinar em Ciências do Mar - Bict Mar | UC: Probabilidade e Estatística - Aula 05

Distribuição a Posteriori no Modelo Normal

Verossimilhança, teorema de Bayes e estimação com Stan

O que construímos na Aula 04

✅ Distribuição Normal como modelo probabilístico para variáveis contínuas

✅ O modelo generativo de três níveis para alturas adultas:

$$Y \sim \text{Normal}(\mu,\, \sigma)$$ $$\mu \sim \text{Normal}(160,\, 20) \qquad \sigma \sim \text{Exponencial}(0{,}1)$$

✅ Distribuição preditiva a priori: verificação de que as distribuições a priori geram alturas biologicamente razoáveis

✅ O modelo ainda não viu nenhum dado real

O que fazemos hoje

Teoria (40 min)
  1. Verossimilhança: medida de compatibilidade entre parâmetros e dados
  2. Teorema de Bayes: combinando priori e verossimilhança
  3. Resumos da distribuição a posteriori e intervalo de credibilidade
  4. Introdução ao Stan e ao rstan
Prática (60 min)
  1. Exploração das alturas adultas e especificação das distribuições a priori
  2. Ajuste com brms: checagem preditiva a priori, distribuição a posteriori e checagem preditiva da posteriori

Dataset: Howell1 (McElreath, 2020)

Verossimilhança e o Teorema de Bayes

Quantificando a compatibilidade entre parâmetros e dados

Dados reais: alturas adultas !Kung San

Dataset Howell1 (Dados: McElreath, 2020)
  • 544 indivíduos !Kung San do deserto de Kalahari
  • Adultos ($\text{age} \geq 18$): $n = 352$ observações
  • Variável de interesse: height (cm)
  • Histograma: forma aproximadamente simétrica em torno de $\bar{y} \approx 154{,}6$ cm

O modelo especificado na Aula 04 assume $Y \sim \text{Normal}(\mu, \sigma)$.

A pergunta agora: quais valores de $\mu$ e $\sigma$ são mais compatíveis com esses 352 dados observados?

A verossimilhança

Definição: a verossimilhança $\mathcal{L}(\mu, \sigma; \mathbf{y})$ é o produto das densidades da Distribuição Normal avaliadas em cada observação, tratado como função dos parâmetros para dados fixos. $$\mathcal{L}(\mu, \sigma; \mathbf{y}) = \prod_{i=1}^n \frac{1}{\sigma\sqrt{2\pi}}\exp\!\left(-\frac{(y_i - \mu)^2}{2\sigma^2}\right)$$

Pergunta que a verossimilhança responde: "Se os parâmetros fossem $(\mu, \sigma)$, qual seria a probabilidade de observar exatamente estes dados?"

  • Pares $(\mu, \sigma)$ que atribuem alta densidade às observações: verossimilhança alta
  • Pares que tornam os dados improváveis: verossimilhança baixa

Log-verossimilhança

Por que usar o logaritmo? O produto de $n$ densidades produz números muito pequenos para $n$ grande. O logaritmo transforma o produto em soma: $$\log \mathcal{L}(\mu, \sigma; \mathbf{y}) = \sum_{i=1}^n \log f(y_i \mid \mu, \sigma)$$

Em R, usando dnorm() com log = TRUE:

alturas <- c(150, 154, 158, 162, 165)
sum(dnorm(alturas, mean = 155, sd = 8, log = TRUE))
# resultado: sempre negativo para variáveis contínuas
Valores menos negativos (mais próximos de zero) indicam maior compatibilidade dos parâmetros com os dados.

🗳️ Enquete

Dois modelos são ajustados com os mesmos dados:

  • Modelo A: log-verossimilhança $= -450$
  • Modelo B: log-verossimilhança $= -420$

Qual modelo é mais compatível com os dados?

(A) Modelo A, porque $-450 < -420$

(B) Modelo B, porque log-verossimilhança mais próxima de zero indica maior compatibilidade

(C) Ambos são igualmente compatíveis

(D) Depende do número de parâmetros de cada modelo

O teorema de Bayes

Combinando distribuição a priori e verossimilhança: $$P(\mu, \sigma \mid \mathbf{y}) \propto \mathcal{L}(\mathbf{y}; \mu, \sigma) \times P(\mu) \times P(\sigma)$$
Distribuição a posteriori
$P(\mu, \sigma \mid \mathbf{y})$
O que sabemos sobre os parâmetros após os dados
Verossimilhança
$\mathcal{L}(\mathbf{y}; \mu, \sigma)$
Compatibilidade dos parâmetros com os dados
Distribuições a priori
$P(\mu) \times P(\sigma)$
O que sabemos antes dos dados

O símbolo $\propto$ ("proporcional a") indica que o produto ainda precisa ser normalizado para integrar 1.

A superfície a posteriori

A distribuição a posteriori $P(\mu, \sigma \mid \mathbf{y})$ é uma superfície bidimensional: para cada par $(\mu, \sigma)$, existe uma densidade proporcional ao produto da verossimilhança pelas distribuições a priori.
Com $n = 352$ adultos (Howell1):
  • Verossimilhança concentrada em $\mu \approx 154{,}6$ cm
  • Verossimilhança concentrada em $\sigma \approx 7{,}7$ cm
  • Dados pesam muito mais do que as distribuições a priori
Distribuição a posteriori:
  • Combina verossimilhança com distribuições a priori
  • Com $n$ grande, distribui-se perto da máxima verossimilhança
  • Com $n$ pequeno, as distribuições a priori têm influência perceptível

Aproximação por Grade

Calculando a distribuição a posteriori célula a célula

A lógica da grade

Cada célula da grade é uma hipótese sobre os parâmetros $(\mu, \sigma)$. A grade calcula quão bem cada hipótese explica os dados (verossimilhança) e pondera pela crença prévia (distribuição a priori).

Os quatro passos da aproximação por grade:
  1. Definir uma grade de valores candidatos para $(\mu, \sigma)$
  2. Calcular a log-verossimilhança em cada ponto da grade
  3. Somar os logaritmos das distribuições a priori à log-verossimilhança
  4. Exponenciar e normalizar para obter probabilidades relativas

Construindo a grade em R

dados <- read.csv("https://raw.githubusercontent.com/FCopf/prob-est-2026/refs/heads/main/_bibliography/mcelreath-2019/mcelreath-2019-data-sets/Howell1.csv",
                 sep = ";")
alturas <- dados$height[dados$age >= 18]  # n = 352 adultos

# Grade bidimensional de parâmetros
mu_grade    <- seq(150, 160, length.out = 200)
sigma_grade <- seq(4, 12,   length.out = 200)
grade       <- expand.grid(mu = mu_grade, sigma = sigma_grade)

# Log-verossimilhança em cada célula
grade$log_lik <- sapply(seq_len(nrow(grade)), function(i) {
  sum(dnorm(alturas, mean = grade$mu[i], sd = grade$sigma[i], log = TRUE))
})

# Log-distribuições a priori
grade$log_prior_mu    <- dnorm(grade$mu,   mean = 160, sd = 20, log = TRUE)
grade$log_prior_sigma <- dexp(grade$sigma, rate = 0.1,           log = TRUE)

# Distribuição a posteriori: normalizada
grade$log_post  <- grade$log_lik + grade$log_prior_mu + grade$log_prior_sigma
grade$prob_post <- exp(grade$log_post - max(grade$log_post))
grade$prob_post <- grade$prob_post / sum(grade$prob_post)

Visualizando a superfície a posteriori

Superfície conjunta $P(\mu, \sigma \mid \mathbf{y})$
  • Eixo $x$: valores candidatos de $\mu$ (cm)
  • Eixo $y$: valores candidatos de $\sigma$ (cm)
  • Cor: probabilidade relativa de cada par
  • Pico: estimativa MAP
Distribuições marginais
  • $P(\mu \mid \mathbf{y})$: projeção da superfície sobre o eixo $\mu$ (soma sobre $\sigma$)
  • $P(\sigma \mid \mathbf{y})$: projeção da superfície sobre o eixo $\sigma$ (soma sobre $\mu$)
A distribuição a posteriori está concentrada em uma região estreita do espaço de parâmetros. Com $n = 352$, os dados são muito informativos sobre $\mu$ e $\sigma$.

Amostrando da distribuição a posteriori

A grade produz uma tabela de probabilidades discretas. Para calcular resumos, amostra-se pares $(\mu, \sigma)$ com probabilidade proporcional à distribuição a posteriori:

set.seed(2026)
idx_amostras   <- sample(nrow(grade), size = 6000, replace = TRUE,
                         prob = grade$prob_post)
amostras_grade <- grade[idx_amostras, c("mu", "sigma")]

# Resumos da distribuição a posteriori de μ
mean(amostras_grade$mu)              # média a posteriori
median(amostras_grade$mu)            # mediana a posteriori
quantile(amostras_grade$mu, c(0.055, 0.945))  # intervalo de credibilidade 89%
As distribuições marginais de $\mu$ e $\sigma$ são extraídas simplesmente ignorando as amostras do outro parâmetro.

Resumos e Introdução ao Stan

Quantificando a incerteza e escalando para modelos complexos

Resumos da distribuição a posteriori

Resumo Definição Código R
Estimativa MAP Par $(\hat\mu, \hat\sigma)$ que maximiza $P(\mu, \sigma \mid \mathbf{y})$ grade[which.max(grade$prob_post), ]
Média a posteriori Média das amostras da distribuição a posteriori mean(amostras_grade$mu)
Mediana a posteriori Mediana das amostras da distribuição a posteriori median(amostras_grade$mu)
Intervalo de credibilidade 89% Intervalo que contém 89% da massa da distribuição a posteriori quantile(amostras_grade$mu, c(0.055, 0.945))

O intervalo de credibilidade de 89%

Interpretação direta: dado o modelo e os dados observados, há 89% de probabilidade de que o parâmetro esteja nesse intervalo.

Esta interpretação não está disponível para o intervalo de confiança frequentista.

Por que 89% e não 95%?
  • 95% carrega conotação frequentista (limiar de significância)
  • 89% é primo: lembra que a escolha do nível é sempre uma decisão do pesquisador
  • Com amostras de tamanho moderado, 89% é mais estável numericamente do que 95%

Pergunta para reflexão

A distribuição a posteriori de $\mu$ está centrada em 154,6 cm com desvio padrão de 0,4 cm.

O que esse desvio padrão de 0,4 cm representa?

(A) A variabilidade típica das alturas individuais na população.

(B) A incerteza residual sobre o valor da média populacional $\mu$, dado o modelo e os dados.

(C) O erro de medição das alturas no dataset.

Limitação da grade: a maldição da dimensionalidade

A aproximação por grade avalia o produto verossimilhança × distribuições a priori em cada célula da grade.

O problema com mais parâmetros:
  • 2 parâmetros, 200 valores cada: $200^2 = 40.000$ avaliações
  • 5 parâmetros, 100 valores cada: $100^5 = 10^{10}$ avaliações
  • 10 parâmetros, 100 valores cada: $100^{10} = 10^{20}$ avaliações
A solução: Stan. O HMC amostra diretamente da distribuição a posteriori, concentrando amostras nas regiões de alta probabilidade sem avaliar a grade completa.

O programa Stan

Estrutura obrigatória do programa Stan:

data {
  int<lower=0> N;   // número de observações
  vector[N] y;       // alturas observadas
}
parameters {
  real mu;           // média populacional
  real<lower=0> sigma; // desvio padrão
}
model {
  // Distribuições a priori
  mu    ~ normal(160, 20);
  sigma ~ exponential(0.1);
  // Verossimilhança
  y     ~ normal(mu, sigma);
}
bloco de dados: declara os dados que chegam do R
bloco de parâmetros: declara os parâmetros a estimar. <lower=0> restringe $\sigma > 0$
bloco do modelo: distribui­ções a priori + verossimilhança. Exatamente as equações da Aula 04

Ajustando o modelo com rstan

library(rstan)
library(posterior)

# Preparar dados no formato esperado pelo Stan
dados_stan <- list(N = length(alturas), y = alturas)

# Ajustar o modelo (model_code: string com o programa Stan)
fit <- stan(model_code = model_code,
            data   = dados_stan,
            chains = 4,
            iter   = 1000,
            warmup = 500,
            refresh = 0)

# Extrair amostras da distribuição a posteriori
amostras_stan <- as_draws_df(fit)

# Resumos
summarise_draws(amostras_stan, mean, median, sd,
                ~quantile(.x, c(0.055, 0.945)))

O rstan compila o programa Stan em C++, amostra com HMC e retorna um objeto stanfit com as amostras de todas as cadeias.

Conexão: teoria e código

Conceito matemático Código R (grade) Código Stan
$\sum \log P(y_i \mid \mu, \sigma)$ sum(dnorm(..., log=TRUE)) y ~ normal(mu, sigma)
$\log P(\mu) = \log \text{Normal}(\mu; 160, 20)$ dnorm(mu, 160, 20, log=TRUE) mu ~ normal(160, 20)
$\log P(\sigma) = \log \text{Exp}(\sigma; 0{,}1)$ dexp(sigma, rate=0.1, log=TRUE) sigma ~ exponential(0.1)
Amostras da distribuição a posteriori sample(..., prob=prob_post) as_draws_df(fit)

Na grade, calculamos a distribuição a posteriori célula a célula. No Stan, o HMC amostra diretamente da distribuição a posteriori. O resultado, em ambos os casos, é um conjunto de amostras de $(\mu, \sigma)$.

Distribuição Preditiva da Posteriori

Predizendo novas observações com incerteza paramétrica

Parâmetros vs. novas observações

Distribuição a posteriori
$P(\mu, \sigma \mid \mathbf{y})$

Descreve a incerteza sobre os parâmetros após observar os dados.

Pergunta: "Onde está $\mu$?"

Distribuição preditiva da posteriori
$P(\tilde{y} \mid \mathbf{y})$

Descreve a distribuição esperada de novas observações, incorporando a incerteza sobre os parâmetros.

Pergunta: "Qual a altura esperada de um novo adulto?"

Formalmente: $P(\tilde{y} \mid \mathbf{y}) = \int\!\int P(\tilde{y} \mid \mu, \sigma)\, P(\mu, \sigma \mid \mathbf{y})\, d\mu\, d\sigma$

Gerando a distribuição preditiva da posteriori

Três passos, repetidos para cada amostra $s$ da distribuição a posteriori:
  1. Extrair $(\mu^{(s)}, \sigma^{(s)})$ da distribuição a posteriori
  2. Gerar $\tilde{y}^{(s)} \sim \text{Normal}(\mu^{(s)}, \sigma^{(s)})$
  3. O conjunto $\{\tilde{y}^{(1)}, \tilde{y}^{(2)}, \ldots\}$ é uma amostra da distribuição preditiva da posteriori
# Distribuição preditiva da posteriori (usando amostras do Stan)
y_post <- rnorm(nrow(amostras_stan),
                mean = amostras_stan$mu,
                sd   = amostras_stan$sigma)

Comparando as três distribuições

Após o ajuste do modelo, três distribuições de alturas podem ser comparadas:

Preditiva a priori
$P(Y)$

Alturas simuladas antes dos dados. Baseada em $\text{Normal}(160, 20)$ e $\text{Exp}(0{,}1)$. Mais larga.
Dados observados
$\mathbf{y}$

$n = 30$ alturas reais do Howell1. Referência empírica.
Preditiva da posteriori
$P(\tilde{y} \mid \mathbf{y})$

Alturas preditas após o ajuste. Centrada nos dados, mais estreita do que a preditiva a priori.

A distribuição preditiva da posteriori é mais larga do que $\text{Normal}(\hat\mu, \hat\sigma)$, pois incorpora a incerteza residual sobre os parâmetros.

Síntese

O que aprendemos hoje

Conceito Definição Ferramenta em R
Verossimilhança Compatibilidade de $(\mu, \sigma)$ com os dados observados sum(dnorm(..., log=TRUE))
Distribuição a posteriori $P(\mu, \sigma \mid \mathbf{y}) \propto \mathcal{L} \times P(\mu) \times P(\sigma)$ grade$prob_post / as_draws_df(fit)
Estimativa MAP Par $(\hat\mu, \hat\sigma)$ que maximiza a distribuição a posteriori which.max(grade$prob_post)
Intervalo de credibilidade 89% Intervalo com 89% da massa da distribuição a posteriori quantile(amostras, c(0.055, 0.945))
Stan / rstan HMC como alternativa escalável à grade stan(...) + as_draws_df()
Distribuição preditiva da posteriori Distribuição de novas observações incorporando incerteza paramétrica rnorm(n, amostras$mu, amostras$sigma)

O fluxo bayesiano completo

Cinco etapas percorridas nas Aulas 01 a 05:
  1. Aulas 01–03: hipóteses contrastantes, atualização de crenças, modelos generativos
  2. Aula 04: modelo probabilístico completo — Distribuição Normal, distribuições a priori, distribuição preditiva a priori
  3. Aula 05: dados reais, verossimilhança, distribuição a posteriori, resumos, Stan

Dúvidas?