Paralelización longitudinal de modelos de regresión bayesiana con muestras grandes

[Víctor de Buen Remiro]

Cada vez es más habitual encontrarse con modelos que sin tener una estructura jerárquica compleja con muchos nodos, ni tampoco un excesivo número de variables sí presentan problemas de tamaño debido a que tienen una longitud muestral de cientos de miles o millones de datos.

Está claro que para problemas con tantos datos hay que montar un sistema de paralelización que permita estimar en tiempos razonables y que evite problemas debidos a la falta de RAM.

[Jorge]

Se ha adjuntado al ticket un paper​ que describe un esquema de estimación muestreando las ecuaciones.

[Víctor de Buen Remiro]

En modelos con residuos independientes como regresión lineal, probit, logit, etc.; en los que no importa el orden de los datos, podría pensarse en una estrategia de tipo divide&conquer muy sencilla:

1. Divide : En lugar de intentar simular el modelo en bloque, se dividirá en S submodelos del siguiente modo:
   1. Se divide aleatoriamente el bloque lineal de la regresión de longitud M en S segmentos de longitud m = [M/S] o bien m+1, pues si la división no es exacta se reparte el resto. Tendremos pues matrices output e input con unos m datos e ncada segmento. De hecho es preferible que la matriz conjunta de input no se cree en ningún momento pues en tal caso podríamos tener problemas de RAM. Por ejemplo, si los datos provienen de una base de datos, deberían extraerse directamente los segmentos por separado.
   2. El resto de bloques del modelo, priors, restricciones, etc.; se replican en cada submodelo y se almacena en disco o base de datos de forma que sea accesible por todas las máquinas que vayan a participar.
   3. Se lanza cada submodelo por separado en la máquina que se le asigne generando cada una de ellas una cadena MCMC que se almacenará localmente para mayor agilidad.
2. Conquer: Ahora hay que mezclar las cadenas en una sola, de forma que la distribución a posteriori sea la correspondiente al modelo conjunto. Y aquí es donde aparece el problema estadístico de verdad. En el [​documento que aporta Jorge se dan distintas aproximaciones para el caso jerárquico lineal puro pero habría que ver hasta qué punto nos sirven en nuestros modelos llenos de restricciones y cosas no lineales.

En modelos con residuos correlados como el caso ARIMA la división tendría que hacerse en segmentos de datos consecutivos y no aleatorios. Además la mezcla sería más complicada porque tendría que tener en cuenta la matriz de covarianzas, pero también sería posible.

Otro problema que surge es que algunas variables input podrían hacerse cero en un subsegmento, o bien podrían darse otro tipo de colinealidades más complejas también coyunturales a un subsegmento dado. Para evitarlo habría que introducir cierta información a priori que resolviera el problema numérico sin aportar demasiada arbitrariedad a no ser que realmente se disponga de información creíble.

[Víctor de Buen Remiro]

En el caso de regresión lineal normal con residuos independientes está claro que, tal como dice ​Gelman en el punto 2.1 de la página 2, las medias ponderadas de las estimaciones de los parámetros en cada submuestra son las medias para la muestra total, tanto para los parámetros del bloque lineal como para la varianza.

Esto quiere decir que si concatenamos las cadenas obtenidas para cada submuestra obtenemos una cadena de la muestra total con el mismo contenido informativo en lo que se refiere a los valores centrales.

Pero nosotros no queremos estimar sólo la media, lo que queremos es construir una MCMC que sea representativa de la distribución conjunta a posteriori. Si el resultado para las medias pudiera extenderse a cualquier cuantil podríamos concluir que la cadena mezcla de las cadenas generadas para las submuestras es una cadena válida para la muestra global.

Ahora mismo no se me ocurre cómo demostrar teóricamente que es así, pero no veo ningún motivo para que no se cumpla y siempre podemos estudiarlo empíricamente en casos reales y simulados.

El único pero a este método de mezclado a toro pasado, es la longitud muestral, pues obviamente si tenemos L simulaciones de cada submuestra obtenemos un total de S*L simulaciones de la muestra global que en realidad darían la misma precisión que una cadena de longitud L. Esto quiere decir, que los criterios de convergencia a tener en cuenta serían los de la submuestra que más simulaciones requiera, no los de la cadena mezclada.

Como estaríamos hablando de miles de simulaciones creo que es aplicable el teorema central del límite que extendería este método de mezclado de cadenas a cualquier regresión lineal generalizada con residuos indendientes (probit, logit, poisson, ...)
También es algo que habría que tratar de comprobar al menos empíricamente.

[Víctor de Buen Remiro]

He comenzado un documento para el tratamiento teórico de este problema:
​Combinación de MCMC independientes