estrutura em R Markdown para explicar o Amostrador de Gibbs passo a passo. Vamos detalhar cada etapa, incluindo a formulação do algoritmo, o cálculo das médias e variâncias para cada parâmetro, e a visualização das cadeias amostradas.

Introdução

Neste documento, vamos realizar uma simulação usando o Amostrador de Gibbs para estimar os parâmetros de um modelo de regressão simples com um termo de erro. O objetivo é estimar os valores dos parâmetros \(\beta_0\), \(\beta_1\), e \(\phi\), utilizando uma abordagem bayesiana.

Simulação dos Dados

Primeiro, simulamos os dados para o modelo, assumindo valores conhecidos para \(\beta_0\), \(\beta_1\), e \(\phi\) (considerado como precisão).

# Tamanho da amostra
n <- 1000  

# Verdadeiro valor dos parâmetros
beta0.real <- 7 
beta1.real <- 1
phi.real <- 2

# Simulação da covariável X
x <- rnorm(n, 0, 1)

# Cálculo da média de Y
mu.real <- beta0.real + (beta1.real * x)

# Simulação dos valores de Y
y <- rnorm(n, mu.real, sqrt(1 / phi.real))

# Histograma de Y
hist(y, col="darkred", main="", xlab=expression(y), ylab="Frequência")

Atribuição dos Hiperparâmetros Escolhemos distribuições a priori para cada parâmetro:

# Hiperparâmetros das distribuições a priori

# Beta_0
a0 <- 6
b0 <- 8

# Beta_1
a1 <- 2
b1 <- 5

# Phi
c <- 6
d <- 3

# Visualização das priors
par(mfrow = c(1, 3))
plot(function(x) dnorm(x, a0, sqrt(b0)), -30, 30, xlab="", ylab=expression(beta[0]), col="darkred")
plot(function(x) dnorm(x, a1, sqrt(b1)), -30, 30, xlab="", ylab=expression(beta[1]), col="darkred")
plot(function(x) dgamma(x, c, d), 0, 10, xlab="", ylab=expression(phi), col="darkred")

Valores iniciais

beta0 <- NULL
beta1 <- NULL
phi <- NULL

beta0[1] <- 5
beta1[1] <- 2
phi[1] <- 4
# Número de iterações
ite <- 1000 
n <- length(y)
for(i in 2:ite) {
  ### Amostrando \(\beta_0\)
  
  # Média e variância para a distribuição condicional de \(\beta_0\)
  media <- ((b0 * sum(y)) - (b0 * beta1[i-1] * sum(x)) + (a0 / phi[i-1])) / ((n * b0) + (1 / phi[i-1]))
  var <- (b0 / phi[i-1]) / ((n * b0) + (1 / phi[i-1]))
  
  # Amostra de \(\beta_0\)
  beta0[i] <- rnorm(1, media, sqrt(var))
  
  ### Amostrando \(\beta_1\)
  
  # Média e variância para a distribuição condicional de \(\beta_1\)
  media1 <- ((b1 * sum(x * y)) - (b1 * beta0[i] * sum(x)) + (a1 / phi[i-1])) / ((sum(x^2) * b1) + (1 / phi[i-1]))
  var1 <- (b1 / phi[i-1]) / (sum(x^2) * b1 + (1 / phi[i-1]))
  
  # Amostra de \(\beta_1\)
  beta1[i] <- rnorm(1, media1, sqrt(var1))
  
  ### Amostrando \(\phi\)
  
  # Parâmetros para a distribuição condicional de \(\phi\)
  alpha <- (n / 2) + c
  beta_phi <- (sum(y^2) / 2) - (beta0[i] * sum(y)) - (beta1[i] * sum(y * x)) +
               ((n * beta0[i]^2) / 2) + (beta0[i] * beta1[i] * sum(x)) +
               (((beta1[i]^2) * sum(x^2)) / 2)
  beta_phi1 <- beta_phi + d
  
  # Amostra de \(\phi\)
  phi[i] <- rgamma(1, alpha, beta_phi1)
}
# Beta0
plot(beta0, type="l", xlab="")
abline(h=beta0.real, lwd=3, col=2)

# Beta1
plot(beta1, type="l", xlab="")
abline(h=beta1.real, lwd=3, col=2)

# Phi
plot(phi, type="l", xlab="")
abline(h=phi.real, lwd=3, col=2)

Tratamento da Cadeia e Diagnóstico Amostra para Aquecimento (Burn-in) Definimos o valor do burn-in para eliminar os valores iniciais.

burnin <- 100
beta0_aq <- beta0[(burnin + 1):ite]
beta1_aq <- beta1[(burnin + 1):ite]
phi_aq <- phi[(burnin + 1):ite]

Autocorrelação Examinamos a autocorrelação das cadeias para avaliar a independência entre amostras sucessivas.

Copiar código

acf(beta0_aq, main = "Autocorrelação Beta 0")

acf(beta1_aq, main = "Autocorrelação Beta 1")

acf(phi_aq, main = "Autocorrelação Phi")

par(mfrow = c(2, 2), oma = c(0, 0, 2, 0))

# Gráfico acumulado para beta0
temp <- beta0_aq - mean(beta0_aq)
temp <- cumsum(temp - mean(temp))
plot(temp, type = "l", xlab = "Iteração", ylab = expression(beta[0]))

# Gráfico acumulado para beta1
temp <- beta1_aq - mean(beta1_aq)
temp <- cumsum(temp - mean(temp))
plot(temp, type = "l", xlab = "Iteração", ylab = expression(beta[1]))

# Gráfico acumulado para phi
temp <- phi_aq - mean(phi_aq)
temp <- cumsum(temp - mean(temp))
plot(temp, type = "l", xlab = "Iteração", ylab = expression(phi))

mtext("Gráficos de Somas Acumuladas", outer = TRUE, cex = 1.5, line = -2)

Conclusão Com o Amostrador de Gibbs, conseguimos estimar os valores de 𝛽 0 β 0 ​ , 𝛽 1 β 1 ​ e 𝜙 ϕ utilizando informações da distribuição condicional de cada parâmetro. As análises de autocorrelação e gráficos acumulados indicam a convergência das cadeias para valores próximos aos valores verdadeiros, demonstrando a eficácia da simulação bayesiana para esse modelo.

Algoritmo para o Amostrador de Gibbs

  1. Defina valores iniciais para os parâmetros \(\beta_0\), \(\beta_1\), e \(\phi\).

  2. Obtenha um novo valor \(\theta^{(i)} = (\beta_0^{(i)}, \beta_1^{(i)}, \phi^{(i)})\), através da geração sucessiva das seguintes distribuições:

    • Distribuição condicional completa de \(\beta_0 | \cdot\)

    • Distribuição condicional completa de \(\beta_1 | \cdot\), utilizando os valores gerados no passo anterior.

    • Distribuição condicional completa de \(\phi | \cdot\), utilizando os valores gerados nos passos anteriores.

  3. Repita o processo a partir do passo 2 até que a cadeia de amostras atinja a convergência.