Métodos de Monte Carlo via Cadeias de Markov

Motivação

O Problema da Intratabilidade

Em diversas áreas, tais como Estatística Bayesiana, Física Estatística e Aprendizado de Máquinas, frequentemente nos deparamos com distribuições de probabilidade \(f(\boldsymbol{x})\) das quais desejamos extrair amostras ou calcular esperanças da forma:

\[\mathbb{E}[g(\boldsymbol{X})] = \int g(\boldsymbol{x}) f(\boldsymbol{x}) \, d\boldsymbol{x}\]

Contudo, na maioria dos problemas complexos de alta dimensão, a densidade \(f(\boldsymbol{x})\) é conhecida apenas a menos de uma constante de normalização:

\[f(\boldsymbol{x}) = \frac{h(\boldsymbol{x})}{C}, \quad \text{onde } C = \int h(\boldsymbol{x}) \, d\boldsymbol{x}\]

O cálculo de \(Z\) pode envolver integrais multidimensionais de altíssima dimensão que não possuem solução analítica e inviabilizam métodos tradicionais de integração numérica (como a regra do trapézio ou a quadratura gaussiana).

Por que utilizar Cadeias de Markov?

Os métodos de Amostragem por Rejeição ou Amostragem de Importância exigem uma função de proposta envelope que cubra adequadamente a distribuição alvo \(f(\boldsymbol{x})\). Em dimensões elevadas, encontrar essa proposta torna-se praticamente impossível.

A ideia central do MCMC (Markov Chain Monte Carlo) é explorar o espaço de estados de forma inteligente: construímos uma Cadeia de Markov discreta no tempo com espaço de estados contínuo/discreto cuja distribuição estacionária (ou invariante) seja exatamente a distribuição alvo \(f(\boldsymbol{x})\).

---

Introdução aos Métodos de Monte Carlo via Cadeias de Markov

Continuamos querendo obter uma amostra da distribuição de interesse para calcular estimativas amostrais (como nos métodos de Monte Carlo usuais); no entanto, em todos os algoritmos de simulação vistos até agora, os valores são gerados de forma independente e não há preocupação com a convergência do algoritmo (estes métodos são chamados não-iterativos).

Este problema se torna cada vez mais complexo, a medida que aumentamos a dimensão da distribuição de interesse. De fato, nem sempre é fácil (ou mesmo possível) encontrar uma distribuição auxiliar boa (no sentido de aproximar bem a distribuição de interesse) e fácil de ser amostrada.

Uma alternativa neste caso é utilizar métodos iterativos de simulação; por exemplo, baseados em Cadeias de Markov, para otimizar o processo de simulação. No entanto, neste caso, os valores gerados não serão mais independentes, e deverão ser tratados com certo cuidado.

Os métodos de Monte Carlo via Cadeias de Markov (MCMC) permitem simular da distribuição de interesse de forma indireta (novamente utilizando distribuições auxiliares). A ideia geral é construir uma cadeia de Markov fácil de ser simulada e com distribuição de equilíbrio igual a distribuição de interesse. Após um número suficientemente grande de iterações, a cadeia converge para a distribuição de interesse.

---

Cadeias de Markov

Definição: Uma cadeia de Markov é uma sequência de variáveis aleatórias \(\{X_1,X_2,...\}\) tal que para todo evento \(A\), \[P(X_t \in A|X_1,...,X_{t-1})=P(X_t \in A|X_{t-1}), \; \forall t,\] ou seja, uma v.a. num certo instante só depende do instante imediatamente anterior.

Algumas propriedades importantes

Métodos de MCMC utilizam cadeias que sejam

• homogêneas \(\rightarrow\) as probabilidades de transição não mudam ao longo das iterações;

• irredutíveis \(\rightarrow\) não existem estados ``isolados”; e

• aperiódicas \(\rightarrow\) sem estados absorventes.

Segundo Robert e Casella (2004), um Método MCMC para simulação de uma distribuição \(f\) é qualquer método que produz uma cadeia de Markov ergótica^\((\star)\) cuja distribuição estacionária é \(f\).

^\((\star)\) Uma cadeia de Markov ergótica é um processo estocástico em que todos os estados são irredutíveis e aperiódicos.

---

Referência importante

Por que o algoritmo produz uma cadeia cuja distribuição estacionária é a de interesse? Veja

Understanding the Metropolis-Hastings Algorithm
S. Chib & E. Greenberg
The American Statistician
1995, 49, 4, 327–335

---

Algoritmo Metropolis-Hastings

Proposto originalmente por Metropolis et al. (1953) e generalizado por Hastings (1970), este algoritmo utiliza a mesma ideia do método de rejeição: gera de uma distribuição auxiliar \(q(\boldsymbol{x})\) e aceita com uma certa probabilidade que depende da distribuição de interesse \(f(\boldsymbol{x})\).

A distribuição auxiliar \(q(\boldsymbol{x})\) é chamada de distribuição proposta e pode depender do estado atual da cadeia (\(\boldsymbol{x}^c\)). A ideia é propor um novo estado para a cadeia, digamos, \(\boldsymbol{x}^p \sim q(.|\boldsymbol{x}^c)\), e aceitá-lo com probabilidade

\[ \alpha = \min\left\{ 1, \frac{f(\boldsymbol{x}^p)}{f(\boldsymbol{x}^c)}\frac{q(\boldsymbol{x}^c|\boldsymbol{x}^p)}{q(\boldsymbol{x}^p|\boldsymbol{x}^c)} \right\} \]

Note na expressão acima que basta conhecer o núcleo de \(f\) para sabermos calcular a probabilidade de aceitação!

Algoritmo

1. Inicialize a cadeia com um valor inicial \(\boldsymbol{x}^{(0)}\) e inicialize um contador \(t=1\);

2. gere um novo valor \(\boldsymbol{x}^p \sim q(.|\boldsymbol{x}^{(t-1)})\);

3. calcule \(\alpha\) e gere \(u \sim U(0,1)\);

4. se \(u<\alpha\), aceite \(\boldsymbol{x}^p\) fazendo \(\boldsymbol{x}^{(t)}=\boldsymbol{x}^p\); caso contrário, rejeite e faça \(\boldsymbol{x}^{(t)}=\boldsymbol{x}^{(t-1)}\);

5. faça \(t=t+1\) e volte ao passo \(2\) até atingir a convergência.

Os valores simulados para a cadeia após a convergência representarão uma amostra da distribuição de interesse \(f(\boldsymbol{x})\).

Observações importantes

1. Os valores iniciais influenciam o comportamento da cadeia, porém, esperamos que conforme o número de iterações aumente, a cadeia gradativamente esqueça os valores iniciais e convirja. Na prática é comum descartar as iterações iniciais (chamadas de amostra de aquecimento, ou, em inglês burn-in).

2. Em uma cadeia de Markov, os valores gerados são por definição correlacionados ao longo das iterações (pois o valor de \(\boldsymbol{x}^{(t)}\) depende de \(\boldsymbol{x}^{(t-1)}\)). Isto significa que pode não haver muito ganho em termos de informação em se armazenar todos os valores simulados da cadeia. Uma abordagem prática para economizar espaço de armazenamento no disco e diminuir a dependência entre os valores gerados é guardar os valores da cadeia a cada certa quantidade de iterações (essa quantidade é conhecida como espaçamento ou thinning). Espera-se que os valores guardados fiquem menos correlacionados.

---

Diagnósticos de Convergência

A análise de convergência em métodos de MCMC é crucial para garantir que as amostras geradas se aproximam da distribuição alvo e, para isso, é necessário que haja boas propriedades das cadeias (estacionariedade - homogeneidade e ergodicidade - irredutibilidade e aperiodicidade).

Os métodos de MCMC consideram que após um certo período (depois de muitas iterações) a cadeia converge para a distribuição alvo. No entanto, a pergunta que fica é: Como avaliar se a convergência foi atingida?

Antes dos diagnósticos formais, podemos realizar uma série de avaliações a partir dos valores gerados que dão indícios sobre a convergência do algoritmo.

1. Análise gráfica: nesta primeira etapa realizamos a análise dos traços das cadeias obtidas, verificando se apresentam comportamento repetitivo ao longo das iterações e se exploram de forma adequada toda a imagem das variáveis de interesse.

