5.2 MCMC
Os métodos de Monte Carlo via Cadeia de Markov (Markov Chain Monte Carlo ou MCMC) foram popularizados por (Gelfand and Smith 1990), quando “apontaram com sucesso para a comunidade estatística em geral que o esquema de amostragem desenvolvido por (Geman and Geman 1984) para distribuições de Gibbs poderia de fato ser usado para uma série de outras distribuições posterioris” (Gamerman and Lopes 2006, 141)6.
Apesar de (Hastings 1970) ter generalizado o algoritmo de (Metropolis et al. 1953) – o que ficou conhecido como algoritmo de Metropolis-Hastings –, (Pettitt and Hincksman 2013) indicam que “a ideia permaneceu sem uso na literatura estatística”. Os autores apontam que as abordagens MCMC possuem “uma clara conexão entre o algoritmo de maximização da expectativa (EM) (…) e a amostragem de Gibbs” e “foram desenvolvidas para modelos onde há variáveis latentes usadas na verossimilhança, como modelos mistos, modelos de mistura e modelos para processos estocásticos”.
Para detalhes dos métodos, veja o Capítulo 7 de (Peng 2022).
5.2.1 Cadeias de Markov
Exemplo 5.6 Pode-se simular o número de clientes em três diferentes categorias.
## [1] 324 233 210## [,1] [,2] [,3]
## [1,] 0.7 0.2 0.1
## [2,] 0.1 0.6 0.3
## [3,] 0.0 0.1 0.9## A B C
## A 0.7 0.2 0.1
## B 0.1 0.6 0.3
## C 0.0 0.1 0.9## A B C
## [1,] 250.1 225.6 291.3## A B C
## [1,] 197.63 214.51 354.86## A B C
## A 0.51 0.27 0.22
## B 0.13 0.41 0.46
## C 0.01 0.15 0.84## A B C
## [1,] 197.63 214.51 354.86Exemplo 5.7 Pode-se simular uma cadeia irredutível, i.e., aquela todos os estados se comunicam.
## [,1] [,2] [,3]
## [1,] 0.50 0.20 0.3
## [2,] 0.15 0.45 0.4
## [3,] 0.25 0.35 0.4# iniciando no estado 2
x0 <- c(0,1,0)
# distribuição de probabilidade dos estados após 1 passo
(x1 <- x0 %*% P) # simbolicamente, x0*P## [,1] [,2] [,3]
## [1,] 0.15 0.45 0.4## [,1] [,2] [,3]
## [1,] 0.2425 0.3725 0.385## [,1] [,2] [,3]
## [1,] 0.2425 0.3725 0.385## [1] 0.385Exemplo 5.8 Crie uma função que encontre a distribuição estacionária (steady state) em função de x0, P e do erro, apresentando o número dw iterações e o vetor de erros. Aplique para o Exemplo 5.7 com e=0.01.
# função
estac <- function(x0,P,e){
z <- x0
Pn <- P
i <- 1
erro <- min(abs(z - x0%*%Pn))
while(erro[i] > e){
z <- x0 %*% Pn
Pn <- Pn %*% Pn
i <- i + 1
erro[i] <- min(abs(z - x0%*%Pn))
}
return(list(i=i,erro=erro,z=z,Pn=Pn))
}
# dados do exemplo
x0 <- c(0,1,0)
(P <- matrix(c(.5,.2,.3, .15,.45,.4, .25,.35,.4),
nrow = 3, byrow = TRUE))## [,1] [,2] [,3]
## [1,] 0.50 0.20 0.3
## [2,] 0.15 0.45 0.4
## [3,] 0.25 0.35 0.4## A B C
## A 0.50 0.20 0.3
## B 0.15 0.45 0.4
## C 0.25 0.35 0.4## $i
## [1] 4
##
## $erro
## [1] 0.150000000 0.015000000 0.012337500 0.001435427
##
## $z
## A B C
## [1,] 0.2832562 0.3440813 0.3726625
##
## $Pn
## A B C
## A 0.2879490 0.3408609 0.3711901
## B 0.2878311 0.3409418 0.3712271
## C 0.2878681 0.3409164 0.37121555.2.1.1 Biblioteca markovchain
A biblioteca markovchain (Spedicato 2017) oferece funções e métodos para criar e gerenciar cadeias de Markov de tempo discreto. A vinheta fornece alguns exemplos aplicados.
## [,1] [,2] [,3]
## [1,] 0.50 0.20 0.3
## [2,] 0.15 0.45 0.4
## [3,] 0.25 0.35 0.4# transformando P em 'markovchain'
Pmc <- new('markovchain', transitionMatrix = P,
states = LETTERS[1:ncol(P)],
name='MarkovChain P')
class(Pmc)## [1] "markovchain"
## attr(,"package")
## [1] "markovchain"## MarkovChain P
## A 3 - dimensional discrete Markov Chain defined by the following states:
## A, B, C
## The transition matrix (by rows) is defined as follows:
## A B C
## A 0.50 0.20 0.3
## B 0.15 0.45 0.4
## C 0.25 0.35 0.4
5.2.2 Metropolis-Hastings
Será considerado o algoritmo de Metropolis-Hastings descrito por (Robert and Casella 2010, 171). Está associado à densidade (alvo) \(f\) e à densidade condicional \(q\), produzindo uma cadeia de Markov \(X^{(t)}\) através do seguinte kernel de transição:
Dado \(x^{(t)}\),
1. Gere \(Y_t \sim q(y|x^{(t)})\).
2. Faça
\[X^{(t+1)} = \left\{
\begin{array}{l l}
Y_t & \text{com probabilidade} & \rho(x^{(t)}, Y_t) \\
x^{(t)} & \text{com probabilidade} & 1-\rho(x^{(t)}, Y_t) \\ \end{array} \right. \]
onde
\[ \rho(x,y) = \min \left\{ \frac{f(y)}{f(x)} \frac{q(x|y)}{q(y|x)}, 1 \right\} \]
A distribuição \(q\) é chamada de distribuição instrumental/proposta/candidata e \(\rho(x,y)\) é a probabilidade de aceitação de Metropolis-Hastings.
Exemplo 5.9 (Robert and Casella 2010, 172) simulam uma \(\mathcal{Beta}(2.7,6.3)\) a partir de uma \(\mathcal{U}(0,1)\).
a <- 2.7; b <- 6.3 # valores iniciais
Nsim <- 5000 # número de simulações
X <- runif(Nsim) # inicia a cadeia
Naccept <- 0
for (i in 2:Nsim){
Y <- runif(1)
rho <- min((dbeta(Y,a,b)/dbeta(X[i-1],a,b)) * (1/1), 1)
accept <- runif(1) < rho
X[i] <- X[i-1] + (Y - X[i-1])*accept
Naccept <- Naccept+accept
}
Naccept/Nsim # taxa de aceitação## [1] 0.452Validando.
bt <- rbeta(Nsim,a,b) # beta via stats::rbeta
ks.test(jitter(X), bt) # teste de aderência de Kolmogorov-Smirnov##
## Asymptotic two-sample Kolmogorov-Smirnov test
##
## data: jitter(X) and bt
## D = 0.0244, p-value = 0.1019
## alternative hypothesis: two-sidedplot(ecdf(X), col = 'blue', main = 'Função de distribuição acumulada empírica')
plot(ecdf(bt), col = 'red', add = TRUE)
legend(0, .9, legend=c('Metropolis-Hastings', 'stats::rbeta'), lty = 1,
col=c("red", "blue"), cex=0.8, bty = 'n')
## [1] 0.3## [1] 0.2966724## [1] 0.021## [1] 0.02053429par(mfrow=c(1,2))
hist(X, freq = FALSE, main = 'Beta(2.7,6.3) via M-H')
curve(dbeta(x,a,b), col = 'blue', add = TRUE)
hist(bt, freq = FALSE, main = 'Beta(2.7,6.3) via stats::rbeta')
curve(dbeta(x,a,b), col = 'red', add = TRUE)
Referências
(Gelfand and Smith 1990) were the first authors to successfully point out to the statistical community at large that the sampling scheme devised by (Geman and Geman 1984) for Gibbs distributions could in fact be used for a host of other posterior distributions. (Gamerman and Lopes 2006, 141)↩︎