Saturday, October 20, 2012

La distribución predictiva en WinBugs

A la hora de modelar un conjunto de datos, uno de los intereses principales del investigador, no sólo está en encontrar la esperanza (o la esperanza condicional, si hay covariables y otros parámetros) de las variables aleatorias, sino también encontrar la distribución condicional de nuevos datos en el modelo para realizar un proceso de predicción.


El ejemplo más claro que se me ocurre es en inferencia de poblaciones finitas. Por ejemplo, nótese que un total poblacional (y tomando como medida de probabilidad la distribución de los datos ignorando el mecanismo de selección de la muestra) se puede escribir de la siguiente forma:


$latex T=sum_{k in s}Y_k+sum_{k notin s}Y_k$


Por supuesto, acá se supone que T es una variable aleatoria y por lo tanto, una realización de T será:


$latex t=sum_{k in s}y_k+sum_{k notin s}y_k$


Entonces, después de la recolección de la muestra, el primer sumando será una cantidad conocida pero la segunda parte seguirá desconocida y deberá predecirse dependiendo del modelo que se haya considerado. Es aquí donde entra la distribución predictiva. Desde el punto de vista Bayesiano, cuando esta se define, es posible simular los valores de la segunda parte de la suma y con esto obtener una predicción del total poblacional. Por el lado frecuentista, es posible obtener una estadística pivote para encontrar la probabilidad de que una nueva observación se encuentre dentro de un intervalo específico. En ambos escenarios lo que uno esperaría es que los intervalos de predicción fuesen más anchos que los intervalos de credibilidad o de confianza.


La función de distribución predictiva está definida por la siguiente integral (no trivial en la mayoría de los casos):


$latex p(y|bold{Y})=int_{theta} p(y|theta)p(theta|bold{Y})$


Sin embargo, esta integral se puede ver como una esperanza condicional de la función de $latex p(y|theta)$. Es decir,


$latex p(y|bold{Y})=E_{theta|bold{Y}}( p(y|theta))approx frac{sum_i p(y|theta_i)}{N}$


Lo anterior permite aproximar esta distribución usando métodos MCMC y evadir el cálculo analítico, que muchas veces es intratable. Este procedimiento se enuncia a continuación:




  1. Generar un valor Θ desde la distribución posterior usando MCMC.

  2. Utilizar ese valor específico y generar un nuevo valor de y desde la verosimilitud p(y|Θ).

  3. Repetir 1 y 2 miles de veces.

  4. Aproximar la función predictiva con todos los valores generados de las verosimilitudes.

  5. Realizar la predicción.


En WinBugs es muy sencillo; aunque parecer ser que no hay mucho escrito en la red de cómo aproximar la función predictiva. Un ejemplo sencillo, es asumir que se tiene una regresión lineal normal. Entonces, el siguiente código hace todo el proceso automáticamente:



model{
for(i in 1:N){
y[i]~dnorm(mu[i], tau)
mu[i]<-b1+b2*x[i]
y.pred[i]~dnorm(mu[i], tau)
             }

b1 ~ dnorm(0, 0.0001)
b2 ~ dnorm(0, 0.0001)
tau ~ dgamma(0.01, 0.01)
}

DATA
list(y=c(51.25, 56.79, 59.32, 63.09, 68.82, 67.84, 69.14, 74.18, 74.76, 76.71, 75.18, 76.23, 79.15),x=c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13), N=13)

Nótese la inclusión de y.pred en el código. Este objeto también debe incluirse en el cuadro de dialogo samples y para obtener las predicciones, basta con hacer clic en el botón stats del cuadro de dialogo samples. Ahora, para obtener los intervalos de predicción, es necesario ir al menú inference y elegir la opción compare. En el cuadro de dialogo emergente escribimos y.pred en node, y en other (para que se muestren los puntos), x en axis. Finalmente clic en model fit. 



Si se quiere, también es posible obtener los intervalos de credibilidad al escribir mu en node. A continuación se muestra la estimación de la media (junto con los intervalos de credibilidad) y la estimación de la predicción (junto con los intervalos de predicción).

No comments:

Post a Comment