2. Retirada do aquecimento: como o estado inicial \(x^{(0)}\) é arbitrário, as primeiras iterações da cadeia podem depender fortemente de \(x^{(0)}\) e ainda não representam amostras da distribuição limite \(f(x)\). Neste caso, para garantir que a amostra não seja afetada pelos valores iniciais das cadeias, é necessário remover as iterações iniciais, também conhecida como amostra de aquecimento. (Tipicamente descarta-se de 10% a 50% das amostras.)

3. Análise da autocorrelação: após eliminar o período de aquecimento, avalia-se a autocorrelação entre os alores da cadeia. Correlações altas devem ser eliminada tomando espaçamento entre valores consecutivos nas cadeias.

4. Tamanho Efetivo de Amostra: por construção, as amostras de uma Cadeia de Markov são correlacionadas (\(X_t\) depende de \(X_{t-1}\)). A variância da média amostral sob dados autocorrelacionados é dada por:

\[\text{Var}(\bar{X}_M) = \frac{\sigma^2}{M} \left( 1 + 2 \sum_{k=1}^{\infty} \rho_k \right)\]

O Tamanho Efetivo de Amostra (ESS - do inglês effective sample size) mede a quantidade equivalente de observações independentes:

\[ESS = \frac{M}{1 + 2 \sum_{k=1}^{\infty} \rho_k}\]

Não existe um valor mínimo universal para o ESS, no entanto, valores baixos desta medida com relação ao número de iterações (por exemplo, \(ESS/M<0,10\)) indicam amostras muito correlacionadas, o que significa que o processo de amostragem não foi eficiente. No entanto, ESS baixo não implica necessariamente falta de convergência. Uma cadeia pode ter convergido para a distribuição estacionária e, ainda assim, ser altamente autocorrelacionada.

5. Múltiplas cadeias: outra boa prática é fazer a análise de múltiplas cadeias, partindo de valores iniciais distintos e verificando se todas apresentam comportmentos semelhantes, dando indícios de convergência.

6. Testes formais: a convergência também pode ser verificada numericamente através de testes formais. Dois testes possíveis são os propostos por

• Geweke (1992): verifica se amostras no início e no fim da cadeia são similares; e

• Gelman e Rubin (1992): utilizado quando existem múltiplas cadeias a serem avaliadas, é um teste formal que compara a variabilidade entre as diferentes cadeias de Markov simuladas com a variabilidade dentro de cada cadeia. De fato, se as cadeias convergiram para a mesma distribuição, a variabilidade intra-cadeia deve ser semelhante à variabilidade inter-cadeia.

---

Exemplo 1: Distribuição Inversa Gama

Gerar valores para \(X\) considerando a distribuição Inversa Gama com parâmetros \(a\) e \(b\), cuja função de densidade de probabilidade é dada por \[f(x) = \dfrac{b^a}{\Gamma(a)}x^{-a-1}\exp\left\{-\frac{b}{x}\right\}, \; x \in (0,\infty),\] usando o algoritmo Metropolis-Hastings.

Vamos implementar o algoritmo utilizando diferentes distribuições auxiliares.

Proposta Uniforme

Inicialmente, vamos supor que sabemos gerar apenas de um gerador uniforme.

Neste caso, utilizando uma variável auxiliar \(U \sim U(0,1)\), podemos fazer a transformação \(Y = \frac{U-a}{b-a}\) para obter amostras de \(Y \sim U(a,b)\).

Note que nosso interesse é gerar de um modelo com imagem infinita (\(Im(X)=(0,\infty)\)); no entanto, na impossibilidade de gerar de uma distribuição com imagem adequada, vamos tomar a aproximação \((0,L)\), com \(L\) grande, de forma a cobrir grande parte da cauda da distribuição de interesse.

  f = function(x,a,b){
    (b^a)/(gamma(a)) * x^(-a-1)*exp(-b/x)*I(x>0)
  }
  a = 2
  b = 3
  curve(f(x,a=a,b=b),0,30,xlab=expression(x),ylab=expression(f(x)),lwd=2)

  f(30,a,b)

## [1] 0.0003016125

Note que a distribuição de interesse possui densidade decrescente a partir de sua moda e que \(f(30)\approx0,0003\); portanto, utilizaremos como ditribuição auxiliar a uniforme \(\text{U}(0,100)\), de forma a desprezar uma cauda com probabilidade bem pequena.

Neste caso, \[ \alpha = \min\left\{ 1, \frac{f(\boldsymbol{x}^p)}{f(\boldsymbol{x}^c)}\frac{q(\boldsymbol{x}^c|\boldsymbol{x}^p)}{q(\boldsymbol{x}^p|\boldsymbol{x}^c)} \right\} = \min\left\{ 1, \frac{\dfrac{b^a}{\Gamma(a)}(x^p)^{-a-1}\exp\left\{-\frac{b}{x^p}\right\}}{\dfrac{b^a}{\Gamma(a)}(x^c)^{-a-1}\exp\left\{-\frac{b}{x^c}\right\}}\frac{\frac{1}{L}}{\frac{1}{L}} \right\} = \min\left\{ 1, \frac{(x^p)^{-a-1}\exp\left\{-\frac{b}{x^p}\right\}}{(x^c)^{-a-1}\exp\left\{-\frac{b}{x^c}\right\}} \right\} \]

  # para garantir reprodutibilidade
  set.seed(1)

  #--------------------------#
  # Usando proposta uniforme #
  #--------------------------#
  
  # criando um vetor para receber a cadeia de valores gerados
  x = NULL
  # inicializando um contador
  cont = 0
  
  # numero de iteracoes
  M = 10000
  
  # inicializando a cadeia
  x[1] = 0.1
  
  for(t in 2:M){
    up = runif(1,0,1)
    x.p = 100*up
    razao.num = f(x.p,a=a,b=b)
    razao.den = f(x[t-1],a=a,b=b)
    razao = razao.num/razao.den
    alpha = min(1,razao)
    u = runif(1,0,1)
    if (u<alpha){
      x[t] = x.p
      cont = cont+1
    }else{
      x[t] = x[t-1]
    }
  }  

  # grafico da cadeia
  ts.plot(x,xlab="iteracoes",ylab=expression(x))

  # tirando o burn-in
  ts.plot(x[1000:M],xlab="iteracoes",ylab=expression(x))

  # grafico da densidade
  hist(x[1000:M],xlab=expression(x),ylab=expression(f(x)),freq=FALSE,nclass=100,main="Histograma dos pontos gerados")
  curve(f(x,a=a,b=b),col=2,lwd=2,add=TRUE)

  # fazendo um recorte da cauda para melhorar a visualização
  hist(x[1000:M],xlab=expression(x),ylab=expression(f(x)),freq=FALSE,nclass=100,main="Histograma dos pontos gerados",xlim=c(0,30))
  curve(f(x,a=a,b=b),col=2,lwd=2,add=TRUE)  

  # taxa de aceitacao
  cont/M    # ficou baixa pois a dist. proposta nao foi boa

## [1] 0.0482

Quando a distribuição proposta é ruim, a cadeia tem dificuldades para se mover, levando a uma taxa de aceitação muito baixa.

Proposta Gama

Vamos melhorar a distribuição proposta utilizando uma curva mais próxima do modelo de interesse.

Neste caso, vamos repetir o exercício utilizando uma distribuição auxiliar/proposta \(\text{Gama}(\gamma,\lambda)\).

  curve(f(x,a=a,b=b),0,30,xlab=expression(x),ylab=expression(f(x)),lwd=2)
  gama = 2/3
  lambda = 1/3
  curve(dgamma(x,gama,lambda),lwd=2,lty=2,col=2,add=TRUE)

Note que a proposta não cobre a distribuição de interesse, mas os comportamentos são semelhantes, no sentido de que regiões de valores pouco/muito prováveis na densidade de interesse também são pouco/muito prováveis na densidade auxiliar.

