No Hackeando o R de hoje, vamos apresentar o método de MCMC, ferramenta de grande importância para a estimação dentro da estatística bayesiana. O exemplo que apresentamos na semana passada era de fácil computação (como pode ser visto pelos códigos utilizados), muito por causa do fato de que a equação encontrada possuía solução fechada. Quando possuímos poucos parâmetros para estimar, e o espaço de possibilidades dos parâmetros é pequeno, também podemos resolver 'manualmente' calculando as probabilidades, pois é questão de percorrer cada uma das possibilidades e encontrar o valor do Teorema de Bayes. Para problemas mais complexos, como os que encontramos na vida real, essas soluções fáceis muitas vezes não estão disponíveis. Pense no caso de uma probabilidade condicional sem fórmula fechada, para a estimação de 6 parâmetros, que podem assumir 1000 valores cada um. Estaríamos olhando para o cálculo de 1000^6 combinações, o que é pesado até para computadores modernos.
A solução para esse problema é o chamado MCMC (Monte-Carlo Markov Chain). Apesar dele ser mais utilizado em problemas com múltiplos parâmetros, vamos mostrar como ele funciona com apenas um, para facilitar a compreensão. O ponto importante do método é que ele só depende da estimação dos termos no numerador do Teorema de Bayes (mostrado abaixo), eliminando qualquer preocupação com a integral do denominador, que é o grande problema para encontrarmos soluções fechadas.
A ideia básica do MCMC é gerar uma estimativa da distribuição a posteriori a partir de amostras dela, sem que nós precisemos construir ela per se. Hoje, vamos mostrar o algoritmo de Metropolis, que gera uma aproximação da distribuição a posteriori através de regras de decisão simples. O algoritmo utiliza apenas a razão entre duas probabilidades a posteriori, logo há estimação da integral acaba sendo desnecessária. Dada uma posição inicial para os parâmetros, os passos do algoritmo são simples:
1: geramos um movimento aleatório dos parâmetros no seu espaço, que será testado;
2: Verificamos se a posteriori na nova posição tem valor maior que na original. Se sim, o movimento ocorre. Caso contrário, o movimento ocorrerá com probabilidade igual à razão entre as posterioris;
3: geramos uma observação uniforme de 0 a 1, e comparamos ela com a razão calculada, de modo a validar o movimento acima.
Com isso, vamos testar o algoritmo para uma distribuição de eventos Bernoulli. Sabemos que essa distribuição possui fórmula fechada, logo iremos utilizá-la para comparar com a distribuição 'aproximada' encontrada de modo empírico.
dados = c(rep(0,6),rep(1,14))
probParcial = function(theta, data){
z = sum(data)
N = length(data)
p_Xi_dado_Theta = theta^z * (1-theta)^(N-z)
p_Xi_dado_Theta[theta > 1 | theta < 0] = 0
pTheta = dbeta(theta , 1, 1)
pTheta[theta > 1 | theta < 0] = 0
parcial = p_Xi_dado_Theta * pTheta
return(parcial)
}
n = 50000
trajetoria = rep(0, n)
trajetoria[1] = 0.5
set.seed(12334)
for ( t in 1:(n-1) ) {
posicao = trajetoria[t]
choque = rnorm(1, mean = 0, sd = 0.2)
probabilidade = min(1,
probParcial(posicao + choque, dados)
/ probParcial(posicao, dados))
if ( runif(1) < probabilidade ) {
trajetoria[t+1] = posicao + choque
} else {
trajetoria[t+1] = posicao
}
}
df = data.frame(passo = 1:50000, valor = trajetoria)
library(ggplot2)
ggplot(df, aes(x=valor)) + geom_histogram(aes(y = stat(density))) +
stat_function(fun = function(x) dbeta(x, 15, 7), color = "red",
size = 1)
________________________
(*) Para entender mais sobre análises estatísticas, confira nosso Curso de Estatística Bayesiana usando o R.
