Package BysInfDiag

Métodos de Montecarlo para la inferencia y la diagnosis sobre cadenas de Markov (MCMC).

Convergencia de cadenas

Los métodos MCMC garantizan la convergencia bajo ciertos supuestos pero no dicen nada de en qué momento se alcanzará la misma. Por este motivo es necesario tener algún método diagnóstico que calcule la longitud inicial de muestras anteriores a la convergencia para eliminarlas de la muestra (burn-in).

Existen dos clases fundamentales de métodos para

• sequencia única: Se basan únicamente en una cadena por lo que sólo pueden determinar en cierta medida que la cadena converge a algo, pero no es posible asegurar que converge a la distribución a posteriori, especialmente con distribuciones multimodales o variables vectoriales muy correlacionadas. Los más utilizados son los presentes en el paquete ​CODA del lenguaje ​R
  • Raftery and Lewis
  • Geweke
  • Heidelberg and Welch
• sequencia múltiple: Estos métodos utilizan varias cadenas generadas desde puntos dispersos de la región factible para comprobar que realmente todas las cadenas no sólo convergen cada una por separado sino que además lo hacen a una misma distribución. El método más empleado, por el mismo motivo de aparecer en el CODA, es el
  • Gelman and Rubin

El problema de éste último método es que sólo sirve para distribuciones normales, o bien para aquellas que pueden ser transformadas en normales de una forma sencilla, de tipo Box-Cox por ejemplo. Sin embargo, en los problemas de la vida real esta es una hipótesis demasiado restrictiva, pues aparte de que se pueden manejar otras distribuciones como la chi-cuadrado, logísticas, Poisson, etc.; es muy usual que las normales se trunquen por motivos de coherencia. En cualquier caso, aunque fuera posible esa transformación de la distribución a posteriori, puede que no sea en absoluto trivial hallarla, y menos de forma general sin intervención directa del analista de datos.

Para ser completamente general se precisaría un método independiente de la distribución, como podrían ser métodos de cociente de verosimilitudes u otros de tipo no paramétrico:

​Barton	Comparisson of two samples by a nonparametric likelihood-ratio test
​Wiel	Exact Distributions of Distribution-free Test Statistics
​Ozturk and Balakrishnan	Exact two-sample nonparametric test for quantile difference between two populations based on ranked set samples
​Kuiper	Kuiper two sample test
​Cramér-von Mises	A C++ Program for the Cramér-von Mises Two-Sample Test
​Dette-Wagener-Volgushev	Nonparametric comparison of quantile curves: a stochastic process approach
​parzen	Statistical methods mining and nonparametric quantile domain data analysis
​Kosorok	Two sample quantile test under general conditions
​Mann-Whitney U-test	Test no paramétrico para la comparación de medianas de muestras de tamaños arbitrarios, sustitutivo del test de la t-student para muestras normales

La principal conclusión tras la lectura de estos documentos es que la comparación no paramétrica de distribuciones es un problema tremendamente complejo desde el punto de vista estadísticoy todavía más en lo que se refiere a su implementación. De hecho sólo he encontrado un método libre implementado en C++, el ​Cramér-von Mises, que es realmente bastante robusto pero tremendamente lento para cadenas de apenas unos cientos de muestras. El ​Mann-Whitney U-test) sí que es rapidísimo pero en principio sólo nos sirve para comparar las medianas. Tiene la ventaja añadida de que ya existe en TOL como la función AlgLib.MannWhitneyUtest. Podría usarse como filtro inicial a cualquier otro test de comparación de muestras, pues si las medianas son significativamente distintas es evidente que las distribuciones no pueden ser las mismas. Pero además de eso puede usar de una forma sencilla, aunque un tanto heurística, para comprobar que las muestras se parecen no sólo en los valores centrales sino en todo el rango de valores posibles

1. Se crea el vector de muestra conjunta mediante la simple concatenación de las muestras a chequear.
2. Se divide el rango de valores de la muestra conjunta en un número dado no, más bien pequeño, de partes equi-frecuentes.
3. En cada intervalo se calcula el test Mann-Whitney U-test sobre las muestras actuales censuradas en dicho intervalo.
4. Si en cualquiera de los intervalos se supera el nivel de significación se rechaza la hipótesis nula.

Habría que estudiar cómo se comporta en casos reales y simulados a ver si es posible parametrizarlo para que resulte eficaz en el ámbito de nuestros objetivos.

Mecanismo de generación de puntos sobre-dispersos

Lo primero que hace falta para poder simular varias cadenas es generar unos cuantos puntos de origen de la cadena lo más dispersos que sea posible dentro de la región factible. El problema es que para poder generar unos pocos puntos dispersos sin favorecer demasiado la zona cercana al punto factible de origen es necesario generar una muestra más o menos uniforme o sobre-disersa que tienda a recorrer todo la región factible y eso no es fácil hacerlo generando sólo unos pocos puntos, más bien hay que generar un paseo aleatorio más o menos largo.