Utilizando a proposta \(\text{Gama}(\gamma,\lambda)\), \[ \alpha = \min\left\{ 1, \frac{f(\boldsymbol{x}^p)}{f(\boldsymbol{x}^c)}\frac{q(\boldsymbol{x}^c|\boldsymbol{x}^p)}{q(\boldsymbol{x}^p|\boldsymbol{x}^c)} \right\} = \min\left\{ 1, \frac{ \dfrac{b^a}{\Gamma(a)}(x^p)^{-a-1}\exp\left\{-\frac{b}{x^p}\right\} }{ \dfrac{b^a}{\Gamma(a)}(x^c)^{-a-1}\exp\left\{-\frac{b}{x^c}\right\} } \frac{ \dfrac{\lambda^\gamma}{\Gamma(\gamma)}(x^c)^{\gamma-1}\exp\left\{-\lambda{(x^c)}\right\} }{ \dfrac{\lambda^\gamma}{\Gamma(\gamma)}(x^p)^{\gamma-1}\exp\left\{-\lambda{(x^p)}\right\} } \right\} \]

  # para garantir reprodutibilidade
  set.seed(1)

  #----------------------#
  # Usando proposta gama #
  #----------------------#
  
  # criando um vetor para receber a cadeia de valores gerados
  x = NULL
  # inicializando um contador
  cont = 0
  
  # numero de iteracoes
  M = 10000
  
  # inicializando a cadeia
  x[1] = 0.1
  
  for(t in 2:M){
    x.p = rgamma(1,gama,lambda)
    razao.num = f(x.p,a=a,b=b)*dgamma(x[t-1],gama,lambda)
    razao.den = f(x[t-1],a=a,b=b)*dgamma(x.p,gama,lambda)
    razao = razao.num/razao.den
    alpha = min(1,razao)
    u = runif(1,0,1)
    if (u<alpha){
      x[t] = x.p
      cont = cont+1
    }else{
      x[t] = x[t-1]
    }
  }  
  
  # grafico da cadeia
  ts.plot(x,xlab="iteracoes",ylab=expression(x))

  # tirando o burn-in
  ts.plot(x[1000:M],xlab="iteracoes",ylab=expression(x))

  # grafico da densidade
  hist(x[1000:M],xlab=expression(x),ylab=expression(f(x)),freq=FALSE,nclass=50,main="Histograma dos pontos gerados")
  curve(f(x,a=a,b=b),col=2,lwd=2,add=TRUE)

  # taxa de aceitacao
  cont/M    # a taxa de aceitacao melhorou

## [1] 0.5855

  # grafico da autocorrelacao
  acf(x,main="Gráfico de autocorrelação")

  # cadeia com lag = 30
  x.lag = x[seq(1000,M,30)]
  ts.plot(x.lag)

  # intervalo de 95%
  quantile(x,0.025)

##      2.5% 
## 0.5408682

  quantile(x,0.975)

##    97.5% 
## 10.89915

  # pacote com medidas de diagnostico
  library(coda)
  
  cadeia = mcmc(x,start=1000,thin=30)
  summary(cadeia)

## 
## Iterations = 1000:300970
## Thinning interval = 30 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean             SD       Naive SE Time-series SE 
##        2.71215        2.91375        0.02914        0.07315 
## 
## 2. Quantiles for each variable:
## 
##    2.5%     25%     50%     75%   97.5% 
##  0.5409  1.1172  1.7709  3.0960 10.8991

  geweke.diag(cadeia) # deve estar entre -1.96 e 1.96 

## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##   var1 
## -1.268

  effectiveSize(cadeia)

##     var1 
## 1586.553

  # Caso haja mais de uma cadeia
  #gelman.diag(cadeias)  

Exemplo 2: Distribuição t-Student com 3 parâmetros

Gerar valores para \(X\) considerando a distribuição t-Student não padronizada com parâmetros \(\mu\), \(\sigma^2\) e \(\nu\), cuja f.d.p. é dada por \[ f(x) = \frac{\Gamma\left(\frac{\nu+1}{2}\right)} {\Gamma\left(\frac{\nu}{2}\right)\sqrt{\pi\nu\sigma^2}} \left( 1+\frac{1}{\nu} \left(\frac{x-\mu}{\sigma}\right)^2 \right)^{-(\nu+1)/2}, \qquad x \in \mathbb{R}. \] usando o algoritmo Metropolis-Hastings.

  dt3p = function(x, mu, sigma2, nu) {
    constante = gamma((nu+1)/2)/(gamma(nu/2)*sqrt(pi*nu*sigma2))
    nucleo = ( 1 + (1/nu)*((x-mu)/sqrt(sigma2))^2 )^(-(nu+1)/2)
    densidade = constante*nucleo
    return(densidade)
  }

  mu = 5
  sigma2 = 1
  nu = 10

  curve(dt3p(x,mu=mu,sigma2=sigma2,nu=nu),-20,20,xlab=expression(x),ylab=expression(f(x)),lwd=2)

Novamente, vamos implementar o algoritmo utilizando diferentes distribuições auxiliares.

Proposta Normal

Inicialmente, utilizaremos como distribuição proposta a distribuição normal com média \(\mu\) e variância 1, ou seja, \(N(\mu,1)\).

\[ \alpha = \min\left\{ 1, \frac{f(\boldsymbol{x}^p)}{f(\boldsymbol{x}^c)}\frac{q(\boldsymbol{x}^c|\boldsymbol{x}^p)}{q(\boldsymbol{x}^p|\boldsymbol{x}^c)} \right\} = \min\left\{ 1, \frac{ \frac{\Gamma\left(\frac{\nu+1}{2}\right)} {\Gamma\left(\frac{\nu}{2}\right)\sqrt{\pi\nu\sigma^2}} \left( 1+\frac{1}{\nu} \left(\frac{x^p-\mu}{\sigma}\right)^2 \right)^{-(\nu+1)/2} }{ \frac{\Gamma\left(\frac{\nu+1}{2}\right)} {\Gamma\left(\frac{\nu}{2}\right)\sqrt{\pi\nu\sigma^2}} \left( 1+\frac{1}{\nu} \left(\frac{x^c-\mu}{\sigma}\right)^2 \right)^{-(\nu+1)/2} } \frac{ \dfrac{1}{\sqrt{2\pi}}\exp\left\{-\dfrac{1}{2}(x^c-\mu)^2\right\} }{ \dfrac{1}{\sqrt{2\pi}}\exp\left\{-\dfrac{1}{2}(x^p-\mu)^2\right\} } \right\} \]

  # para garantir reprodutibilidade
  set.seed(1)
  
  #------------------------#
  # Usando proposta normal #
  #------------------------#
  
  # criando um vetor para receber a cadeia de valores gerados
  x = NULL
  # inicializando um contador
  cont = 0
  
  # numero de iteracoes
  M = 10000
  
  # inicializando a cadeia
  x[1] = 0.1
  
  for(t in 2:M){
    x.p = rnorm(1,mu,1)
    razao.num = dt3p(x.p,mu=mu,sigma2=sigma2,nu=nu)*dnorm(x[t-1],mu,1)
    razao.den = dt3p(x[t-1],mu=mu,sigma2=sigma2,nu=nu)*dnorm(x.p,mu,1)
    razao = razao.num/razao.den
    alpha = min(1,razao)
    u = runif(1,0,1)
    if (u<alpha){
      x[t] = x.p
      cont = cont+1
    }else{
      x[t] = x[t-1]
    }
  }  
  
  # grafico da cadeia
  ts.plot(x,xlab="iteracoes",ylab=expression(x))

  # tirando o burn-in
  ts.plot(x[1000:M],xlab="iteracoes",ylab=expression(x))

  # grafico da densidade
  hist(x[1000:M],xlab=expression(x),ylab=expression(f(x)),freq=FALSE,nclass=50,main="Histograma dos pontos gerados")
  curve(dt3p(x,mu=mu,sigma2=sigma2,nu=nu),col=2,lwd=2,add=TRUE)

  # taxa de aceitacao
  cont/M    # a taxa de aceitacao melhorou

## [1] 0.9314

  # grafico da autocorrelacao
  acf(x,main="Gráfico de autocorrelação")

  # cadeia com lag = 50
  x.lag = x[seq(1000,M,50)]
  acf(x.lag,main="Gráfico de autocorrelação")

  # intervalo de 95%
  quantile(x,0.025)

## 2.5% 
##  0.1

  quantile(x,0.975)

##    97.5% 
## 7.187872

  # pacote com medidas de diagnostico
  library(coda)
  
  cadeia = mcmc(x,start=1000,thin=50)
  summary(cadeia)

