rm(list = ls())
#Algoritmo MCMC Metropolis-Hastings

# Algoritmo
metropolis <- function(log_posterior,
                       n_iter = 10000,
                       theta0 = 1,
                       proposal_sd = 1) {
  
  theta <- numeric(n_iter)
  theta[1] <- theta0
  
  n_accept <- 0
  
  for (t in 2:n_iter) {
    
    # Proposta
    theta_star <- rnorm(
      1,
      mean = theta[t - 1],
      sd = proposal_sd
    )
    
    # Log da razão de aceitação
    if(theta_star < 0){
      theta[t] <- theta[t-1]
    }
    else{
      log_alpha <- log_posterior(theta_star) - log_posterior(theta[t - 1])
      
      # Aceitação/rejeição
      if(log(runif(1)) < min(0, log_alpha)){
        theta[t] <- theta_star
        n_accept <- n_accept + 1
      } 
      else{
        theta[t] <- theta[t - 1]
      }
    }
  }
  list(
    samples = theta,
    acceptance_rate = n_accept / (n_iter - 1)
  )
}

# Para o caso de uma priori theta ~ Gamma(a,b) com a = 1 e b=1/5, 
# e Y|theta ~ exp(theta), com dados n = 20 e sum_y = 10.5 (assim como na 
# questão 1 da Prova 1), temos uma posteriori theta|y ~ Gamma(a + n, b + sum_y):

a <- 1
b <- 1/5
n <- 20
sum_y <- 10.5

# definindo a log(posteriori)
log_posterior <- function(theta) {
  log(theta^(a+n-1)*exp(-(b+sum_y)*theta))
}

#gerando as amostras por MCMC Metropolis-Hastings
set.seed(123)

resultado <- metropolis(
  log_posterior = log_posterior,
  n_iter = 100000,
  theta0 = 2,
  proposal_sd = 1
)

# Avaliação da taxa de aceitação 
# (em teoria taxas entre 0.3 e 0.5 são consideradas ideais)
resultado$acceptance_rate

#visualização das amostras (algumas, para ter uma ideia da cadeia)
resultado$samples[1:15]

# histograma da amostra gerada por MCMC M-H
hist(resultado$samples,
     breaks = 50,
     probability = TRUE, 
     main = "Histograma posteriori",
     xlab = expression(theta))

# Curva da posteriori calculada analiticamente (para comparação com o resultado 
# do histograma)
curve(
  dgamma(x, a+n, b+sum_y),
  add = TRUE,
  lwd = 2
)

theta <- resultado$samples

# Avaliação da cadeia (observar se houve convergência 
# e avaliar possíveis diagnósticos)
plot(
  theta,
  type = "l",
  xlab = "Iteração",
  ylab = expression(theta),
  main = "Traço da cadeia"
)


# Burn-in -----------------------------------


theta_burn <- theta[-(1:5000)]

hist(theta_burn,
     breaks = 50,
     probability = TRUE, 
     main = "Histograma posteriori",
     xlab = expression(theta))

curve(
  dgamma(x, a+n, b+sum_y),
  add = TRUE,
  lwd = 2
)

# Tamanho Efetivo de Amostra (ESS)

library(coda)

# Converte o vetor de amostras para um objeto mcmc
cadeia_mcmc <- mcmc(theta_burn)

# Calcula o Tamanho Efetivo de Amostra (ESS)
ess <- effectiveSize(cadeia_mcmc)

cat("Tamanho atual da cadeia:", length(theta_burn), "\n")
cat("Tamanho Efetivo de Amostra (ESS):", round(ess, 2), "\n")


# Autocorrelação

acf(theta_burn, type = "correlation")

# Mantem como amostra apenas 1 a cada N valores, para termos uma amostra
# com menor autocorrelação (o que é chamado de thinning).
N <- 10
theta_final <- theta_burn[seq(1, length(theta_burn), by = N)]


acf(theta_final, type = "correlation")


hist(theta_final,
     breaks = 50,
     probability = TRUE, 
     main = "Histograma posteriori",
     xlab = expression(theta))

curve(
  dgamma(x, a+n, b+sum_y),
  add = TRUE,
  lwd = 2
)

# Após tratarmos a autocorrelação:

# Converte o vetor de amostras para um objeto mcmc
cadeia_mcmc <- mcmc(theta_final)

# Calcula o Tamanho Efetivo de Amostra (ESS)
ess <- effectiveSize(cadeia_mcmc)

cat("Tamanho atual da cadeia:", length(theta_final), "\n")
cat("Tamanho Efetivo de Amostra (ESS):", round(ess, 2), "\n")

# --------------------------------------------------------------------------
# Testar com um chute inicial mais distante e um proposal_sd menor,
# para mostrar como a cadeia leva um certo número de iterações para
# chegar à distribuição estacionária. 

#Testar com proposal_sd = 0.01 e theta0 = 12
#Testar com proposal_sd = 10 e theta0 = 2

resultado2 <- metropolis(
  log_posterior = log_posterior,
  n_iter = 100000,
  theta0 = 12,
  proposal_sd = 0.01
)

resultado3 <- metropolis(
  log_posterior = log_posterior,
  n_iter = 100000,
  theta0 = 2,
  proposal_sd = 10
)

hist(resultado2$samples,
     breaks = 50,
     probability = TRUE, 
     main = "Histograma posteriori",
     xlab = expression(theta))
curve(
  dgamma(x, a+n, b+sum_y),
  add = TRUE,
  lwd = 2
)

hist(resultado3$samples,
     breaks = 50,
     probability = TRUE, 
     main = "Histograma posteriori",
     xlab = expression(theta))
