f_x <- function(x) ifelse(x < 0, 0, x*exp(-x))
x_grid <- seq(0, 15, by=0.01)Simulação de dados - MMC e MCMC
Universidade Federal do Amazonas
Introdução
O objetivo deste trabalho é gerar dados simulados utilizando o Método de Monte Carlo (MMC) e os Métodos de Monte Carlo via Cadeias de Markov (MCMC).
Inicialmente, será abordado o Método de Monte Carlo (MMC).
Questão 1
Utilizar o método de Monte Carlo para obter 1000 amostras de uma variável aleatória com densidade de probabilidade proporcional a
\[ f(x) \propto x e^{-x}, \quad x \ge 0. \] A partir dessas amostras, calcular numericamente a distribuição acumulada correspondente.
Visão Geral de MMC
O que é o Método de Monte Carlo (MMC)?
- O Método de Monte Carlo é uma técnica estatística que utiliza amostragem aleatória para aproximar valores esperados, integrais e probabilidades. Em Inferência Bayesiana, ele é usado para aproximar integrais do tipo:
\[\mathbb{E}[g(\theta)] = \int g(\theta), p(\theta \mid y), d\theta, \]onde \(p(\theta \mid y)\) é a distribuição a posteriori.
- Quando conseguimos gerar amostras independentes \(\theta^{(1)}, \ldots, \theta^{(N)} \sim p(\theta \mid y)\), o valor esperado pode ser estimado por:
\[ \widehat{\mathbb{E}}[g(\theta)] = \frac{1}{N}\sum_{i=1}^{N} g(\theta^{(i)}). \] - Esse é o princípio fundamental do Monte Carlo: substituir integrais analíticas por médias amostrais.
Resolução
Geração de amostras de \(f(x) = x e^{-x}\)
Passo 1: Definir a função alvo e calcular a CDF numericamente Explicação: Aqui definimos a função alvo \(f(x)\). A grade de valores x_grid será usada para calcular a CDF em vários pontos, numericamente.
Passo 2: Gerar amostras via transformação inversa Explicação:
integratecalcula a integral da função de 0 até cada ponto da grade \(\rightarrow\) isso nos dá a CDF acumulada.approxfuncria uma função interpoladora que representa a inversa da CDF, usada para gerar amostras.
cdf_values <- sapply(x_grid, function(val) integrate(f_x, 0, val)$value)
inv_cdf <- approxfun(cdf_values, x_grid)Passo 3: Gerar amostras uniformes e transformá-las Explicação:
runif(1000)gera 1000 números uniformes entre 0 e 1.inv_cdftransformamos essas amostras uniformes em amostras da nossa função alvo.- Isso é o princípio do método de Monte Carlo via inversão da CDF.
set.seed(123)
u <- runif(1000)
monte_carlo_samples <- inv_cdf(u)Histograma e Curva Teórica
hist(monte_carlo_samples,
freq = FALSE,
breaks = 30,
main = "Histograma da Simulação de Monte Carlo",
xlab = "Valores de x",
ylab = "Densidade",
col = "lightblue",
xlim = c(0, 11))
curve(f_x,
from = 0,
to = 11,
add = TRUE,
col = "red",
lwd = 2)
legend("topright",
legend = c("Amostras Simuladas", "Curva Teórica"),
col = c("lightblue", "red"),
lwd = c(10, 2),
bty = "n")
Análise do Gráfico e do Método MMC
Análise Estatística:
- O histograma das amostras geradas pelo Método de Monte Carlo mostra a distribuição dos valores simulados, que segue o formato da função alvo \(f(x) \propto x e ^{- x}\).
- O uso da transformação inversa da CDF permitiu gerar amostras que respeitam a probabilidade relativa de cada valor de \(x\).
- O gráfico ficou bem próximo da curva teórica, validando o método.
Questão 2
Simule por MCMC a distribuição de probabilidades \[ f(x) \propto x e^{-x}, \quad x \ge 0. \]
Faça simulações com 10000 iterações. Use uma função de aceitação gaussiana com \(\sigma= 0.1, 1, 10\). Faça figuras relevantes e comente como a escolha de \(\sigma\) afeta o algoritmo. Comente sobre o burn-in nesses casos.
Visão Geral do MCMC
A ideia central dos métodos de Monte Carlo via Cadeias de Markov (MCMC) é a seguinte:
- Construir uma cadeia de Markov cuja distribuição estacionária seja a distribuição alvo \(f\).
- Rodar a cadeia por um período longo (muitas iterações) até que a cadeia converja para essa distribuição estacionária.
- Após a convergência (e um período de burn-in), as amostras obtidas podem ser usadas para estimar quantidades de interesse da distribuição \(f\).
De forma geral: O uso do MCMC é para gerar amostras de uma distribuição de probabilidade complexa da qual não conseguiríamos amostrar de outra forma. Nós aceitamos a dependência (autocorrelação) como uma consequência do método. Depois, nós lidamos com ela através de técnicas como:
Análise de ACF: Para medir o nível de correlação.
Thinning: Para selecionar amostras com uma distância grande entre si, criando um conjunto final que é praticamente independente.
Resolução
Simulação por MCMC: Código Base
set.seed(132)
# Função alvo
func_alvo <- function(x) {
if (x < 0) return(0)
return(x * exp(-x))
}Vamos utilizar o método de Metropolis-Hastings para a resolução.
# Algoritmo de Metropolis-Hastings
metropolis_hastings <- function(sigma, niter, xini) {
cadeia <- numeric(niter)
x_atual <- xini
n_aceito <- 0
for (i in 1:niter) {
# 1. Propor novo estado
x_prop <- rnorm(1, mean = x_atual, sd = sigma)
# 2. Calcular razão de aceitação
if (x_prop < 0) {
alpha <- 0
} else {
alpha <- min(1, func_alvo(x_prop) / func_alvo(x_atual))
}
# 3. Aceitar ou rejeitar
if (runif(1) < alpha) {
x_atual <- x_prop
n_aceito <- n_aceito + 1
}
cadeia[i] <- x_atual
}
taxa_aceite <- n_aceito / niter
return(list(cadeia = cadeia, taxa_aceite = taxa_aceite))
}Parâmetros e Execução
# Parâmetros
niter <- 10000
xini <- 24
burn_in <- 7000 # Descartar as primeiras 7000 iterações
sigmas <- c(0.1, 1, 10)
# Rodar para cada sigma
set.seed(42) # Para reprodutibilidade
result1 <- metropolis_hastings(sigma = sigmas[1], niter = niter, xini)
result2 <- metropolis_hastings(sigma = sigmas[2], niter = niter, xini)
result3 <- metropolis_hastings(sigma = sigmas[3], niter = niter, xini)
# Descartar burn-in
cadeia_pos_burn_in1 <- result1$cadeia[-(1:burn_in)]
cadeia_pos_burn_in2 <- result2$cadeia[-(1:burn_in)]
cadeia_pos_burn_in3 <- result3$cadeia[-(1:burn_in)]Gráfico de Trajetória (Trace Plot)
Cadeias Iniciais (Burn-in)
Um certo número de iterações iniciais, denominado burn-in, deve ser removido da análise, pois a cadeia ainda não atingiu sua distribuição estacionária. O gráfico a seguir mostra um zoom em um intervalo da cadeia para ilustrar o comportamento após a convergência.
Gráfico de Densidades
Autocorrelação (ACF) e Taxa de Aceitação
Com \(\sigma = 0.1\), os passos são muito curtos. A taxa de aceitação muito alta (0.97), pois a nova posição é sempre próxima da anterior. Já o trace-plot mostra que a cadeia se move muito lentamente. A autocorrelação (ACF) decai vagarosamente, indicando amostras altamente correlacionadas e ineficientes. Já o burn-in precisa ser muito longo, pois a convergência é lenta.
Com \(\sigma = 1\), a cadeia é eficiente. A taxa de aceitação razoável (0.73), permitindo explorar o espaço amostral de forma eficiente. A autocorrelação decai rapidamente, o que é ideal. O histograma das amostras se ajusta bem à FDP teórica. O burn-in pode ser curto, pois a cadeia converge rapidamente.
Com \(\sigma = 10\), os passos são muito grandes. A taxa de aceitação muito baixa (0.14), pois a maioria das propostas é rejeitada. O trace-plot mostra que a cadeia fica travada no mesmo valor por longos períodos. O burn-in pode ser difícil de avaliar, e as amostras são de baixa qualidade.
Considerações Finais: Ajuste Fino e Redução da Autocorrelação
É importante esclarecer que o objetivo do MCMC é gerar amostras identicamente distribuidas (id).
No caso de \(\sigma=10\), embora a taxa de aceitação seja baixa, podemos usar a autocorrelação a nosso favor. Analisando o gráfico de ACF, notamos que a correlação se torna insignificante (cai para dentro do intervalo de confiança) por volta do lag=50. Isso nos dá uma estratégia clara: se selecionarmos apenas uma amostra a cada 50 iterações, o thinning, obteremos uma sub-amostra com correlação praticamente nula.
O código abaixo demonstra exatamente isso: primeiro, visualizamos a alta correlação na cadeia original e, em seguida, o resultado após aplicar o thinning.
# Define o layout para 2 linhas e 1 coluna
par(mfrow = c(1,2))
# Gráfico 1: ACF da cadeia original (sigma=10)
# Aumentamos o lag.max para ver claramente onde a correlação desaparece
acf(cadeia_pos_burn_in3,
lag.max = 60,
main = "ACF Original (sigma = 10)")
mtext(paste("Taxa de Aceitação:", round(result3$taxa_aceite, 2)),
side = 3, line = 0.5)
# Gráfico 2: ACF da cadeia após o thinning
# Criamos a nova amostra selecionando 1 valor a cada 60
cadeia_final <- cadeia_pos_burn_in3[seq(from = 1, to = length(cadeia_pos_burn_in3), by = 60)]
acf(cadeia_final,
main = "ACF após Thinning (lag=60)")
mtext(paste("Nº de amostras:", length(cadeia_final)),
side = 3, line = 0.5)
Como o segundo gráfico mostra, a autocorrelação na cadeia final foi praticamente eliminada. O resultado é um conjunto de amostras muito mais eficiente para realizar as análises a posteriori, como calcular médias, variâncias ou intervalos de credibilidade.