## 
## Iterations = 1000:500950
## Thinning interval = 50 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean             SD       Naive SE Time-series SE 
##        4.82854        1.39414        0.01394        0.12780 
## 
## 2. Quantiles for each variable:
## 
##  2.5%   25%   50%   75% 97.5% 
## 0.100 4.215 4.942 5.666 7.188

  geweke.diag(cadeia) # deve estar entre -1.96 e 1.96 

## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##   var1 
## -1.145

  # Caso haja mais de uma cadeia
  #gelman.diag(cadeias)  
  
  # tamanho efetovo da amostra
  effectiveSize(cadeia)

##     var1 
## 118.9975

Vamos rodar novamente, agora com mais iterações.

  # para garantir reprodutibilidade
  set.seed(1)

  #------------------------#
  # Usando proposta normal #
  #------------------------#
  
  # criando um vetor para receber a cadeia de valores gerados
  x = NULL
  # inicializando um contador
  cont = 0
  
  # numero de iteracoes
  M = 100000
  
  # inicializando a cadeia
  x[1] = 0.1
  
  for(t in 2:M){
    x.p = rnorm(1,5,1)
    razao.num = dt3p(x.p,mu=mu,sigma2=sigma2,nu=nu)*dnorm(x[t-1],5,1)
    razao.den = dt3p(x[t-1],mu=mu,sigma2=sigma2,nu=nu)*dnorm(x.p,5,1)
    razao = razao.num/razao.den
    alpha = min(1,razao)
    u = runif(1,0,1)
    if (u<alpha){
      x[t] = x.p
      cont = cont+1
    }else{
      x[t] = x[t-1]
    }
  }  
  
  # grafico da cadeia
  ts.plot(x,xlab="iteracoes",ylab=expression(x))

  # tirando o burn-in
  ts.plot(x[1000:M],xlab="iteracoes",ylab=expression(x))

  # grafico da densidade
  hist(x[1000:M],xlab=expression(x),ylab=expression(f(x)),freq=FALSE,nclass=50,main="Histograma dos pontos gerados")
  curve(dt3p(x,mu=mu,sigma2=sigma2,nu=nu),col=2,lwd=2,add=TRUE)

  # taxa de aceitacao
  cont/M    # a taxa de aceitacao melhorou

## [1] 0.9608

  # grafico da autocorrelacao
  acf(x,main="Gráfico de autocorrelação")

  # cadeia com lag = 100
  x.lag = x[seq(1000,M,100)]
  acf(x.lag,main="Gráfico de autocorrelação")

  # intervalo de 95%
  quantile(x,0.025)

##    2.5% 
## 2.70069

  quantile(x,0.975)

##    97.5% 
## 7.187745

  # pacote com medidas de diagnostico
  library(coda)
  
  cadeia = mcmc(x,start=1000,thin=100)
  summary(cadeia)

## 
## Iterations = 1000:10000900
## Thinning interval = 100 
## Number of chains = 1 
## Sample size per chain = 1e+05 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean             SD       Naive SE Time-series SE 
##       4.977254       1.132451       0.003581       0.013105 
## 
## 2. Quantiles for each variable:
## 
##  2.5%   25%   50%   75% 97.5% 
## 2.701 4.290 4.997 5.693 7.188

  # tamanho efetivo da amostra
  effectiveSize(cadeia)

##    var1 
## 7467.42

Proposta Passeio Aleatório

Agora utilizaremos como distribuição proposta um passeio aleatório; mais especificamente, uma distribuição normal centrada no ponto corrente e variabilidade fixa, da forma \(N(x^c,\tau^2)\).

Note que neste caso a razão de aceitação se simplifica devido a simetria da distribuição normal em torno da média:

\[ \alpha = \min\left\{ 1, \frac{f(\boldsymbol{x}^p)}{f(\boldsymbol{x}^c)}\frac{q(\boldsymbol{x}^c|\boldsymbol{x}^p)}{q(\boldsymbol{x}^p|\boldsymbol{x}^c)} \right\} = \min\left\{ 1, \frac{ \frac{\Gamma\left(\frac{\nu+1}{2}\right)} {\Gamma\left(\frac{\nu}{2}\right)\sqrt{\pi\nu\sigma^2}} \left( 1+\frac{1}{\nu} \left(\frac{x^p-\mu}{\sigma}\right)^2 \right)^{-(\nu+1)/2} }{ \frac{\Gamma\left(\frac{\nu+1}{2}\right)} {\Gamma\left(\frac{\nu}{2}\right)\sqrt{\pi\nu\sigma^2}} \left( 1+\frac{1}{\nu} \left(\frac{x^c-\mu}{\sigma}\right)^2 \right)^{-(\nu+1)/2} } \frac{ \dfrac{1}{\sqrt{2\pi\tau^2}}\exp\left\{-\dfrac{1}{2\tau^2}(x^c-x^p)^2\right\} }{ \dfrac{1}{\sqrt{2\pi\tau^2}}\exp\left\{-\dfrac{1}{2\tau^2}(x^p-x^c)^2\right\} } \right\} \]

  # para garantir reprodutibilidade
  set.seed(1)
  
  #--------------------------#
  # Usando passeio aleatório #
  #--------------------------#
  
  # variância da distribuição proposta
  tau2 = 1
  
  # criando um vetor para receber a cadeia de valores gerados
  x = NULL
  # inicializando um contador
  cont = 0
  
  # numero de iteracoes
  M = 10000
  
  # inicializando a cadeia
  x[1] = rnorm(1)
  
  for(t in 2:M){
    x.p = rnorm(1,x[t-1],sqrt(tau2))
    razao.num = dt3p(x.p,mu=mu,sigma2=sigma2,nu=nu)
    razao.den = dt3p(x[t-1],mu=mu,sigma2=sigma2,nu=nu)
    razao = razao.num/razao.den
    alpha = min(1,razao)
    u = runif(1,0,1)
    if (u<alpha){
      x[t] = x.p
      cont = cont+1
    }else{
      x[t] = x[t-1]
    }
  }  
  
  # grafico da cadeia
  ts.plot(x,xlab="iteracoes",ylab=expression(x))

  # tirando o burn-in
  ts.plot(x[1000:M],xlab="iteracoes",ylab=expression(x))

  # grafico da densidade
  hist(x[1000:M],xlab=expression(x),ylab=expression(f(x)),freq=FALSE,nclass=50,main="Histograma dos pontos gerados")
  curve(dt3p(x,mu=mu,sigma2=sigma2,nu=nu),col=2,lwd=2,add=TRUE)

  # taxa de aceitacao
  cont/M    

## [1] 0.7185

  # grafico da autocorrelacao
  acf(x,main="Gráfico de autocorrelação")

  # cadeia com lag = 40
  x.lag = x[seq(1000,M,40)]
  acf(x.lag,main="Gráfico de autocorrelação")

  # intervalo de 95%
  quantile(x,0.025)

##    2.5% 
## 2.51416

  quantile(x,0.975)

##    97.5% 
## 7.308525

  # pacote com medidas de diagnostico
  library(coda)
  
  cadeia = mcmc(x,start=1000,thin=100)
  summary(cadeia)

## 
## Iterations = 1000:1000900
## Thinning interval = 100 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean             SD       Naive SE Time-series SE 
##         4.9751         1.1703         0.0117         0.0425 
## 
## 2. Quantiles for each variable:
## 
##  2.5%   25%   50%   75% 97.5% 
## 2.514 4.284 4.981 5.682 7.309

  geweke.diag(cadeia) # deve estar entre -1.96 e 1.96 

## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##    var1 
## -0.7899

  effectiveSize(cadeia)

##     var1 
## 758.3629

Exemplo 3: Verificando o efeito da variância do passeio aleatório

Gerar valores para \(X\) considerando a distribuição t-Student padrão (caso particular em que \(\mu=0\) e \(\sigma^2=1\)) usando o algoritmo Metropolis-Hastings.

  dt3p = function(x, mu, sigma2, nu) {
    constante = gamma((nu+1)/2)/(gamma(nu/2)*sqrt(pi*nu*sigma2))
    nucleo = ( 1 + (1/nu)*((x-mu)/sqrt(sigma2))^2 )^(-(nu+1)/2)
    densidade = constante*nucleo
    return(densidade)
  }

  mu = 0
  sigma2 = 1
  nu = 5

  curve(dt3p(x,mu=mu,sigma2=sigma2,nu=nu),-10,10,xlab=expression(x),ylab=expression(f(x)),lwd=2)
  curve(dt(x,df=nu),col=2,lwd=2,lty=2,add=TRUE)

