I3_informe.knit

Introducción

En la inferencia Bayesiana, normalmente se conoce \(f(y|\theta)\) y \(p(\theta)\) donde \(y\) son las observaciones y \(\theta\) los parámetros, en muchas ocasiones por la complejidad que puede adquirir la distribución a posteriori \(p(\theta|y)\) es dificil de muestrearla, por lo que se requiere de métodos iterativos para poder obtener la forma de la distribución y poder resumir sus características. En estos casos, la aproximación se puede realizar a través de un Gibbs Sampling, ya que este algoritmo se construye en base a una secuencia que depende de valores de parámetros provenientes de una distribucón que sabemos que converge.

En este informe trataremos de demostrar que el Gibbs Sampling, es efectivo bajo ciertas condiciones, ya que, demostraremos que para ciertos valroes de la varianza el algoritmo no converge o le cuesta más, lo cual se debe a que no se tiene un parámetro mínimo suficiente, lo que no asegura la convergencia del proceso. Por lo que, se analizarán los casos en los que es efectivo y cuales no, a través de diversas simulaciones.

Definción posteriori empírica

Tal como se nos introdujo en clase, se ha estado trabajando en el modelo Bayesiano, el cual esta compuesto por:

Distribución de muestreo:

\[ (Y_i| a,b) \sim N(a + b, \sigma^2), \hspace{0.5 cm} i = 1,...,n; \hspace{1 cm} \perp\perp_{1 \leq i \leq n} Y_i | (a, b)\]

Especificación a priori de los parámetros

\[ a \sim N(\mu_a, \sigma^2_a),\hspace{0.5 cm} b \sim N(\mu_b, \sigma_b^2), \hspace{0.5 cm} a ⊥⊥ b \]

Si queremos obtener la distribución a posteriori empírica de los parámetros, podemos trabajar bajo el supuesto de que como todas las \(Yes\) distribuyen Normal y los parámetros a priori también lo hacen, se puede trabajar bajo el supuesto que a posteriori seguirán distribuyendo Normal, por lo que tenemos:

\[E(a | Y_1) = E(a) + Cov(a,Y_1)Var[Y_1]^{-1}(Y_1-E(Y_1))\] \[E(a| Y_1) = \mu_a + \frac{\sigma_a^2}{\sigma_a^2+\sigma_b^2+\sigma^2} (Y_1 - \mu_a -\mu_b)\] \[E(a| Y_1) = \frac{ (\sigma_a^2 + \sigma^2)\mu_a + \sigma_a^2(Y_1 -\mu_b) }{\sigma_a^2+\sigma_b^2+\sigma^2}\] Luego para la varianza tenemos:

\[Var(a | Y_1) = Var(a) + Cov(a,Y_1)Var[Y_1]^{-1}Cov(Y_1,a)\]

\[Var(a| Y_1) = \sigma_a + \frac{\sigma_a^2\sigma_b^2}{\sigma_a^2+\sigma_b^2+\sigma^2}\] \[Var(a| Y_1) = \frac{\sigma_a^2(\sigma_b^2 + \sigma^2)}{\sigma_a^2+\sigma_b^2+\sigma^2}\] Por lo que podemos concluir que la distribución es:

\[(a|Y_1) \sim N\bigg( \frac{ (\sigma_a^2 + \sigma^2)\mu_a + \sigma_a^2(Y_1 -\mu_b) }{\sigma_a^2+\sigma_b^2+\sigma^2},\frac{\sigma_a^2(\sigma_b^2 + \sigma^2)}{\sigma_a^2+\sigma_b^2+\sigma^2}\bigg)\]

Analogamente para el caso de \((b|Y_1)\) y \((a+b|Y_1)\) tenemos:

\[(b|Y_1) \sim N\bigg( \frac{ (\sigma_b^2 + \sigma^2)\mu_b + \sigma_b^2(Y_1 -\mu_a) }{\sigma_a^2+\sigma_b^2+\sigma^2},\frac{\sigma_b^2(\sigma_a^2 + \sigma^2)}{\sigma_a^2+\sigma_b^2+\sigma^2}\bigg)\] y

