4.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]5.

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).

4.2.1 Cadeias de Markov

Exemplo 4.6 Pode-se simular o número de clientes em três diferentes categorias.

(v0 <- c(324,233,210))
## [1] 324 233 210
(P <- matrix(c(.7,.2,.1, .1,.6,.3, 0,.1,.9),
             nrow = 3, byrow = TRUE))
##      [,1] [,2] [,3]
## [1,]  0.7  0.2  0.1
## [2,]  0.1  0.6  0.3
## [3,]  0.0  0.1  0.9
rownames(P) <- colnames(P) <- LETTERS[1:nrow(P)]
P
##     A   B   C
## A 0.7 0.2 0.1
## B 0.1 0.6 0.3
## C 0.0 0.1 0.9
(v1 <- t(v0) %*% P)
##          A     B     C
## [1,] 250.1 225.6 291.3
(v2 <- v1 %*% P)
##           A      B      C
## [1,] 197.63 214.51 354.86
(P2 <- P %*% P)
##      A    B    C
## A 0.51 0.27 0.22
## B 0.13 0.41 0.46
## C 0.01 0.15 0.84
t(v0) %*% P2
##           A      B      C
## [1,] 197.63 214.51 354.86

Exemplo 4.7 Pode-se simular uma cadeia irredutível, i.e., aquela todos os estados se comunicam.

(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
# 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
# distribuição de probabilidade dos estados após 2 passos
(x2 <- x1 %*% P)        # simbolicamente, x1*P
##        [,1]   [,2]  [,3]
## [1,] 0.2425 0.3725 0.385
(x2 <- x0 %*% P %*% P)  # simbolicamente, x0*P²
##        [,1]   [,2]  [,3]
## [1,] 0.2425 0.3725 0.385
# a probabilidade de atingir o estado 3 em 2 passos é
# Pr(X_2 = s_3 | X_0 = s_2)
x2[3]
## [1] 0.385

Exemplo 4.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 4.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
rownames(P) <- colnames(P) <- LETTERS[1:nrow(P)]
P
##      A    B   C
## A 0.50 0.20 0.3
## B 0.15 0.45 0.4
## C 0.25 0.35 0.4
e <- 0.01

estac(x0,P,e)
## $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.3712155

4.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.

library(markovchain)

(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
# 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"
Pmc
## 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
# gráfico
plot(Pmc)

plot(P)

# estados após n-passos
initialState <- c(0,1,0)
steps <- 4
finalState <- initialState*Pmc^steps # using power operator
finalState
##              A         B         C
## [1,] 0.2832562 0.3440813 0.3726625
estac(initialState,P,e=0.01)$z
##           [,1]      [,2]      [,3]
## [1,] 0.2832562 0.3440813 0.3726625
estac(initialState,P,e=0.01)$i
## [1] 4
steadyStates(Pmc) # S4 method
##              A         B         C
## [1,] 0.2878788 0.3409091 0.3712121
estac(initialState, P, e=0.0001)
## $i
## [1] 5
## 
## $erro
## [1] 1.500000e-01 1.500000e-02 1.233750e-02 1.435427e-03 1.495009e-05
## 
## $z
##           [,1]      [,2]      [,3]
## [1,] 0.2878311 0.3409418 0.3712271
## 
## $Pn
##           [,1]      [,2]      [,3]
## [1,] 0.2878788 0.3409091 0.3712121
## [2,] 0.2878788 0.3409091 0.3712121
## [3,] 0.2878788 0.3409091 0.3712121

Exercício 4.9 Veja o artigo de (N. Friedman and Koller 2003) e a Seção 9.16.

4.2.2 Metropolis-Hastings

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 4.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 <- rep(runif(1), Nsim)  # inicia a cadeia
Naccept <- 0
for (i in 2:Nsim){
  Y <- runif(1)
  rho <- (dbeta(Y,a,b)/dbeta(X[i-1],a,b)) * (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.4564

Validando.

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.0154, p-value = 0.5936
## alternative hypothesis: two-sided
plot(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')

a/(a+b) # média teórica
## [1] 0.3
mean(X) # média empírica
## [1] 0.297259
(a*b)/((a+b)^2*(a+b+1)) # varância teórica
## [1] 0.021
var(X) # variância empírica
## [1] 0.02069883
par(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)

Exercício 4.10 Considerando o Exemplo 4.5.
a. Faça a simulação de uma exponencial de parâmetro \(\lambda=1\) utilizando o algoritmo de Metropolis-Hastings.
b. Compare os métodos da rejeição e de Metropolis-Hastings quanto ao tempo de processamento. O que você percebe?

Sugestão: Capítulo 12.

Referências

Friedman, Nir, and Daphne Koller. 2003. “Being Bayesian about Network Structure. A Bayesian Approach to Structure Discovery in Bayesian Networks.” Machine Learning 50: 95–125. https://link.springer.com/content/pdf/10.1023/A:1020249912095.pdf.
Gamerman, Dani, and Hedibert F Lopes. 2006. Markov Chain Monte Carlo: Stochastic Simulation for Bayesian Inference. CRC press.
Gelfand, Alan E, and Adrian FM Smith. 1990. “Sampling-Based Approaches to Calculating Marginal Densities.” Journal of the American Statistical Association 85 (410): 398–409. https://www.tandfonline.com/doi/abs/10.1080/01621459.1990.10476213.
Geman, Stuart, and Donald Geman. 1984. “Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images.” IEEE Transactions on Pattern Analysis and Machine Intelligence, no. 6: 721–41. https://cs.uwaterloo.ca/~mannr/cs886-w10/GemanandGeman84.pdf.
Hastings, W Keith. 1970. “Monte Carlo Sampling Methods Using Markov Chains and Their Applications.” https://www.jstor.org/stable/pdf/2334940.pdf?casa_token=ab0bhhoQvbcAAAAA:1ksPSl_D_2e_zY3ijO30RLGwUFF3KYeF0b6gBH5cJemjei8tDnceRX7CC0h-sREBZFTjZGwmQHCjObmgge83_a-UpFBV9rS_vl36pq5MHvQF6DlrS7o.
Metropolis, Nicholas, Arianna W Rosenbluth, Marshall N Rosenbluth, Augusta H Teller, and Edward Teller. 1953. “Equation of State Calculations by Fast Computing Machines.” The Journal of Chemical Physics 21 (6): 1087–92. https://doi.org/10.1063/1.1699114.
Peng, Roger D. 2022. Advanced Statistical Computing. Work in Progress. https://bookdown.org/rdpeng/advstatcomp/.
Pettitt, Anthony N, and Candice M Hincksman. 2013. “Introduction to MCMC.” Case Studies in Bayesian Statistical Modelling and Analysis, 17–29.
Robert, Christian P, and George Casella. 2010. Introducing monte carlo methods with R. Springer.
Spedicato, Giorgio Alfredo. 2017. “Discrete Time Markov Chains with r.” The R Journal, July. https://journal.r-project.org/archive/2017/RJ-2017-036/index.html.

  1. (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)↩︎