Vamos implementar o algoritmo utilizando como proposta um passeio aleatório com diferentes valores para a variância.

Proposta Passeio Aleatório com variância pequena

Agora utilizaremos como distribuição proposta normal centrada no ponto corrente e variabilidade fixa, da forma \(N(x^c,\tau^2)\).

Note que neste caso a razão de aceitação se simplifica devido a simetria da distribuição normal em torno da média:

\[ \alpha = \min\left\{ 1, \frac{f(\boldsymbol{x}^p)}{f(\boldsymbol{x}^c)}\frac{q(\boldsymbol{x}^c|\boldsymbol{x}^p)}{q(\boldsymbol{x}^p|\boldsymbol{x}^c)} \right\} = \min\left\{ 1, \frac{ \dfrac{b^a}{\Gamma(a)}(x^p)^{-a-1}\exp\left\{-\frac{b}{x^p}\right\} }{ \dfrac{b^a}{\Gamma(a)}(x^c)^{-a-1}\exp\left\{-\frac{b}{x^c}\right\} } \frac{ \dfrac{1}{\sqrt{2\pi\tau^2}}\exp\left\{-\dfrac{1}{2\tau^2}(x_c-x_p)^2\right\} }{ \dfrac{1}{\sqrt{2\pi\tau^2}}\exp\left\{-\dfrac{1}{2\tau^2}(x_p-x_c)^2\right\} } \right\} \]

Fixando \(\tau^2 = 1\)

  # para garantir reprodutibilidade
  set.seed(1)
  
  #--------------------------#
  # Usando passeio aleatório #
  #--------------------------#
  
  # variância da distribuição proposta
  tau2 = 1
  
  # criando um vetor para receber a cadeia de valores gerados
  x = NULL
  # inicializando um contador
  cont = 0
  
  # numero de iteracoes
  M = 10000
  
  # inicializando a cadeia
  x[1] = rnorm(1)
  
  for(t in 2:M){
    x.p = rnorm(1,x[t-1],sqrt(tau2))
    razao.num = dt3p(x.p,mu=mu,sigma2=sigma2,nu=nu)
    razao.den = dt3p(x[t-1],mu=mu,sigma2=sigma2,nu=nu)
    razao = razao.num/razao.den
    alpha = min(1,razao)
    u = runif(1,0,1)
    if (u<alpha){
      x[t] = x.p
      cont = cont+1
    }else{
      x[t] = x[t-1]
    }
  }  
  
  # grafico da cadeia
  ts.plot(x,xlab="iteracoes",ylab=expression(x))

  # tirando o burn-in
  ts.plot(x[1000:M],xlab="iteracoes",ylab=expression(x))

  # grafico da densidade
  hist(x[1000:M],xlab=expression(x),ylab=expression(f(x)),freq=FALSE,nclass=50,main="Histograma dos pontos gerados")
  curve(dt3p(x,mu=mu,sigma2=sigma2,nu=nu),col=2,lwd=2,add=TRUE)

  # taxa de aceitacao
  cont/M    

## [1] 0.7271

  # grafico da autocorrelacao
  acf(x,main="Gráfico de autocorrelação")

  # cadeia com lag = 40
  x.lag = x[seq(1000,M,40)]
  acf(x.lag,main="Gráfico de autocorrelação")

  # intervalo de 95%
  quantile(x,0.025)

##     2.5% 
## -2.78516

  quantile(x,0.975)

##    97.5% 
## 2.692719

  # pacote com medidas de diagnostico
  library(coda)
  
  cadeia = mcmc(x,start=1000,thin=100)
  summary(cadeia)

## 
## Iterations = 1000:1000900
## Thinning interval = 100 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean             SD       Naive SE Time-series SE 
##      -0.008242       1.419336       0.014193       0.068879 
## 
## 2. Quantiles for each variable:
## 
##     2.5%      25%      50%      75%    97.5% 
## -2.78516 -0.74231 -0.02009  0.68442  2.69272

  geweke.diag(cadeia) # deve estar entre -1.96 e 1.96 

## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##    var1 
## -0.8377

  effectiveSize(cadeia)

##     var1 
## 424.6107

Proposta Passeio Aleatório com variância grande

Modificando a variância do passeio aleatório:

Agora vamos fixar \(\tau^2 = 10\)

  # para garantir reprodutibilidade
  set.seed(1)

  #--------------------------#
  # Usando passeio aleatório #
  #--------------------------#
  
  # variância da distribuição proposta
  tau2 = 10
  
  # criando um vetor para receber a cadeia de valores gerados
  x = NULL
  # inicializando um contador
  cont = 0
  
  # numero de iteracoes
  M = 10000
  
  # inicializando a cadeia
  x[1] = rnorm(1)
  
  for(t in 2:M){
    x.p = rnorm(1,x[t-1],sqrt(tau2))
    razao.num = dt3p(x.p,mu=mu,sigma2=sigma2,nu=nu)
    razao.den = dt3p(x[t-1],mu=mu,sigma2=sigma2,nu=nu)
    razao = razao.num/razao.den
    alpha = min(1,razao)
    u = runif(1,0,1)
    if (u<alpha){
      x[t] = x.p
      cont = cont+1
    }else{
      x[t] = x[t-1]
    }
  }  
  
  # grafico da cadeia
  ts.plot(x,xlab="iteracoes",ylab=expression(x))

  # tirando o burn-in
  ts.plot(x[1000:M],xlab="iteracoes",ylab=expression(x))

  # grafico da densidade
  hist(x[1000:M],xlab=expression(x),ylab=expression(f(x)),freq=FALSE,main="Histograma dos pontos gerados")
  curve(dt3p(x,mu=mu,sigma2=sigma2,nu=nu),col=2,lwd=2,add=TRUE)

  # taxa de aceitacao
  cont/M    

## [1] 0.3907

  # grafico da autocorrelacao
  acf(x,main="Gráfico de autocorrelação")

  # cadeia com lag = 40
  x.lag = x[seq(1000,M,40)]
  acf(x.lag,main="Gráfico de autocorrelação")

  # intervalo de 95%
  quantile(x,0.025)

##      2.5% 
## -2.702669

  quantile(x,0.975)

##    97.5% 
## 2.489775

  # pacote com medidas de diagnostico
  library(coda)
  
  cadeia = mcmc(x,start=1000,thin=100)
  summary(cadeia)

## 
## Iterations = 1000:1000900
## Thinning interval = 100 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean             SD       Naive SE Time-series SE 
##      -0.001628       1.256091       0.012561       0.027975 
## 
## 2. Quantiles for each variable:
## 
##     2.5%      25%      50%      75%    97.5% 
## -2.70267 -0.67334  0.03337  0.70827  2.48977

  geweke.diag(cadeia) # deve estar entre -1.96 e 1.96 

## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##    var1 
## -0.8553

  effectiveSize(cadeia)

##     var1 
## 2016.031

**Note que o valor de $\tau^2$ tem grande impacto na taxa de aceitação do algoritmo**. 

---

Outras observações

1. Distribuições propostas do tipo \(q(.|\boldsymbol{x}^{(t-1)})\) que geram valores tais que \[\boldsymbol{x}^{p} = \boldsymbol{x}^{(t-1)} + \epsilon,\] onde \(\epsilon\) é uma perturbação aleatória simétrica em torno de 0, são conhecidas como passeios aleatórios.

2. Na prática, é bem comum utilizar passeios aleatórios supondo \(\epsilon \sim N(0,\tau^2)\); sendo necessário ajustar a variância da distribuição:

• se a variância for pequena, aceitamos mais movimentos, porém os movimentos são lentos no domínio das variáveis;

• por outro lado, se a variância for grande, aceitamos menos movimentos, porém visitamos mais rapidamente o dominio de interesse.

---

Exemplo 4: Mistura de Normais Bivariadas

