O que é inferência aproximada (MCMC e Gibbs)?
Inferência aproximada por MCMC é uma técnica para amostrar de distribuições complexas sem fórmula fechada. MCMC significa Markov Chain Monte Carlo, que gera uma cadeia de amostras dependentes. A cadeia é construída para convergir para a distribuição alvo (posteriori). Gibbs sampling é um caso especial do MCMC onde cada variável é amostrada condicionalmente. A cada passo, atualiza-se uma variável por vez, usando as demais fixas. Isso é feito quando a distribuição condicional completa é conhecida. Diferentemente da inferência exata, MCMC é escalável para redes com muitas variáveis. Ele não sofre com explosão combinatória, mas produz amostras correlacionadas. Portanto, ele é amplamente usado em problemas de alta dimensão.
Características fundamentais do MCMC
O MCMC possui três características principais que o definem. Primeiro, ele é assintoticamente exato: converge para a verdadeira posteriori com muitas amostras. Segundo, as amostras iniciais (burn-in) são descartadas para eliminar o efeito do início. Terceiro, a autocorrelação entre amostras reduz a eficiência, exigindo thinning. A convergência é monitorada por diagnósticos como o fator de escala de Gelman-Rubin. O MCMC não requer conhecimento da constante de normalização da posteriori. Isso é uma grande vantagem sobre métodos analíticos.
Vantagens e aplicações típicas
A principal vantagem é a capacidade de lidar com modelos hierárquicos e não-lineares. Ele é usado em aprendizado de máquina, econometria e genética populacional. Também é aplicado em análise de imagens e processamento de sinais. Contudo, a convergência pode ser lenta e difícil de diagnosticar.
O algoritmo de Metropolis-Hastings é o MCMC mais geral. Ele propõe um novo estado a partir de uma distribuição de proposta. Aceita-se o novo estado com probabilidade dada pela razão de verossimilhanças. Gibbs sampling é mais eficiente quando as condicionais são fáceis de amostrar. Por exemplo, em modelos com prioris conjugados, as condicionais são distribuições conhecidas. A cadeia de Gibbs amostra cada variável sequencialmente, mantendo as outras fixas. Após um número suficiente de iterações, as amostras representam a posteriori. A média das amostras estima a esperança a posteriori. Intervalos de credibilidade são obtidos dos quantis empíricos das amostras. O MCMC pode ser implementado em bibliotecas como PyMC3 e Stan. Essas ferramentas automatizam a escolha do amostrador e diagnósticos. O usuário só precisa definir o modelo e os dados. A inferência aproximada é a espinha dorsal da estatística bayesiana moderna. Assim, o MCMC é uma ferramenta indispensável para problemas reais.
Um exemplo clássico é estimar a média e variância de uma distribuição normal com dados observados. A posteriori conjunta não tem forma simples, mas as condicionais são normais e gamma inversa. Gibbs sampling amostra média dado variância, e variância dado média, alternadamente. Após algumas milhares de iterações, as amostras convergem para a posteriori verdadeira.
Enunciado do exemplo clássico
Implemente Gibbs sampling para estimar os parâmetros (μ, σ²) de uma normal com dados sintéticos. Gerar 50 dados de N(μ=5, σ²=4). Use prioris: μ ~ N(0, 1000) e σ² ~ Inverse-Gamma(3, 2). As condicionais completas são: μ | dados, σ² ~ N( média ponderada, variância ) e σ² | dados, μ ~ Inverse-Gamma( α + n/2, β + 0.5*Σ(xᵢ-μ)² ). Execute 5000 iterações, descarte os primeiros 1000 (burn-in). Plote as cadeias e os histogramas das amostras para μ e σ².
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 |
import numpy as np import matplotlib.pyplot as plt from scipy.stats import norm, invgamma # Gerar dados sintéticos np.random.seed(42) n = 50 mu_real = 5.0 sigma2_real = 4.0 dados = np.random.normal(mu_real, np.sqrt(sigma2_real), n) # Prioris: # μ ~ N(0, 1000) mu_prior_mean = 0.0 mu_prior_var = 1000.0 # σ² ~ Inv-Gamma(alpha=3, beta=2) alpha_prior = 3.0 beta_prior = 2.0 # Gibbs sampling n_iter = 5000 burn_in = 1000 amostras_mu = np.zeros(n_iter) amostras_sigma2 = np.zeros(n_iter) # Inicialização mu_atual = 0.0 sigma2_atual = 1.0 for i in range(n_iter): # Amostrar μ | dados, σ² # Variância condicional var_mu = 1.0 / (1.0/mu_prior_var + n/sigma2_atual) # Média condicional mean_mu = var_mu * (mu_prior_mean/mu_prior_var + np.sum(dados)/sigma2_atual) mu_atual = np.random.normal(mean_mu, np.sqrt(var_mu)) # Amostrar σ² | dados, μ alpha_post = alpha_prior + n/2.0 beta_post = beta_prior + 0.5 * np.sum((dados - mu_atual)**2) sigma2_atual = invgamma.rvs(alpha_post, scale=beta_post) amostras_mu[i] = mu_atual amostras_sigma2[i] = sigma2_atual # Descartar burn-in mu_amostras = amostras_mu[burn_in:] sigma2_amostras = amostras_sigma2[burn_in:] # Estatísticas media_mu = np.mean(mu_amostras) ic_mu = np.percentile(mu_amostras, [2.5, 97.5]) media_sigma2 = np.mean(sigma2_amostras) ic_sigma2 = np.percentile(sigma2_amostras, [2.5, 97.5]) print(f"μ real = {mu_real}, estimado = {media_mu:.3f} (IC 95%: [{ic_mu[0]:.3f}, {ic_mu[1]:.3f}])") print(f"σ² real = {sigma2_real}, estimado = {media_sigma2:.3f} (IC 95%: [{ic_sigma2[0]:.3f}, {ic_sigma2[1]:.3f}])") # Gráficos plt.figure(figsize=(12, 10)) # Cadeias de μ plt.subplot(2, 2, 1) plt.plot(mu_amostras, 'b-', linewidth=0.5, alpha=0.7) plt.axhline(mu_real, color='red', linestyle='--', label=f'μ real = {mu_real}') plt.axhline(media_mu, color='green', linestyle='-', label=f'Média amostral = {media_mu:.3f}') plt.title('Cadeia de Gibbs - μ') plt.xlabel('Iteração (após burn-in)') plt.ylabel('μ') plt.legend() plt.grid(True) # Histograma de μ plt.subplot(2, 2, 2) plt.hist(mu_amostras, bins=40, density=True, color='blue', alpha=0.6, edgecolor='black') x = np.linspace(mu_real-2, mu_real+2, 200) # Aproximar posteriori por normal (teórica) plt.plot(x, norm.pdf(x, media_mu, np.std(mu_amostras)), 'r-', linewidth=2, label='Aprox. normal') plt.axvline(mu_real, color='red', linestyle='--', label=f'μ real') plt.title('Distribuição posteriori de μ') plt.xlabel('μ') plt.ylabel('Densidade') plt.legend() plt.grid(True) # Cadeias de σ² plt.subplot(2, 2, 3) plt.plot(sigma2_amostras, 'g-', linewidth=0.5, alpha=0.7) plt.axhline(sigma2_real, color='red', linestyle='--', label=f'σ² real = {sigma2_real}') plt.axhline(media_sigma2, color='darkgreen', linestyle='-', label=f'Média amostral = {media_sigma2:.3f}') plt.title('Cadeia de Gibbs - σ²') plt.xlabel('Iteração (após burn-in)') plt.ylabel('σ²') plt.legend() plt.grid(True) # Histograma de σ² plt.subplot(2, 2, 4) plt.hist(sigma2_amostras, bins=40, density=True, color='green', alpha=0.6, edgecolor='black') x2 = np.linspace(0, 10, 200) # Aproximação por Inv-Gamma (teórica) plt.plot(x2, invgamma.pdf(x2, alpha_prior + n/2, scale=beta_prior + 0.5*np.var(dados)*n), 'r-', linewidth=2, label='Aprox. Inv-Gamma') plt.axvline(sigma2_real, color='red', linestyle='--', label=f'σ² real') plt.title('Distribuição posteriori de σ²') plt.xlabel('σ²') plt.ylabel('Densidade') plt.legend() plt.grid(True) plt.tight_layout() plt.show() |
Este código implementa Gibbs sampling com prioris conjugados. As cadeias mostram a convergência rápida para os valores reais. Os histogramas exibem as distribuições posteriores aproximadas. Os intervalos de credibilidade contêm os valores verdadeiros. Para iniciantes, este exemplo demonstra o poder do MCMC prático. A inferência aproximada é, portanto, essencial para modelos complexos.