Isso é um R Markdown Notebook depositado no Bioforum. Quando executar o código, os resultados aparecerão abaixo dele!

Neste R Markdown, vamos usar o recurso do software Tracer para visualizar e realizar os diagnósticos das amostras de MCMC!

De acordo com a aula de hoje, vamos criar uma rotina em que possamos desenvolver o algoritmo de Metropolis-Hastings via MCMC:

Etapa 1: fazer as f.d.p. das distribuições normal, gama e uniforme:

\(\color{red}{\text{f.d.p. da distribuição normal:}}\) \[f(x| \mu, \sigma) = \frac{1}{\sqrt{2 \pi \sigma^2}} e^{-\frac{(x - \mu)^2}{2 \sigma^2}}\]

Função de densidade de probabilidade da distribuição normal logaritimizada

log_fdp_normal <- function(x, mu, sig){
  fdp = -0.5 * log(2 * pi * sig^2) - ((x - mu)^2) / (2 * sig^2)
    return(fdp)
} # retorna o log da verossimilhança de um valor

\(\color{red}{\text{A f.d.p. da distribuição gamma é:}}\) \[f(x| \alpha, \beta) = \frac{\beta^ \alpha}{\Gamma(\alpha)} x^{\alpha - 1} e^{-\beta x}\]

Função de densidade de probabilidade da distribuição gama logaritimizada

log_fdp_gamma <- function(x, alpha, beta){
  lik = (alpha * log(beta)) - log(gamma(alpha)) + (alpha - 1) * log(x) - beta * x
  return(lik)
}

\(\color{red}{\text{A f.d.p. da distribuição uniforme:}}\) \[f(x| min, max) = \frac{1}{max - min}\]

Função de densidade de probabilidade da distribuição uniforme logaritimizada

log_fdp_uniform <- function(x, inferior, superior){ # assumindo-se que x é um único número 
    if (x > inferior & x < superior){
      return(-log(superior - inferior))
    } else {# se x estiver fora da amplitude, será -infinito porque é log(0)
      return(-Inf)
    }
}

Etapa 2: Estimar as funções de proposição de valores

# proposal functions (symmetrical: uniforme)
sliding_window <- function(x, size){
    low <- x - 0.5 * size
    up <- x + 0.5 * size # porque é metade da distribuição para baixo e para cima
    i <- runif(n = 1, min = low, max = up)
    hastings_atual <- 0 # probabilidades de ir e voltar são as mesmas e o log10 de 1 é zero!
    return(c(i, hastings_atual))
}

Etapa 3: Oferecendo valores e calculando os prioris. Esses são os valores atuais de mu e sigma!

# população de baratas d'água:
barata_size <- c(25.0, 32.5, 43.8, 84.74, 84.88, 94.60, 96.37, 102.93, 109.11, 125.76) #dez baratas d'água
mu_baratas_atual <- mean(barata_size) # média que será o valor escolhido inicial de mu
sig_baratas_atual <- sd(barata_size) # desvio que será o valor escolhido inicial de sigma

# priori de mu:
priori_med_atual <- log_fdp_uniform(x = mu_baratas_atual, inferior = min(barata_size), superior = max(barata_size))

# priori de sigma:
shape_priori <- 2
rate_priori <- 0.0001
priori_sig_atual <- log_fdp_gamma(x = sig_baratas_atual, alpha = shape_priori, beta = rate_priori)

# soma de todos os prioris: o do valor da média e o do valor de sigma
total_priori_atual <- priori_med_atual + priori_sig_atual

Etapa 4: Verossimilhanças (atuais)

# soma de todas as verossimilhanças:
total_vero_atual <- sum(log_fdp_normal(x = barata_size, mu = mu_baratas_atual, sig = sig_baratas_atual))

Etapa 5: Probabilidades posteriores (atuais)

total_posterior_atual <- total_vero_atual + total_priori_atual