Gerar uma amostra da mistura de normais bivariadas, cuja função de densidade de probabilidade conjunta é dada por \[f(\boldsymbol{x}) = 0,4\;N_2(\mu_1,\Sigma_1) + 0,6\;N_2(\mu_2,\Sigma_2),\] com \[ \mu_1 = \left(\begin{array}{c}0\\6\end{array}\right), \quad \Sigma_1 = \left(\begin{array}{cc}4&0,5\\0,5&2\end{array}\right),\quad \mu_2 = \left(\begin{array}{c}-5\\4\end{array}\right), \quad \Sigma_2 = \left(\begin{array}{cc}1&-0,2\\-0,2&1\end{array}\right). \]

Vamos considerar como densidade auxiliar a distribuição Normal Bivariada centrada no ponto anterior e com matriz de covariâncias igual a matriz identidade.

Note que pela simetria da distribuição normal em torno do vetor de médias, a razão de aceitação será simplificada.

\[ \alpha = \min\left\{ 1, \frac{f(\boldsymbol{x}^p)}{f(\boldsymbol{x}^{(t-1)})}\right\} \]

  # para garantir reprodutibilidade
  set.seed(1)
  
  #------------------------------------------#
  # Geracao de Mistura de Normais Bivariadas #
  #------------------------------------------#
    
  library(mvtnorm)

  # Exemplo de normal bivariada
  x = seq(-10, 10, length= 100)
  y = x
  mu = c(0,5)
  Sigma = matrix(c(4,0.2,0.2,1),2,2)
  zfun  = function(x,mu,Sigma) { 
    x1 = x[1]
    x2 = x[2]
    dmvnorm(c(x1,x2),mean=mu,sigma=Sigma) }
  nx = length(x)
  ny = length(y)
  z = matrix(NA,nx,ny)
  for(i in 1:nx)
  {
    for(j in 1:ny)
    {
      z[i,j] = zfun(c(x[i],y[j]),mu,Sigma)
    }
  }  
  persp(x,y,z,theta = 45, phi = 30, expand = 1, col = "lightblue")

  contour(x,y,z)

  # Mistura de normais bivariadas
  x = seq(-10, 10, length= 100)
  y = x
  mu1 = c(0,6)
  Sigma1 = matrix(c(4,0.5,0.5,2),2,2)
  mu2 = c(-5,4)
  Sigma2 = matrix(c(1,-0.2,-0.2,1),2,2)  
  pi = 0.4
  zfun  = function(vec,mu1,Sigma1,mu2,Sigma2) { 
    x = vec[1]
    y = vec[2]
    pi*dmvnorm(c(x,y),mean=mu1,sigma=Sigma1) + (1-pi)*dmvnorm(c(x,y),mean=mu2,sigma=Sigma2) }
  nx = length(x)
  ny = length(y)
  z = matrix(NA,nx,ny)
  for(i in 1:nx)
  {
    for(j in 1:ny)
    {
      z[i,j] = zfun(c(x[i],y[j]),mu1,Sigma1,mu2,Sigma2)
    }
  }  
  persp(x,y,z,theta = 215, phi = 30, expand = 1, col = "lightblue")

  contour(x,y,z)

  z = z*100
  zlim = range(z)
  zlen = zlim[2] - zlim[1] + 1
  colorlut = topo.colors(zlen) 
  col = colorlut[ z-zlim[1]+1 ]
  require(rgl)

## Carregando pacotes exigidos: rgl

  open3d()

## wgl 
##   1

  surface3d(x, y, z, color=col, back="lines") 
  
  # Inicializando o Metropolis
  
  cont = 0
  
  # numero de iteracoes
  M = 100000
  
  # inicializando a cadeia
  x.vec = matrix(NA,2,M)
  x.vec[,1] = c(0,0)
  sigma = 1
  
  for(t in 2:M){
    x.vec.p = rmvnorm(1,x.vec[,t-1],diag(sigma,2))
    razao = zfun(x.vec.p,mu1,Sigma1,mu2,Sigma2)/zfun(x.vec[,t-1],mu1,Sigma1,mu2,Sigma2)
    alpha = min(1,razao)
    u = runif(1,0,1)
    if (u<alpha){
        x.vec[,t] = x.vec.p
        cont = cont+1
    }else{
        x.vec[,t] = x.vec[,t-1]
    }
  }
  
  # graficos das cadeias
  ts.plot(x.vec[1,],xlab="iteracoes",ylab=expression(x[1]))

  ts.plot(x.vec[2,],xlab="iteracoes",ylab=expression(x[2]))

  # tirando o burn-in
  ts.plot(x.vec[1,10000:M],xlab="iteracoes",ylab=expression(x[1]))

  ts.plot(x.vec[2,10000:M],xlab="iteracoes",ylab=expression(x[2]))

  # graficos das densidades marginais
  hist(x.vec[1,10000:M],freq=F,xlab=expression(x[1]),ylab=expression(f[X[1]](x[1])),main="Histograma dos pontos gerados")

  hist(x.vec[2,10000:M],freq=F,xlab=expression(x[2]),ylab=expression(f[X[2]](x[2])),main="Histograma dos pontos gerados")

  # taxa de aceitacao
  cont/M 

## [1] 0.63044

  # graficos da autocorrelacao
  acf(x.vec[1,],main="Gráfico de autocorrelação")

  acf(x.vec[2,],main="Gráfico de autocorrelação")

  # cadeia com lag = 40
  x.vec1.lag = x.vec[1,seq(10000,M,100)]
  acf(x.vec1.lag,main="Gráfico de autocorrelação")

  ts.plot(x.vec1.lag,xlab="iteracoes",ylab=expression(x[1]))

  x.vec2.lag = x.vec[2,seq(10000,M,100)]
  acf(x.vec2.lag,main="Gráfico de autocorrelação")

  ts.plot(x.vec2.lag,xlab="iteracoes",ylab=expression(x[2]))  

  # pacote com medidas de diagnostico
  library(coda)
  
  cadeia = mcmc(t(x.vec),start=10000,thin=100)
  summary(cadeia)

## 
## Iterations = 10000:10009900
## Thinning interval = 100 
## Number of chains = 1 
## Sample size per chain = 1e+05 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##        Mean    SD Naive SE Time-series SE
## [1,] -3.022 2.853 0.009023        0.08657
## [2,]  4.801 1.535 0.004854        0.03293
## 
## 2. Quantiles for each variable:
## 
##        2.5%    25%    50%     75% 97.5%
## var1 -6.730 -5.220 -4.111 -0.6913 3.031
## var2  2.267  3.697  4.588  5.7591 8.177

  geweke.diag(cadeia)  

## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##   var1   var2 
## 0.6935 0.4801

---

Algoritmo Amostrador de Gibbs (Gibbs Sampler)

