Métodos Métodos de simulación basados Lección 2: Variables Aleatorias en Cadenas de Markov Prof. Prof. Daisy Da isy Arroyo Ar royo Universidad de Concepción Departamento de Estadística Facultad Facultad de Ciencias Físicas y Matemáticas Junio 2017
Contenido
•
•
•
Cadenas de Markov Algoritmo de Metrópolis-Hast Metrópolis-Hastings ings Muestreo de Gibbs
Ejemplo Algoritmo Metrópolis
Generar la distribución -Student con libertad, usando la distribución .
grados de propuesta
Ejemplo Veremos el efecto de las diferentes elecciones de varianza de la distribución propuesta, repitiendo la simulación con diferentes elecciones de . La densidad Entonces:
es proporcional a
Código Parámetros: , Valor inicial: Longitud de la cadena: Metropolis <- function(n, sigma, x0, N) { x <- numeric(N) x[1] <- x0 u <- runif(N) k <- 0 for (i in 2:N) { y <- rnorm(1, x[i-1], sigma) if (u[i] <= (dt(y, n) / dt(x[i-1], n))) x[i] <- y else { x[i] <- x[i-1] k <- k + 1 } } return(list(x=x, k=k)) }
Convergencia La convergencia del algoritmo de Metrópolis a menudo es sensible a la elección del parámetro de escala. •
•
Cuando la varianza de los incrementos es muy grande, muchos de los puntos candidatos son rechazados y el algoritmo es muy ineficiente. Si la varianza del incremento es muy pequeño, los puntos candidatos son casi todos aceptados, el algoritmo genera una cadena que es casi como un verdadero, lo cual también es ineficiente.
Un enfoque para seleccionar el parámetro de escala es monitorear las tasas de aceptación, los cuales deberían estar en el rango
Convergencia La varianza de la distribución
es
Si consideramos , entonces la desviación estándar de la distribución objetivo es .
Ejemplo
Comparando cuantiles En general, los problemas MCMC no tienen los cuantiles teóricos de la distribución objetivo disponible para comparar. Vamos a comparar los resultados del ejemplo con los cuantiles teóricos de la distribución objetivo. Descartamos las primeras 500 filas de cada cadena.
Comparando cuantiles Los cuantiles son calculados por medio de la función apply (aplicando quantile a las columnas de la matriz). a <- c(0.05, seq(0.1,0.9,0.1), 0.95) Q <- qt(a,n) w <- c(w1$x, w2$x, w3$x, w4$x) mc <- w[501:N, ] Qw <- apply(mc, 2, function(x) quantile(x,a)) print(round(cbind(Q, Qw),3))
Contenido
•
•
•
Cadenas de Markov Algoritmo de Metrópolis-Hastings Muestreo de Gibbs
Introducción El muestreador de Gibbs fue descrito por los hermanos Geman & Geman (1984) en la aplicación del análisis de Gibbs en distribuciones lattice. El muestreador Gibbs se aplica a menudo cuando el objetivo es una distribución multivariadas . Supongamos que todas las densidades condicionales univariadas están completamente especificadas y es razonablemente fácil tomar muestras de ellas. La cadena se genera por muestreo a partir de las distribuciones marginales de la distribución objetivo y, por lo tanto, se acepta cada punto candidato.
Introducción Sea •
Definamos vectores aleatorios ()
•
un vector aleatorio.
Denotamos la correspondiente densidad condicional univariada de dado ()
•
dimensionales:
()
El muestreador de Gibbs genera la cadena por muestreo desde cada una de las densidades condicionales () .
Algoritmo 1. Inicializamos por
en el tiempo
2. Para cada iteración indexada (a) Fijar
repetir
.
(b) Para cada coordenada i.
Generar
ii.
Actualizar
∗
(c) Fijar (d) Incrementar
desde
()
∗ (cada candidato es aceptado)
Ejemplo Generar una distribución normal bivariada con vector media varianzas , y correlación usando el muestreador de Gibbs.
Ejemplo En el caso bivariado:
Las densidades condicionales de una distribución normal bivariada son univariadas normales con parámetros
Ejemplo La cadena es generada por muestreo desde
Para una distribución bivariada Gibbs: 1. 2. 3. 4. 5.
Fija Genera ∗ Actualiza Genera ∗ Fija
desde
∗
desde
∗
∗
,
en cada iteración el muestreador de
Ejemplo
Cadena normal bivariada generada con el Muestreador de Gibbs
Muestreador de Gibbs en Geoestadística Principio
Se busca simular un vector aleatorio con componentes. La transición de un estado al estado siguiente se define como: 1) Seleccionar una componente al azar en 2) Plantear , salvo para la -ésima componente, la cual se simula condicionalmente a las otras componentes. Esto asegura que la distribución sea invariante por el núcleo de transición. El núcleo de transición también es aperiódico (existe una probabilidad no nula de seleccionar varias veces consecutivas la misma componente , luego no puede haber periodicidad). Falta entonces comprobar la propiedad de irreducibilidad para asegurar la convergencia.
Muestreador de Gibbs en Geoestadística Si es un vector Gaussiano de media y matriz de varianza-covarianza , la transición del muestreador de Gibbs se define como: 1) Seleccionar una componente al azar en 2) Realizar el kriging simple de la -ésima componente de a partir de las otras componentes. Se obtiene un valor estimado y una varianza de error. 3) Plantear salvo para la n-ésima componente, la cual se simula como una variable Gaussiana de media igual al valor estimado por kriging y varianza igual a la varianza de error de kriging.
Muestreador de Gibbs en Geoestadística Dado que la densidad de probabilidad Gaussiana es positiva sobre , se puede pasar de cualquier estado a cualquier otro en un número finito de transiciones. Por lo tanto, el núcleo de transición es irreducible y tiene la distribución multi-Gaussiana deseada como su límite ergódico.
Ejemplo
Ejemplo Grilla
Ejemplo El muestreador de Gibbs se utiliza para generar valores Gaussianos correspondientes a litotipos/facies en los puntos de muestreo (pozos, perforaciones).
Datos
Simulado
Real
Referencias •
•
Geof H. Givens and Jennifer Hoeting (2005). Computational Statistics . Wiley series in Probability and Statistics. Maria L. Rizzo (2008). Statistical Computing with R . Chapman & Hall / CRC ISBN 9781584885450