Metropolis-Hastings MCMC

# MCMC settings
n_iter  <- 1000
sampling_freq <- 10
print_freq    <- 10

##### iniciando o log file:
setwd('~/Dropbox/aulas/Bayes/aulas_&_Rmark/')
# dir.create('resultados_aula4') # cria um diretório dentro do diretório de interesse
setwd('~/Dropbox/aulas/Bayes/aulas_&_Rmark/resultados_aula4/')
logfile <- "MCMC_baratas1.txt"
cat(c("it", "post", "priori", "veros", "mu", "sig", "\n"), file = logfile, sep="\t")

#### MCMC loop #####
for (m in 1:n_iter){
  # fazer a renovação de valores com uma frequência de 50%, para mudarmos uma hora um parâmetro e outra hora outro parâmetro
  moeda <- runif(n = 1) # jogando a moeda...
  if (moeda < 0.5){
    mu_baratas_novo <- sliding_window(x = mu_baratas_atual, size = 0.5)[1]
    Hastings_novo <- sliding_window(x = mu_baratas_atual, size = 0.5)[2]
    sig_baratas_novo <- sig_baratas_atual # mantendo o mesmo valor
  } else {
    mu_baratas_novo <- mu_baratas_atual #mantendo o mesmo valor
    sig_baratas_novo <- abs(sliding_window(x = sig_baratas_atual, size = 2)[1]) # sig não pode ser negativo (reflection at the boundary)!
    Hastings_novo <- sliding_window(x = sig_baratas_atual, size = 2)[2]
  }
  
  ## VEROSSIMILHANÇA E PRIORIS NOVOS:
  # soma de todas as verossimilhanças:
  total_vero_novo <- sum(log_fdp_normal(x = barata_size, mu = mu_baratas_novo, sig = sig_baratas_novo))
  
  # soma de todos os prioris: o do valor da média e o do valor de sigma (novas propostas)
  priori_med_novo <- log_fdp_uniform(x = mu_baratas_novo, inferior = min(barata_size), superior = max(barata_size)) # mu novo
  priori_sig_novo <- log_fdp_gamma(x = sig_baratas_novo, alpha = shape_priori, beta = rate_priori) # sigma novo
  
  # priori total novo:
  total_priori_novo <- priori_med_novo + priori_sig_novo
  
  ## Acceptance ratio (r):
  # posterior ratio = likelihood ratio (LR) * prior ratio (PR) * hastings ratio (HR)
  # Lembre-se que o log de uma razão é a diferença entre os logs de cada termo!
  # Lembre-se que o log de uma multiplicação é a soma dos logs de cada termo!
  
  LR <- total_vero_novo - total_vero_atual # likelihood ratio
  PR <- total_priori_novo - total_priori_atual
  
  erre <- LR + PR + 0 # lembre-se que o Hastings ratio aqui é nulo porque o log (1) = 0!
  
  # aceitar ou rejeitar? 
  u <- runif(1) # simulando a escolha de um valor entre zero e um, se a proposta não for melhor que a atual!
  if (log(u) <= erre){ #aceitando...
    mu_baratas_atual <- mu_baratas_novo
    sig_baratas_atual <- sig_baratas_novo
    total_vero_atual <- total_vero_novo
    total_priori_atual <- total_priori_novo
  }
  
  # print to screen
  if (m %% print_freq == 0){
    print(c(m, total_vero_atual, total_priori_atual, mu_baratas_atual, sig_baratas_atual))
  }
  
  # save to file
    if (m %% sampling_freq == 0){
      cat(c(m, total_vero_atual + total_priori_atual, total_priori_atual, total_vero_atual, mu_baratas_atual,
            sig_baratas_atual, "\n"), sep = "\t", file = logfile, append = T)
    }
}

Aula 4: estimativas de probabilidade posterior usando Metropolis-Hastings MCMC

José Ricardo Inacio Ribeiro

05/11/2020