Suponha \(\boldsymbol{X} = (\boldsymbol{X}_1,\boldsymbol{X}_2,\cdots,\boldsymbol{X}_{i-1},\boldsymbol{X}_i,\boldsymbol{X}_{i+1},\cdots,\boldsymbol{X}_k)'\) um vetor aleatório.

A ideia do amostrador de Gibbs é dividir o vetor \(\boldsymbol{X}\) de forma que cada pedaço de \(\boldsymbol{X}\) (\(\boldsymbol{X}_i\) uni ou multidimensional) possa ser gerado de sua \(^{\star}\) \(f(\boldsymbol{x}_i|\boldsymbol{x}_{-i})\), \(\boldsymbol{x}_{-i}=(\boldsymbol{x}_1,\boldsymbol{x}_2,\cdots,\boldsymbol{x}_{i-1},\boldsymbol{x}_{i+1},\cdots,\boldsymbol{x}_k)'\).

\(^{\star}\) a distribuição condicional completa é a distribuição da \(i\)-ésima componente de \(\boldsymbol{X}\) condicionada em todas as outras componentes de \(\boldsymbol{X}\), exceto \(\boldsymbol{X}_i\).

Exemplo

Suponha um vetor aleatório contínuo \((X,Y,Z,W)\), de forma que sua função de densidade de probabilidade conjunta possa ser decomposta da seguinte forma (supondo as devidas independências condicionais):

\[p(x,y,z,w) = p(x|y,z,w)p(y|z,w)p(z|w)p(w) = p(x|y,w)p(y|z)p(z)p(w)\] Neste caso, as distribuições condicionais completas da componente do vetor serão dadas por

• \(f(x|y,z,w) \propto p(x|y,w)\)
• \(f(y|x,z,w) \propto p(x|y,w)p(y|z)\)
• \(f(z|x,y,w) \propto p(y|z)p(z)\)
• \(f(w|x,y,z) \propto p(x|y,w)p(w)\)

Algoritmo

Supondo que sabemos amostrar de cada distribuição condicional completa \(f(\boldsymbol{x}_i|\boldsymbol{x}_{-i})\), \(i=1,\cdots,k\), o algoritmo pode ser descrito pelos seguintes passos.

1. Inicialize a cadeia com um valor inicial \(\boldsymbol{x}^{(0)}\) e inicialize um contador \(t=1\);

2. para \(i=1,...,k\), gere um novo valor \(\boldsymbol{x}_i^{(t)} \sim f(\boldsymbol{x}_i|\boldsymbol{x}_{-i}^{(t-1)})\), \(\boldsymbol{x}_{-i}^{(t)}=(\boldsymbol{x}_1^t,\boldsymbol{x}_2^t,\cdots,\boldsymbol{x}_{i-1}^t,\boldsymbol{x}_{i+1}^{(t-1)},...,\boldsymbol{x}_k^{(t-1)})'\);

3. faça \(t=t+1\) e volte ao passo \(2\) até atingir a convergência.

Os valores simulados para a cadeia após a convergência representarão uma amostra da distribuição conjunta de interesse \(f(\boldsymbol{x})\).

Exemplo 5: Normal Bivariada

Suponha que queremos gerar valores de um vetor aleatório \((X,Y)\) com distribuição normal bivariada, cuja função de densidade conjunta é dada por

\[ \begin{pmatrix} X \\ Y \end{pmatrix} \sim N \left( \begin{pmatrix} 0 \\ 0 \end{pmatrix}, \begin{pmatrix} 1 & \rho \\ \rho & 1 \end{pmatrix} \right) \]

  library(mvtnorm)

  # visualização da distribuição de interesse

  # parâmetro de correlação
  rho = -0.8

  x = seq(-5, 5, length= 100)
  y = x
  mu = c(0,0)
  Sigma = matrix(c(1,rho,rho,1),2,2)
  zfun  = function(x,mu,Sigma) { 
    x1 = x[1]
    x2 = x[2]
    dmvnorm(c(x1,x2),mean=mu,sigma=Sigma) }
  nx = length(x)
  ny = length(y)
  z = matrix(NA,nx,ny)
  for(i in 1:nx)
  {
    for(j in 1:ny)
    {
      z[i,j] = zfun(c(x[i],y[j]),mu,Sigma)
    }
  }  
  persp(x,y,z,theta = 170, phi = 30, expand = 1, col = "lightblue")

  contour(x,y,z)

Relembrando um resultado visto em Probabilidade 2:

Sejam \(\boldsymbol{X}\) e \(\boldsymbol{Y}\) dois vetores aleatórios de dimensão \(p\) e \(q\), respectivamente, tais que \[ \displaystyle \left(\begin{array}{c} \boldsymbol{X} \\ \boldsymbol{Y} \end{array}\right) \sim N_{p+q}\left( \boldsymbol{\mu} = \left(\;\begin{array}{c} \boldsymbol{\mu}_X \\ \boldsymbol{\mu}_Y \end{array}\right)\;,\; \boldsymbol{\Sigma} = \left(\begin{array}{cc} \boldsymbol{\Sigma}_X & \boldsymbol{\Sigma}_{XY} \\ \boldsymbol{\Sigma}_{YX} & \boldsymbol{\Sigma}_Y \end{array}\right) \;\right). \] Neste caso, \[{{\boldsymbol{X}|\boldsymbol{Y}=\boldsymbol{y} \sim N_p\left(\boldsymbol{\mu}_{X|Y},\boldsymbol{\Sigma}_{X|Y}\right)}},\] \[ \displaystyle \boldsymbol{\mu}_{X|Y} = \boldsymbol{\mu}_X+\boldsymbol{\Sigma}_{XY}\boldsymbol{\Sigma}_Y^{-1}(\boldsymbol{y}-\boldsymbol{\mu}_Y) \qquad \qquad \mbox{e} \qquad \qquad \boldsymbol{\Sigma}_{X|Y} = \boldsymbol{\Sigma}_X-\boldsymbol{\Sigma}_{XY}\boldsymbol{\Sigma}_Y^{-1}\boldsymbol{\Sigma}_{YX}. \]

Utilizando as propriedades da distribuição normal multivariada, obtemos as seguintes distribuições condicionais completas exatas:

\[X \mid Y=y \sim N(\rho y, 1 - \rho^2) \qquad \text{e} \qquad Y \mid X=x \sim N(\rho x, 1 - \rho^2)\]

  # para garantir reprodutibilidade
  set.seed(1)
  
  # numero de iteracoes
  M = 10000
  
  # criando os vetores para receber as cadeias de valores gerados
  x = NULL
  y = NULL
  
  # inicializando a cadeia
  x[1] = 10
  y[1] = -5
  
  # parâmetro de correlação
  rho = -0.8
  
  # desvio padrão condicional
  sd_cond = sqrt(1 - rho^2)
  
  for(t in 2:M){
    # 1. atualizando X dado o Y anterior
    x[t] = rnorm(1, mean = rho * y[t-1], sd = sd_cond)
    
    # 2. atualizando Y dado o X recém-atualizado
    y[t] = rnorm(1, mean = rho * x[t], sd = sd_cond)
  }
  
  # graficos das cadeias
  ts.plot(x,xlab="iteracoes",ylab=expression(x))

  ts.plot(y,xlab="iteracoes",ylab=expression(y))

  # tirando o burn-in
  ts.plot(x[1000:M],xlab="iteracoes",ylab=expression(x))

  ts.plot(y[1000:M],xlab="iteracoes",ylab=expression(y))

  # graficos das densidades marginais
  hist(x[1000:M],freq=F,xlab=expression(x),ylab=expression(f[X](x)),main="Histograma dos pontos gerados")

  hist(y[1000:M],freq=F,xlab=expression(y),ylab=expression(f[Y](y)),main="Histograma dos pontos gerados")

  # graficos da autocorrelacao
  acf(x,main="Gráfico de autocorrelação")

  acf(y,main="Gráfico de autocorrelação")

  # cadeia com lag = 20
  x.lag = x[seq(1000,M,20)]
  acf(x.lag,main="Gráfico de autocorrelação")

  ts.plot(x.lag,xlab="iteracoes",ylab=expression(x))

  y.lag = y[seq(1000,M,20)]
  acf(y.lag,main="Gráfico de autocorrelação")

  ts.plot(y.lag,xlab="iteracoes",ylab=expression(y))  

  # pacote com medidas de diagnostico
  library(coda)
  
  cadeia = mcmc(cbind(x,y),start=1000,thin=20)
  summary(cadeia)

## 
## Iterations = 1000:200980
## Thinning interval = 20 
## Number of chains = 1 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##       Mean     SD Naive SE Time-series SE
## x -0.03438 1.0016 0.010016        0.02121
## y  0.03161 0.9956 0.009956        0.02108
## 
## 2. Quantiles for each variable:
## 
##     2.5%     25%      50%    75% 97.5%
## x -1.951 -0.7093 -0.05067 0.6502 1.908
## y -1.913 -0.6447  0.02756 0.7176 1.946

  geweke.diag(cadeia)    

## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##       x       y 
## -0.9497  0.9854

  library(ggplot2)
  library(tibble)

  amostras_normal = tibble(iteracao = 1:M, x = x, y = y)
  # Gráfico de dispersão mostrando o caminho inicial (Trace de convergência)
  ggplot(amostras_normal, aes(x = x, y = y)) +
    geom_point(alpha = 0.3, color = "darkblue") +
    geom_path(data = head(amostras_normal, 50), aes(color = iteracao), size = 0.8) +
    scale_color_gradient(low = "red", high = "yellow") +
    labs(title = "Cadeia do Gibbs para Normal Bivariada",
         subtitle = "Linha colorida mostra os primeiros 50 passos do algoritmo",
         x = "Variavel X", y = "Variavel Y") +
    theme_minimal()

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Observação: Algoritmo de Gibbs com passos de Metropolis-Hastings

Se \(\boldsymbol{X}\) é multivariado e não sabemos gerar das condicionais completas das componentes de \(\boldsymbol{X}\), podemos usar o algoritmo Metropolis-Hastings a cada \(i\) (\(i=1,\cdots,k\)) no passo 2 do algoritmo para gerar \(\boldsymbol{X}_{i}\) a cada iteração.

Exemplo 6: Revisitando o Exemplo da Mistura de Normais Bivariadas

Gerar uma amostra da mistura de normais bivariadas, cuja função de densidade de probabilidade conjunta é dada por \[f(\boldsymbol{x}) = 0,4\;N_2(\mu_1,\Sigma_1) + 0,6\;N_2(\mu_2,\Sigma_2),\] com \[ \mu_1 = \left(\begin{array}{c}0\\6\end{array}\right), \quad \Sigma_1 = \left(\begin{array}{cc}4&0,5\\0,5&2\end{array}\right),\quad \mu_2 = \left(\begin{array}{c}-5\\4\end{array}\right), \quad \Sigma_2 = \left(\begin{array}{cc}1&-0,2\\-0,2&1\end{array}\right); \]

desta vez, através do Algoritmo Amostrador de Gibbs considerando um passeio aleatório para cada componente.

Note que pela simetria da distribuição normal em torno do vetor de médias, cada razão de aceitação será simplificada.

\[ \begin{array}{lll} \alpha_1 = \min\left\{ 1, \frac{f(x_1^p,x_2^{(t-1)})}{f(x_1^{(t-1)},x_2^{(t-1)})}\right\} &\quad& \alpha_2 = \min\left\{ 1, \frac{f(x_1^{(t)},x_2^p)}{f(x_1^{(t)},x_2^{(t-1)})}\right\} \end{array} \]

  # para garantir reprodutibilidade
  set.seed(1)
  
  #------------------------------------------#
  # Geracao de Mistura de Normais Bivariadas #
  #------------------------------------------#
    
  library(mvtnorm)

  mu1 = c(0,6)
  Sigma1 = matrix(c(4,0.5,0.5,2),2,2)
  mu2 = c(-5,4)
  Sigma2 = matrix(c(1,-0.2,-0.2,1),2,2)  
  pi = 0.4
  zfun  = function(vec,mu1,Sigma1,mu2,Sigma2) { 
    x = vec[1]
    y = vec[2]
    pi*dmvnorm(c(x,y),mean=mu1,sigma=Sigma1) + (1-pi)*dmvnorm(c(x,y),mean=mu2,sigma=Sigma2) }

  x = seq(-10,5,by=0.1)
  y = seq(0,10,by=0.1)

  nx = length(x)
  ny = length(y)
  z = matrix(NA,nx,ny)
  for(i in 1:nx){
    for(j in 1:ny){
      z[i,j] = zfun(c(x[i],y[j]),mu1,Sigma1,mu2,Sigma2)
    }
  }
  
  persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = "lightblue",
        ltheta = 120, shade = 0.75, ticktype = "detailed",
        xlab = "x", ylab = "y", zlab = "")

  contour(x,y,z,main="curvas de contorno")  

  # Inicializando o Metropolis
  
  cont1 = cont2 = 0
  
  # numero de iteracoes
  M = 50000
  
  # definindo os vetores que receberao as cadeias
  x.vec = matrix(NA,2,M)
  
  # inicializando as cadeias (definindo valores iniciais)    
  x.vec[,1] = c(0,0)
  
  tau2 = 1
  
  for(t in 2:M){

    # Cadeia para x.vec1
    x.vec1.p = rnorm(1,x.vec[1,t-1],sqrt(tau2))
    razao = zfun(c(x.vec1.p,x.vec[2,t-1]),mu1,Sigma1,mu2,Sigma2)/zfun(c(x.vec[1,t-1],x.vec[2,t-1]),mu1,Sigma1,mu2,Sigma2)
    alpha = min(1,razao)
    u = runif(1,0,1)
    if (u<alpha){
        x.vec[1,t] = x.vec1.p
        cont1 = cont1+1
    }else{
        x.vec[1,t] = x.vec[1,t-1]
    }

    # Cadeia para x.vec2
    x.vec2.p = rnorm(1,x.vec[2,t-1],sqrt(tau2))
    razao = zfun(c(x.vec[1,t],x.vec2.p),mu1,Sigma1,mu2,Sigma2)/zfun(c(x.vec[1,t],x.vec[2,t-1]),mu1,Sigma1,mu2,Sigma2)
    alpha = min(1,razao)
    u = runif(1,0,1)
    if (u<alpha){
      x.vec[2,t] = x.vec2.p
      cont2 = cont2+1
    }else{
      x.vec[2,t] = x.vec[2,t-1]
    }   
    
  }
  
  # graficos das cadeias
  ts.plot(x.vec[1,],xlab="iteracoes",ylab=expression(x[1]))

  ts.plot(x.vec[2,],xlab="iteracoes",ylab=expression(x[2]))

  # tirando o burn-in
  ts.plot(x.vec[1,5000:M],xlab="iteracoes",ylab=expression(x[1]))

  ts.plot(x.vec[2,5000:M],xlab="iteracoes",ylab=expression(x[2]))

  # graficos das densidades marginais
  hist(x.vec[1,5000:M],freq=F,xlab=expression(x[1]),ylab=expression(f[X[1]](x[1])),main="Histograma dos pontos gerados")

  hist(x.vec[2,5000:M],freq=F,xlab=expression(x[2]),ylab=expression(f[X[2]](x[2])),main="Histograma dos pontos gerados")

  # taxas de aceitacao
  cont1/M 

## [1] 0.783

  cont2/M 

## [1] 0.74086

  # graficos da autocorrelacao
  acf(x.vec[1,],main="Gráfico de autocorrelação")

  acf(x.vec[2,],main="Gráfico de autocorrelação")

  # cadeia com lag = 40
  x.vec1.lag = x.vec[1,seq(5000,M,100)]
  acf(x.vec1.lag,main="Gráfico de autocorrelação")

  ts.plot(x.vec1.lag,xlab="iteracoes",ylab=expression(x[1]))

  x.vec2.lag = x.vec[2,seq(5000,M,100)]
  acf(x.vec2.lag,main="Gráfico de autocorrelação")

  ts.plot(x.vec2.lag,xlab="iteracoes",ylab=expression(x[2]))  

  # pacote com medidas de diagnostico
  library(coda)
  
  cadeia = mcmc(t(x.vec),start=5000,thin=100)
  summary(cadeia)

## 
## Iterations = 5000:5004900
## Thinning interval = 100 
## Number of chains = 1 
## Sample size per chain = 50000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##        Mean    SD Naive SE Time-series SE
## [1,] -2.816 2.888 0.012917        0.11961
## [2,]  4.891 1.568 0.007014        0.04833
## 
## 2. Quantiles for each variable:
## 
##        2.5%    25%    50%     75% 97.5%
## var1 -6.709 -5.160 -3.880 -0.3184 2.999
## var2  2.296  3.748  4.684  5.9041 8.322

  geweke.diag(cadeia)  

## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
##    var1    var2 
## -0.6395 -0.7053

---

Referências

CHIB, S.; GREENBERG, E. Understanding the Metropolis-Hastings Algorithm. The American Statistician, 49, 4, p. 327-335, 1995.

GELFAND, A. E.; SMITH, A. F. M. Sampling-Based Approaches to Calculating Marginal Densities. Journal of the American Statistical Association, 85, 410, p. 398-409, 1990.

GELMAN, A.; RUBIN, D. B. Inference from Iterative Simulation Using Multiple Sequences. Staff Report, Federal Reserve Bank of Minneapolis, 148, 1991.

GEWEKE, J. Evaluating the accuracy of sampling-based approaches to the calculation of posterior moments. Journal of the American Statistical Association, 85, 410, p. 398-409, 1990.

HASTINGS, W. K. Monte Carlo Sampling Methods Using Markov Chains and Their Applications. Biometrika, 57, 1, p. 97-109, 1970.

METROPOLIS, N.; ROSENBLUTH, A. W.; ROSENBLUTH, M. N.; TELLER, A. H.; TELLER, E. Equation of State Calculations by Fast Computing Machines. The Journal of Chemical Physics, 21, 6, p. 1087-1092, 1953.

ROBERT, Christian P.; CASELLA, George. Monte Carlo Statistical Methods. 2. ed. New York: Springer, 2004. (Springer Texts in Statistics).

ROSS, Sheldon M. Simulation. 5. ed. San Diego: Academic Press, 2012.