\[(a+b|Y_1) \sim N\bigg( \frac{ \sigma^2(\mu_a^2 + \mu_b^2) + Y_1(\sigma_a^2 -\sigma_b^2)}{\sigma_a^2+\sigma_b^2+\sigma^2}, \frac{\sigma^2(\sigma_a^2 + \sigma_b^2)}{\sigma_a^2+\sigma_b^2+\sigma^2}\bigg)\]

Claramente aquí estamos basando nuestro análisis solo para el caso univariado, es decir, solo para \(Y_1\), por lo que si lo generalizamos para \(n\) \(Yes\) tenemos:

\[(a|Y_i) \sim N_i\bigg( \frac{ (\sigma_a^2 + \sigma^2)\mu_a + \sigma_a^2(Y_i -\mu_b) }{\sigma_a^2+\sigma_b^2+\sigma^2},\frac{\sigma_a^2(\sigma_b^2 + \sigma^2)}{\sigma_a^2+\sigma_b^2+\sigma^2}\bigg)\] ,

\[(b|Y_i) \sim N_i\bigg( \frac{ (\sigma_b^2 + \sigma^2)\mu_b + \sigma_b^2(Y_i -\mu_a) }{\sigma_a^2+\sigma_b^2+\sigma^2},\frac{\sigma_b^2(\sigma_a^2 + \sigma^2)}{\sigma_a^2+\sigma_b^2+\sigma^2}\bigg)\] y

\[(a+b|Y_i) \sim N_i\bigg( \frac{ \sigma^2(\mu_a^2 + \mu_b^2) + Y_i(\sigma_a^2 -\sigma_b^2)}{\sigma_a^2+\sigma_b^2+\sigma^2}, \frac{\sigma^2(\sigma_a^2 + \sigma_b^2)}{\sigma_a^2+\sigma_b^2+\sigma^2}\bigg)\]

De modo, que cada distribución será una coordenada del vector y tendremos \(n\) entradas.

Definición de posteriori con Gibbs Sampleing

Ahora que ya sabemos como distribuye empíricamente nuestros parámetros, mediante el proceso de Gibbs Sampling intentaremos explicar si esto es consistente, para ellos revisaremos distintos casos.

Lo primero a hacer es definir la distribución que tendrán los parámetros a posteriori, pero cuando son simulados desde un Gibbs Sampler, serán:

\[ (a^{(k)}|Y_i,b^{(k-1)}) \sim N\bigg(\frac{\mu_a + \sigma_a^2(Y_i - b^{(k-1)})}{ \sigma_a^2 + \sigma^2},\frac{\sigma_a^2}{\sigma_a^2+\sigma^2} \bigg) \]

\[ (b^{(k)}|Y_i,a^{(k-1)}) \sim N\bigg(\frac{\mu_b + \sigma_b^2(Y_i - a^{(k-1)})}{ \sigma_b^2 + \sigma^2},\frac{\sigma_b^2}{\sigma_b^2+\sigma^2} \bigg) \]

\[ (a^{(k)}+b^{(k)}|Y_i,a^{(k-1)},b^{(k-1)}) = (a^{(k)}|Y_i,b^{(k-1)})+(b^{(k)}|Y_i,a^{(k-1)}) \]

Con esto podemos analizar la convergencia del algoritmo, para cuando las medias son 0 y las desviaciones estándar 1:

# Definimos el tamaño del vector Y

n = 10000

# Definimos los hiperparámetros mu_a, mu_b, sigma_a, sigma_b y sigma

mu_a = 0
mu_b = 0 
sigma_a = 1
sigma_b = 1
sigma = 1 

# Definimos los valores iniciales de a y b
# 
a_0 = rnorm(1,mu_a,sigma_a)
b_0 = rnorm(1,mu_b,sigma_b)

# Definimos el vector de las observaciones 

Y_obs <- rnorm(n,a_0 + b_0, sigma)

# Con ello empezamos el algoritmo de Gibbs para obtener la dist a posteriori
# de la (a|Y_obs),(b|Y_obs) y (a+b|Y_obs)

a_sim = vector(mode = "numeric", length = n)
b_sim = vector(mode = "numeric", length = n)
a_b_sim = vector(mode = "numeric", length = n)

a_sim[1] = a_0 
b_sim[1] = b_0
a_b_sim[1] = a_0 + b_0

i = 2
k = 21000

