O algoritmo Metropolis-Hastings

O que é?

O algoritmo Metropolis-Hastings é um Método de Monte Carlo via cadeias de Markov (MCMC).

Cadeias de Markov

Seja (T_n)_{n \in \mathbb{N}} uma sequência de variáveis discretas. Dizemos que (T_n)_{n \in \mathbb{N}} é uma Cadeia de Markov se, para todo n \in \mathbb{N},

\begin{align*} &\mathbb{P}(T_{n+1} = t_{n+1} | T_n = t_n, T_{n-1} = t_{n-1}, \dots, T_1 = t_1) = \\ &\mathbb{P}(T_{n+1} = t_{n+1} | T_n = t_n) \end{align*}

Cadeia homogênea

Seja (T_n)_{n \in \mathbb{N}} uma cadeia de Markov. Dizemos que (T_n)_{n \in \mathbb{N}} é homogênea se, \mathbb{P}(T_{n+1} = j | T_n = i) = \mathbb{P}(T_{1} = j | T_0 = i), \forall n,i,j.

A partir deste ponto, consideraremos apenas cadeias de Markov homogêneas.

Matriz de transição

Seja (T_n)_{n \in \mathbb{N}} uma cadeia de Markov homogênea. \mathbb{P} é a matriz de transição para (T_n)_{n \in \mathbb{N}} se

\mathbb{P}_{i,j} = \mathbb{P}(T_{n+1} = j | T_n = i) = p(j|i).

Exemplo

Considere uma cadeia de Markov com estados \{0,1\} com a seguinte matriz de transição:

\mathbb{P} = \begin{bmatrix} \frac{1}{3} & \frac{2}{3} \\ \frac{1}{2} & \frac{1}{2} \end{bmatrix}

Distribuição estacionária

Seja (T_n)_{n \in \mathbb{N}} uma cadeia de Markov com matriz de transição, \mathbb{P}. Para cada estado, i, definimos

\mu(i) = \lim_{n \to \infty} \mathbb{P}(T_n = i | T_0 = t_0) = \lim_{n \to \infty} (\mathbb{P}^n)_{t_0, i}.

Se o limite for bem definido, dizemos que \mu é a distribuição estacionária da Cadeia.

Distribuição invariante

Seja (T_n)_{n \in \mathbb{N}} uma cadeia de Markov com matriz de transição, \mathbb{P}. Seja \mu um vetor não-negativo tal que \sum_i \mu(i) = 1. \mu é a distribuição invariante para \mathbb{P} se

\mu \mathbb{P} = \mu.

Distribuição reversível

Seja (T_n)_{n \in \mathbb{N}} uma cadeia de Markov com matriz de transição \mathbb{P}. Seja \mu um vetor positivo tal que \sum_i \mu(i) = 1. \mu é a distribuição reversível para \mathbb{P} se, para todos estados i e j,

\mu(i) \mathbb{P}_{i,j} = \mu(j) \mathbb{P}_{j,i}.

Lei dos Grandes Números para Cadeias de Markov

Seja (T_n)_{n \in \mathbb{N}} uma cadeia de Markov irredutível e aperiódica com distribuição estacionária \mu. Seja g uma função tal que \mathbb{E}_\mu[|g(T)|] < \infty. Então, para qualquer estado inicial T_0,

\mathbb{E}_\mu[g(T)] = \int g(\theta)\mu(\theta)d\theta\approx \frac{1}{B}\sum_{i=1}^B g(T_i).

Cadeias de Markov e Posteriori

Se você construir uma cadeia de Markov, T_1, \dots, T_B, com distribuição estacionária f(\theta|x), então

\frac{\sum_{i=1}^{B} g(T_i)}{B}

é um estimador consistente para

\int g(\theta)f(\theta|x)d\theta = E[g(\theta)|X].

O Rei Markov (A Intuição)

Um rei governa 6 ilhas dispostas em anel, onde a população da ilha i é Pop_i = i.

O Objetivo: A frequência de longo prazo das visitas deve ser proporcional à população de cada ilha.

O Rei Markov (A Intuição)

A Estratégia de Transição:

1. Proposta: O rei joga uma moeda para decidir se propõe ir à ilha vizinha (esquerda ou direita).
2. Avaliação: - Se Pop_{proposta} > Pop_{atual}, ele aceita e se move.
   • Se Pop_{proposta} < Pop_{atual}, ele move-se com probabilidade R = \frac{Pop_{proposta}}{Pop_{atual}}.
3. Permanência: Caso não aceite o movimento, ele permanece na ilha atual por mais um dia.

O Algoritmo de Metropolis-Hastings: passo 1

Para um parâmetro contínuo \theta com densidade alvo (posteriori) f(\theta|x):

1. Defina um valor inicial T_1.

O Algoritmo de Metropolis-Hastings: passo 2

2. Para i = 2, \dots, B:
   1. Gere um candidato T_i^* de uma distribuição de proposta h(t_i^* | T_{i-1}).
   2. Calcule a razão de aceitação: R_i = \frac{f(T_i^* | x) \cdot h(T_{i-1} | T_i^*)}{f(T_{i-1} | x) \cdot h(T_i^* | T_{i-1})}
   3. Gere U_i \sim \text{Uniforme}(0, 1).
   4. Se U_i \le R_i, defina T_i = T_i^*. Caso contrário, T_i = T_{i-1}.

Exemplo 1: Normal-Normal

Considere o modelo: X|\theta \sim N(\theta, 1) com x=10 e prior \theta \sim N(0, 1).

Função de Proposta: Utilizaremos um Random Walk Gaussiano. O candidato \theta^* é gerado a partir do estado atual \theta^{(i-1)}: \theta^* | \theta^{(i-1)} \sim N(\theta^{(i-1)}, \sigma^2_{prop})

No R: proposta <- rnorm(1, mean = estado_atual, sd = 1)

Exemplo 1: Nota sobre Simetria no caso Normal-Normal

Como h(\theta^* | \theta) = h(\theta | \theta^*), a razão de Hastings simplifica-se para 1 (ou 0 em log). O que acontece com a taxa de aceitação se alterarmos o sd da proposta para 0.1 ou 100?

Exemplo 2: Beta-Binomial

Considere o modelo de verossimilhança Binomial: X | \theta \sim \text{Binomial}(n=40, \theta) Dados: Observamos x = 26 sucessos.

Priori: \theta \sim \text{Beta}(\alpha=1, \beta=10).

Função de Proposta (Independente): Desta vez, o candidato \theta^* não depende do estado atual: \theta^* \sim \text{Beta}(2, 2)

Exemplo 2: Beta-Binomial: Razão de Hastings

Como h(\theta^* | \theta) não depende de \theta, a razão de aceitação deve compensar a probabilidade de termos sorteado o candidato: R = \frac{f(\theta^*|x) \cdot h(\theta^{(i-1)})}{f(\theta^{(i-1)}|x) \cdot h(\theta^*)}

No R: proposta <- rbeta(1, 2, 2)