Kriging con Stan

Este mes de julio, cuórum mediante, impartiré en la UPC un curso que he maltitulado, mor de brevedad, Estadística Bayesiana Aplicada.

Los cursos de estadística bayesiana son teoría, mucha teoría, y unos ejemplos tontos que quieren justificarla. Del tipo: hagamos lo que ya sabemos hacer de otra manera más; busquemos una alternativa molona al p-valor (y usémosla como usaríamos un p-valor, por supuesto), etc.

Mi curso debería haberse titulado algo así como: Problemas reales (aunque simplificados por motivos estrictamente pedagógicos) resueltos con tecnología bayesiana porque, si no, dígame Vd. cómo lo haría: ¿con optim? Jajajajaja… Es un curso en el que aplicaciones reales justifican la tecnología y esta, a su vez, la teoría. Y sí, enlazará con teoría de la decisión propiamente formulada, remitirá (a lo más) a libracos, y enseñará a resolver los problemas con Stan.

Uno de los problemas avanzados pero a la mano que quiero plantear es el del kriging. Por simplificar, en una dimensión. Hay muchas cosas de contar eso del kriging, muchas de las cuales no he entendido jamás, pero la más sencilla es generalizando los modelos lineales generalizados. Por simplificar, desgeneralicemos los GLMs y partamos del modelo lineal,

Y_i | X_i \,\, \sim \,\, N(\mu_i, \sigma)

o bien,

y_i = \mu(x_i) + \epsilon

donde \epsilon se ajusta a la especificaión previa. Si conocemos \mu, podremos hacer predicciones para nuevos valores x.

Aquí las observaciones son independientes. Pero en nuestro problema observamos medidas en determinados puntos (en un terreno (2D), en una línea (1D), etc.) y sabemos que existe cierta regularidad: los puntos próximos están correlacionados entre sí.

El objetivo del kriging es poder realizar estimaciones en puntos no observados utilizando no solo nuestro conocimiento de \mu, sino también aquello que nos aporten los puntos x conocidos próximos al que nos interesa.

Por concretar, dados los valores

podríamos querer obtener estimaciones de la función subyacente en alguno intermedio. O en una rejilla determinada. Sí, como si usásemos loess o similar. Pero, en nuestro caso, usando el siguiente modelo:

Y | X \,\, \sim \,\, N(\mu, \Sigma)

donde \mu es un vector de medias y \Sigma es una matriz de covarianzas en la que \Sigma_{ij} es una función decreciente de la distancia entre x_i y x_j.

Al modelo básico se le pueden añadir cascabeles variados:

  • Sumarle a \Sigma una matriz diagonal que recoja el error de medición.
  • Variables independientes y sus coeficientes en la determinación de la media \mu.
  • Funciones de enlace para que $Y$ sea no el valor deseado sino \log(\lambda) en un modelo de Poisson o la invlogis de una probabilidad.

En nuestro caso, se puede estimar la función subyacente en una rejilla,

con el código

que tira de

En el código anterior \Sigma es una matriz que tiene valores \sigma^2 + r^2 en la diagonal (r es la desviación estándar del ruido de observación) y \sigma^2 \exp(-\phi d_{ij}^{1.5}) fuera de ella.

Así que, en resumen, el kriging es conceptualmente sencillo y computacionalmente asequible. Y si te interesan estas y más cosas y a primeros de julio estás por Barcelona con tu portátil, planteáte aprender y practicar estas cosas por las mañanicas de una semana.

3 comentarios sobre “Kriging con Stan

  1. Jose Luis Cañadas Reche 3 marzo, 2018 18:25

    Ya me gustaría. Que tengo oxidado lo que aprendí de modelos espacio temporales.

  2. Guillermo Santamaria 3 marzo, 2018 19:55

    Hola Carlos, me gustaría mucho asistir pero lamentablemente tengo que estar por Madrid, tienes pensado impartirlo por aquí, aunque sea parcialmente? Muchas gracias!

  3. Carlos J. Gil Bellosta 3 marzo, 2018 22:09

    Pues si sale medianamente bien, lo repito en Madrid… Si hay cuórum, claro.

Los comentarios están desabilitados.