En ​https://www.tol-project.org/ticket/1306#comment:14 se presenta el algoritmo actualmente usado que funciona estupendamente con pocas dimensiones, pero que hay que mejorar cunado aumentan las variables porque se atasca un poco y hay que generar muchos miles de muestras para recorrer medianamente la región factible. Tampoco es que tarde mucho en generar porque sólo tiene que evaluar si el punto es factible, aquí no hay densidades ni Gibbs ni nada de eso, algún producto matriz x vector y cosas así bastante lineales. Comparado con lo que cuesta simular luego la distribución a posteriori no es demasiado, aunque lógicamente sobrecarga la habrá.

Se trata de un algoritmo completamente genérico pero implica la existencia de una función g(x):R^n->R^m tal que la región factible quede definida por las inecuaciones g(x)<=0. Tengo algunas ideas sobre cómo se puede mejorar el rendimiento

Integración de métodos de evaluación de la factibilidad en BSR

Para poder aplicar el algoritmo del punto anterior, o cualquier variante futura del mismo, es necesario introducir en BSR un método que me devuelva g(x) para un punto arbitrario x. Como BSR maneja un modelo dividido en bloques de Gibbs, es cada uno de los bloques el que es capaz de evaluar su parte correspondiente de g(x), a la que podemos llamar g[i](x[i]) para el bloque i-ésimo. Luego el master de BSR se encarga de dividir la entrada x en los bloques x[i] y de concatenar los g[i] de salida.

Esta parte está hecha aunque hace falta mucho chequeo para darlo por bueno, tanto en robustez como en eficacia, especialmente en los bloques no lineales definidos por el usuario que es la parte más complicada y a lo peor es necesario que el filtro incorpore algún método extra para mayor eficacia.

Hay que tener en cuenta que la generación de múltiples cadenas es incompatible con algunas características de BSR

1. La continuación de cadenas es evidentemente imposible si queremos forzar distintos puntos de inicio
   

   Real config::do.resume == 0;
   

2. La simulación parcial no sería en principio incompatible al 100% pero la implementación de los métodos de restricción g(x) serían demasiado complicados de implementar pues dependerían de que partes de x están fijadas y cuáles son libres. Por otra parte no es demasiado problema pues la simulación parcial se usa para reestimar modelos estables o generar previsiones, por lo que se da por supuesto que ya no hay problemas de convergencia.

Generación de múltiples cadenas en BSR

Hasta ahora los pasos para la simulación con BSR eran

1. Crear la definición del modelo como un Struct @BSR.ModelDef modelDef
2. Parametrizar la configuración como una instancia de BysMcmc::@Config config
3. Crear el gestor del modelo BSR como una instancia de la clase BysMcmc::@Cycler cycler
4. Crear el gestor de la estimación BSR como una instancia de la clase BysMcmc::@Estim estim
5. Generar la cadena MCMC llamando al método estim::Run

Si queremos generar múltiples cadenas, los 4 primeros pasos siguen siendo los mismos, salvo que en el paso 2 hay que guardar una copia de config para recuperar la configuración tras el proceso de mezcla

1. Generar unos pocos puntos iniciales muy alejados entre sí
   1. Resetear el modelo llamando al método cycler::Initialize
   2. Generar un paseo aleatorio sobre-disperso en la región factible, es decir, con más probabilidad cerca de la frontera que en el interior.
   3. Extraer unos pocos puntos iniciales lo más alejados entre sí que sea posible partiendo del centro de masas.
   4. Cambiar los parámetros de configuración de las estimaciones locales para que
      1. no haga burn-in ni thinning estándar
         

         Real config::mcmc.burnin    :=  0;
         Real config::mcmc.thinning  :=  0;
         

      2. no haga la diagnosis ni los informes de resultados usuales
         

         Real config::do.report          :=  False;
         Real config::do.eval            :=  False;
         Real config::do.linear.effects  :=  False;
         

2. Para cada punto inicial x0[k] generar de forma independiente una cadena MCMC de la distribución a posteriori
   1. Resetear el modelo llamando al método cycler::Initialize
   2. Asignar el punto inicial con la sentencia Real cycler::_.sampler::setStore(x0[k]);
3. Crear el k-esimo gestor de la estimación BSR como una instancia local de la clase BysMcmc::@Estim local_estim
4. Generar la k-esima cadena local MCMC llamando al método local_estim::Run
5. Devolver la cadena llamando a cycler::loadFullMcmc(?)

Generación en paralelo con BSR