curve(
  dgamma(x, a+n, b+sum_y),
  add = TRUE,
  lwd = 2
)


# Comparação do traço de cadeias com diferentes theta0 e proposal_sd
x11()
par(mfrow = c(3,1))

plot(
  resultado$samples,
  type = "l",
  xlab = "Iteração",
  ylab = expression(theta),
  main = "Traço da cadeia com theta0 = 2 e sd = 1"
)

plot(
  resultado2$samples,
  type = "l",
  xlab = "Iteração",
  ylab = expression(theta),
  main = "Traço da cadeia com theta0 = 12 e sd = 0.01"
)

plot(
  resultado3$samples,
  type = "l",
  xlab = "Iteração",
  ylab = expression(theta),
  main = "Traço da cadeia com theta0 = 2 e sd = 10"
)

par(mfrow = c(1,1))
# --------------------------------------------------------------------

# Calculos de esperança, variância, moda, e probabilidades intervalares 
# por simulação e por cálculo de valores tabelados/analíticos:
 
# Esperança:
media_analitica <- (a+n)/(b+sum_y)

media_sim <- mean(theta_final)

c(Media_analitica = media_analitica, media_simulada = media_sim)

# Variância:
var_analitica <- (a+n)/(b+sum_y)^2

var_sim <- var(theta_final)

c(Variancia_analitica = var_analitica, Variancia_simulada = var_sim)

# Moda: 
moda_analitica <- (a+n-1)/(b+sum_y)

densidade <- density(theta_final)

moda_sim <- densidade$x[which.max(densidade$y)]

c(Moda_analitica = moda_analitica, Moda_simulada = moda_sim)


# Mediana:

mediana_analitica <- qgamma(0.5,shape = a+n, rate = b+sum_y)

mediana_sim <- median(theta_final)

c(Mediana_analitica = mediana_analitica, Mediana_simulada = mediana_sim)

# Probabilidade de x pertencer ao intervalo lim_inf e lim_sup

lim_inf <- 1.5
lim_sup <- 2

P_analitica <- diff(pgamma(c(lim_inf,lim_sup),shape = a+n, rate = b+sum_y)) 

P_sim <- mean(theta_final > lim_inf & theta_final < lim.sup)

c(Prob_analitica = P_analitica, Prob_simulada = P_sim)

## ------------------------------------------------------------------------------

# Anexo:

# Código com a reparametrização de theta (phi = log(theta))
# Torna o algoritmo mais eficaz, e não presimamos nos preocupar 
# com o domínio de theta (theta > 0).

metropolis_par <- function(log_posterior,
                       n_iter = 10000,
                       theta0 = 1, # Importante: theta0 deve ser > 0 se log_par = TRUE
                       proposal_sd = 1, 
                       log_par = FALSE) {
  
  phi <- numeric(n_iter)
  
  # Ponto inicial (se log_par for TRUE, inicializamos no espaço log)
  if (log_par) {
    if (theta0 <= 0) stop("theta0 deve ser estritamente positivo quando log_par = TRUE")
    phi[1] <- log(theta0)
  } else {
    phi[1] <- theta0
  }
  
  n_accept <- 0
  
  for (t in 2:n_iter) {
    # 1. Proposta na escala atual (Normal em torno de phi[t - 1])
    phi_star <- rnorm(1, mean = phi[t - 1], sd = proposal_sd)
    
    if (!log_par) {
      # Caso Padrão (Sem transformação log) 
      if (phi_star < 0) {
        log_alpha <- -Inf
        phi[t] <- phi[t - 1]
      } else {
        log_alpha <- log_posterior(phi_star) - log_posterior(phi[t - 1])
        
        if (log(runif(1)) < min(0, log_alpha)) {
          phi[t] <- phi_star
          n_accept <- n_accept + 1
        } else {
          phi[t] <- phi[t - 1]
        }
      }
      
    } else {
      # Caso Transformado (Trabalhando na escala Log)
      # Adicionamos o termo do Jacobiano (phi_star e phi[t - 1])
      log_post_star <- log_posterior(exp(phi_star)) + phi_star
      log_post_curr <- log_posterior(exp(phi[t - 1])) + phi[t - 1]
      
      log_alpha <- log_post_star - log_post_curr
      
      # Aceitação/Rejeição
      if (log(runif(1)) < min(0, log_alpha)) {
        phi[t] <- phi_star
        n_accept <- n_accept + 1
      } else {
        phi[t] <- phi[t - 1]
      }
    }
  }
  
  # Se rodou na escala log, converte o vetor final de volta para a escala original (theta)
  samples <- if (log_par) exp(phi) else phi
  
  list(
    samples = samples,
    acceptance_rate = n_accept / (n_iter - 1)
  )
}

log_posterior_estavel <- function(theta) {
  if (theta <= 0) return(-Inf)
  
  # log(priori Gamma(a, b)): (a - 1) * log(theta) - b * theta
  log_priori <- (a - 1) * log(theta) - b * theta
  
  # log(verossimilhança Exponential(theta)): n * log(theta) - sum_y * theta
  log_verossimilhanca <- n * log(theta) - sum_y * theta
  
  return(log_priori + log_verossimilhanca)
}

set.seed(123)

# Gerando as amostras por MCMC Metropolis-Hastings 
# com log_par = TRUE (reparametrizando).
resultado4 <- metropolis_par(
  log_posterior = log_posterior_estavel,
  n_iter = 100000,
  theta0 = 2,
  proposal_sd = 0.6, log_par = TRUE
)

# Os resultados serão praticamente os mesmos. A diferença está na estabilidade
# do algoritmo.