METRÓPOLIS HASTING

CASO

Se obtuvieron los siguientes valores considerando una muestra aleatoria de tamaño 7: 1.327, 4.266, 7.148, 6.422, 3.301, 7.816 y 5.076. Además, la muestra se obtuvo desde la distribución T de Student con 1 grado de libertad \(t_1(\theta, 1)\). Hallar la distribución posterior si \(θ\) tiene como distribución a priori una distribución uniforme.

Resolución:

Distribución Fundamental: Distribución \(t_1(\theta, 1)\)

\[ f(y|\theta, \sigma) = \frac{\Gamma((v+1)/2)}{\Gamma(v/2)\sqrt{v\pi }\sigma }\left (1+\frac{1}{v}\left ( \frac{y-\theta}{\sigma} \right )^{2}\right )^{-(v+1)/2}\]

\[ f(y|\theta) = \frac{\Gamma(1)}{\Gamma(1/2)\sqrt{1\pi } (1) }\left (1+\frac{1}{1}\left ( \frac{y-\theta}{1} \right )^{2}\right )^{-1}\]

\[ f(y|\theta) = \frac{0!}{\sqrt{\pi }\sqrt{\pi } }\left (1+ ({y-\theta})^{2}\right )^{-1}\]

\[ f(y|\theta) = \frac{1}{\pi\left (1+ ({y-\theta})^{2}\right )}\]

\[ \theta \sim CAUCHY(\theta, 1)\]

Distribución A priori: Distribución Uniforme

\[f(\theta) \quad \overline{\alpha }\qquad k\]

Función de Verosimilitud:

\[L(\theta, \underline{y}) = \prod_{i=1}^{n}\frac{1}{\pi\left (1+ ({y_i-\theta})^{2}\right )}\]

\[L(\theta, \underline{y}) = \prod_{i=1}^{7}\frac{1}{\pi\left (1+ ({y_i-\theta})^{2}\right )}\]

\[L(\theta, \underline{y}) =\frac{1}{\pi^7}\times \frac{1}{\prod_{i=1}^{7} \left (1+ ({y_i-\theta})^{2}\right )}\] Distribución Posterior:

\[f(\theta, \underline{y}) \quad \overline{\alpha }\qquad L(\theta, \underline{y})f(\theta)\]

\[f(\theta, \underline{y}) \quad \overline{\alpha }\qquad \frac{1}{\prod_{i=1}^{7} \left (1+ ({y_i-\theta})^{2}\right )}\]

\[f(\theta, \underline{y}) = k \times \frac{1}{\prod_{i=1}^{7} \left (1+ ({y_i-\theta})^{2}\right )} \qquad -\infty\leq \theta\leq \infty \]

# Función objetivo
g <- function(theta){
  1/((1 + (1.327 - theta)^2)*
       (1 + (4.266 - theta)^2)*
       (1 + (7.148 - theta)^2)*
       (1 + (6.422 - theta)^2)*
       (1 + (3.301 - theta)^2)*
       (1 + (7.816 - theta)^2)*
       (1 + (5.076 - theta)^2))
}
k <- 1/integrate(g, lower = -Inf, upper = Inf)$value
k
## [1] 5152.568

\[f(\theta, \underline{y}) = 5152.568 \times \frac{1}{\prod_{i=1}^{7} \left (1+ ({y_i-\theta})^{2}\right )} \qquad -\infty\leq \theta\leq \infty \]

Aplicando regla de normalización de las distribuciones

\[\int_{-\infty }^{\infty} \quad 5152.568 \times \frac{1}{\prod_{i=1}^{7} \left (1+ ({y_i-\theta})^{2}\right )} = 1\]

f.post <- function(theta){
  k/((1 + (1.327 - theta)^2)*
       (1 + (4.266 - theta)^2)*
       (1 + (7.148 - theta)^2)*
       (1 + (6.422 - theta)^2)*
       (1 + (3.301 - theta)^2)*
       (1 + (7.816 - theta)^2)*
       (1 + (5.076 - theta)^2))
}
round(integrate(f.post, lower = -Inf, upper = Inf)$value,3)
## [1] 1
## Gráfica de la distribución posterior
curve(f.post, from = 0, to = 10, xlab = "theta",ylim = c(0, 0.5),
      ylab = "Distribución posterior", col = "blue")

Asimismo, aplicar el algoritmo de Metropolis-Hastings para generar una muestra de la distribución posterior usando la siguiente semilla \(set.seed(20241)\), con 4500 simulaciones de la distribución propuesta \(Logistic(\theta,\sigma^2=1)\) y el estado inicial \(\theta^{(0)}\) = 2.

Resolución:

## Algoritmo Metrópolis Hasting
MetropolisHasting <- function(m, theta.init, sd){
  ## Paso 1. Inicializar
  theta.acept <- numeric(m)
  acept <- 0
  theta.actual <- theta.init
  ## Paso 2. Repetir
  for (i in 1:m){
    ## Paso 2.1. Simular un valor candidato a partir de la distribución propuesta (Dist. Logistica)
    theta.cand <- rlogis(1, location   = theta.actual, scale   = sd)
    ## Paso 2.2. Calcular la probabilidad de aceptar el valor simulado
    g.num <- g(theta.cand)
    g.den <- g(theta.actual)
    q.num <- dlogis(theta.actual, location  = theta.cand, scale   = sd)
    q.den <- dlogis(theta.cand, location  = theta.actual, scale   = sd)
    p <- (g.num*q.num)/(g.den*q.den)
    alpha <- min(1, p)
    ## Paso 2.3. Obtener u
    u <- runif(1, min = 0, max = 1)
    ## Paso 2.4.
    if (u < alpha){
      theta.actual <- theta.cand         ## Se acepta el candidato
      acept <- acept + 1                 ## Se actualiza el número de candidatos aceptados
      theta.acept[acept] <- theta.actual ## Se guarda el valor generado
    }
  }
  list(theta = theta.acept[1:acept], aceptados = acept/m)
}
RNGkind(sample.kind = "Rejection")
set.seed(20241)
f.post.MH <- MetropolisHasting(m = 4500, theta.init = 2, sd = 1)

Proporción de valores aceptados

f.post.MH$aceptados
## [1] 0.5484444

Histograma para los valores simulados

hist(f.post.MH$theta, freq = FALSE, xlab = "theta", ylim = c(0, 0.5), main = " ",xlim = c(0, 10))

plot(density(f.post.MH$theta), ylim = c(0, 0.5), col = "red", main = " ")
curve(f.post, from = 0, to = 10, col = "blue", add = TRUE)
legend("topright", legend = c("Dist. Posterior MH", "Dist. Posterior"), 
       col = c("red", "blue"), lty = 1)

Comportamiento de los valores simulados

library(coda)
traceplot(as.mcmc(f.post.MH$theta), ylim = c(0, 11))

Repetir el algoritmo considerando ahora como estado inicial (\(\theta^{(0)}\)) el nivel en donde los valores simulados alcanzan la estabilidad, 9000 simulaciones y \(set.seed(20242)\).

Proporción de valores aceptados

## Cambiando el valor inicial theta.init = 5
RNGkind(sample.kind = "Rejection")
set.seed(20242)
f.post.MH <- MetropolisHasting(m = 9000, theta.init = 5, sd = 1)
f.post.MH$aceptados
## [1] 0.5602222

Histograma para los valores simulados

hist(f.post.MH$theta, freq = FALSE, xlab = "theta", ylim = c(0, 0.5), main = " ",xlim = c(0, 10))

plot(density(f.post.MH$theta), ylim = c(0, 0.5), col = "red", main = " ")
curve(f.post, from = 0, to = 10, col = "blue", add = TRUE)
legend("topright", legend = c("Dist. Posterior MH", "Dist. Posterior"), 
       col = c("red", "blue"), lty = 1)

Comportamiento de los valores simulados

library(coda)
traceplot(as.mcmc(f.post.MH$theta), ylim = c(0, 11))

Se utilizó el algoritmo de Metropolis-Hastings para generar una muestra de la distribución posterior, usando una distribución propuesta logística con parámetro de locación θ y parámetro de escala igual a 1. Al ejecutar el algoritmo, se obtuvo una proporción de valores aceptados igual a 0.5484444, es decir, aproximadamente el 55% de los valores candidatos fueron aceptados. Además, la gráfica de densidad de los valores aceptados se ajustó bien a la distribución posterior. También se logró graficar el comportamiento de los valores simulados a lo largo de las iteraciones, y se observó una cierta tendencia a estabilizarse alrededor de un valor medio de 5.

Posteriormente, se repitió el algoritmo considerando ahora como estado inicial el nivel en donde los valores simulados alcanzan la estabilidad. Para ello, cambiamos el valor inicial a 5, que es el valor aproximado donde se estabilizan los valores simulados. Al ejecutar la función con los nuevos parámetros, se obtuvo una proporción de valores aceptados igual a 0.5602222, lo que significa que aproximadamente el 56% de los valores candidatos fueron aceptados. Además, al graficar la densidad de los valores aceptados y el comportamiento de los valores simulados, como antes. Se pudo apreciar resultados muy similares a los obtenidos con el valor inicial anterior. Sin embargo, para este caso los valores aceptados se ajustaron bastante bien a la distribución posterior, lo que indica que el algoritmo ha convergido a la distribución estacionaria.

  1. Universidad Nacional Agraria La Molina, ↩︎

  2. Universidad Nacional Agraria La Molina, ↩︎

  3. Universidad Nacional Agraria La Molina, ↩︎

  4. Universidad Nacional Agraria La Molina, ↩︎

  5. Universidad Nacional Agraria La Molina, ↩︎