while (i<=k) {
  a_sim[i] <- rnorm(1,mean = (mu_a + sigma_a^2*(Y_obs - b_sim[i-1]))/(sigma_a^2+sigma^2),
                    sd = sqrt((sigma_a^2)/(sigma_a^2+sigma^2)))
  b_sim[i] <- rnorm(1,mean = (mu_b + sigma_b^2*(Y_obs - a_sim[i-1]))/(sigma_b^2+sigma^2),
                    sd = sqrt((sigma_b^2)/(sigma_b^2+sigma^2)))
  a_b_sim[i] <- a_sim[i] + b_sim[i]
  i=i+1  
}

a_sim_20000 <- a_sim[1001:21000]
b_sim_20000 <- b_sim[1001:21000]
a_b_sim_20000 <- a_b_sim[1001:21000]

pos <- seq(1,20000, by = 10)

a_sim_2000 <- c()
b_sim_2000 <- c()
a_b_sim_2000 <- c()

for (i in pos) {
  aux1 <- a_sim_20000[i]
  aux2 <- b_sim_20000[i]
  aux3 <- a_b_sim_20000[i]
  a_sim_2000 <- c(a_sim_2000,aux1)
  b_sim_2000 <- c(b_sim_2000,aux2)
  a_b_sim_2000 <- c(a_b_sim_2000,aux3)
}

Una vez obtenidos los valores de las parámetros a través de la simulación de Gibbs, podemos analizar su convergencia mediante un traceplot():

medias = 0 y sigmas = 1

Aquí podemos ver que el parámetro de \(a\) converge a 0.3259672 y \(b\) converge a 0.3360609, mientras que el mínimo suficiente \(a+b\), converge a 0.6620281 que coincide exactamente con la suma de los parámetros \(a\) y \(b\) (0.6620281) , por lo que hasta el minuto podemos ver que todo va bien encaminado. Ahora que tenemos las distribuciones a posteriori mediante un qqplot() y el ks.test() las podemos comparar en base a la distribución empírica original:

Para el parámetro \(a\), tenemos que el qqplot() se acerca mucho a la Normal, como se puede ver en el siguiente gráfico:

y analizando el resultado que nos da el ks.test() es claro que el parámetro si distribuye Normal, ya que su valor es demasiado pequeño, como se muestra a continuación:


    Asymptotic one-sample Kolmogorov-Smirnov test

data:  datos_ordenados_a
D = 0.60929, p-value < 2.2e-16
alternative hypothesis: two-sided

Para el parámetro \(b\), tenemos que el qqplot() se acerca mucho a la Normal:

Es claro observar que practicamente todos los valores son cercanos a la Normal, nuevamente si analizamos de manera cuantitativa usando el ks.test(), el valor es demasiado pequeño lo que hace validar nuestro supuesto.


    Asymptotic one-sample Kolmogorov-Smirnov test

data:  datos_ordenados_b
D = 0.59964, p-value < 2.2e-16
alternative hypothesis: two-sided

Para el parámetro \(a+b\), tenemos que el qqplot() se acerca mucho a la Normal, mucho más que los parámetros de forma individual, esto tambien se debe principalmente a que es el mínimo suficiente:

Nuevamente en este caso la gráfica es practicamente perfecta, y el ks.test() nos confirma esta idea de que el parámetro distribuye Normal, dado su tamaño extremadamente pequeño:


    Asymptotic one-sample Kolmogorov-Smirnov test

data:  datos_ordenados_a_b
D = 0.29368, p-value < 2.2e-16
alternative hypothesis: two-sided

Finalmente, si nos interesa revisar si las distribuciones a posteriori de \((a|Y_1,...,Y_n)\) y de \((b|Y_1,...,Y_n)\) son función o se parecen en algo a la distribución a posteriori del parámetro mínimo suficiente podemos usar un histograma para revisar su comportamiento:

Con lo anterior, podemos ver que el algoritmo de Gibbs aplicado a una cadena de 10^{4} observaciones, cuyas medias son \(0_s\) y sus desviaciones estándar son \(1_s\), todo resulta perfecto. Sin embargo, ¿Qué ocurre cuando la varianza comienza a aumentar?. A continuación, ampliaremos el análisis para distintos casos donde la desviación estándar comienza a cambiar (para revisar todos los casos ver Anexos):