[[PageOutline]] = 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 el chequeo de la convergencia * Métodos de '''cadena ú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 [http://www.stat.auckland.ac.nz/~millar/Bayesian/BUGS/coda.pdf CODA] del lenguaje [http://www.r-project.org/ R] * ''Raftery and Lewis'' * ''Geweke'' * ''Heidelberg and Welch'' * Métodos de '''múltiples cadenas''': 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: || [http://uknowledge.uky.edu/cgi/viewcontent.cgi?article=1097&context=gradschool_diss Barton] || Comparisson of two samples by a nonparametric likelihood-ratio test || || [http://alexandria.tue.nl/extra2/200012964.pdf Wiel] || Exact Distributions of Distribution-free Test Statistics || || [http://www.google.es/url?sa=t&source=web&cd=2&ved=0CDUQFjAB&url=http%3A%2F%2Fwww.ism.ac.jp%2Feditsec%2Faism%2Fpdf%2F061_1_0235.pdf&ei=tH6NTsqENPPE4gSD3sSGAQ&usg=AFQjCNEHCQEM5Qc8BZSgU2v9Es6Il7YB4w&sig2=WtNhwn3_dy_wTNCcnpn6YQ Ozturk and Balakrishnan] || Exact two-sample nonparametric test for quantile difference between two populations based on ranked set samples || || [http://www.ge.infn.it/statisticaltoolkit/gof/deployment/userdoc/statistics/documents/Kuiper.ps Kuiper] || Kuiper two sample test || || [http://www.google.es/url?sa=t&rct=j&q=a%20c%2B%2B%20program%20for%20the%20cram%C2%B4er-von%20mises&source=web&cd=2&ved=0CCwQFjAB&url=http%3A%2F%2Fwww.jstatsoft.org%2Fv17%2Fi08%2Fpaper&ei=0B2lTpCkDZGWOvvY2K4N&usg=AFQjCNFaDAewPd7AcY7N3I8dmtQZAPMZIQ&sig2=ySS0WwmbqcButNNDMTDByQ Cramér-von Mises] || A C++ Program for the Cramér-von Mises Two-Sample Test || || [http://www.ruhr-uni-bochum.de/imperia/md/content/mathematik3/publications/dewagvol7.pdf Dette-Wagener-Volgushev] || Nonparametric comparison of quantile curves: a stochastic process approach || || [http://www.stat.tamu.edu/~eparzen/smethmin.pdf Parzen] || Statistical methods mining and nonparametric quantile domain data analysis || || [http://www.bios.unc.edu/~kosorok/909.pdf Kosorok] || Two sample quantile test under general conditions || || [http://www.alglib.net/hypothesistesting/mannwhitneyu.php 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ístico, y 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 [http://www.google.es/url?sa=t&rct=j&q=a%20c%2B%2B%20program%20for%20the%20cram%C2%B4er-von%20mises&source=web&cd=2&ved=0CCwQFjAB&url=http%3A%2F%2Fwww.jstatsoft.org%2Fv17%2Fi08%2Fpaper&ei=0B2lTpCkDZGWOvvY2K4N&usg=AFQjCNFaDAewPd7AcY7N3I8dmtQZAPMZIQ&sig2=ySS0WwmbqcButNNDMTDByQ Cramér-von Mises], que es realmente bastante robusto pero tremendamente lento para cadenas de apenas unos cientos de muestras. El [http://www.alglib.net/hypothesistesting/mannwhitneyu.php 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 [wiki:TolGuiaDelUsuarioSet#FunciónAlgLib.MannWhitneyUtest AlgLib.MannWhitneyUtest] y 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, se 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. 1. 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. 1. En cada intervalo se calcula el test Mann-Whitney U-test sobre las muestras actuales censuradas en dicho intervalo. 1. 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 cuando 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, por lo que habrá que asegurarse de que la configuración es correcta: [[BR]] {{{ #!java Real config::do.resume == 0; }}} 1. 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. 1. Recuperación de problemas numéricos en el bloque lineal cuando la densidad es prácticamente nula que dan lugar a puntos no factibles. Actualmente se reinicializa la cadena con {{{Polytope::FindFeasible}}} lo cual impide aplicar los criterios de convergencia sobre múltiples cadenas. O bien hay que resolver estos problemas o bien hay que utilizar el método MTM que no tiene ese problema. == Generación de múltiples cadenas en BSR == Hasta ahora los pasos para la simulación con BSR eran 1. Preparación del modelo 1. Creación de la definición del modelo como un {{{Struct @BSR.ModelDef modelDef}}} 1. Creación de la configuración como una instancia de {{{BysMcmc::@Config config}}} 1. Creación del gestor del modelo BSR como una instancia de la clase {{{BysMcmc::@Cycler cycler}}} 1. Creación del gestor de la estimación BSR como una instancia de la clase {{{BysMcmc::@Estim estim}}} 1. Generación de la cadena MCMC llamando al método {{{estim::Run}}} Si queremos generar múltiples cadenas, la preparación del modelo sigue siendo igual, salvo que en el paso 1.2 hay que guardar una copia de config para recuperar la configuración tras el proceso de mezcla. 1. Preparación del modelo ... 1. Generación de unos pocos puntos iniciales muy alejados entre sí 1. Reseteo del modelo llamando al método {{{cycler::Initialize}}} 1. Generación de candidatos a puntos iniciales mediante un paseo aleatorio sobre-disperso en la región factible, es decir, con más probabilidad cerca de la frontera que en el interior. 1. Extracción de {{{K}}} puntos iniciales lo más alejados entre sí que sea posible partiendo del centro de masas. 1. Reconfiguración de las estimaciones locales para que 1. no hagan ''burn-in'' ni ''thinning'' estándar[[BR]] {{{ #!java Real config::mcmc.burnin := 0; Real config::mcmc.thinning := 0; }}} 1. no hagan la diagnosis ni los informes de resultados usuales[[BR]] {{{ #!java Real config::do.report := False; Real config::do.eval := False; Real config::do.linear.effects := False; }}} 1. Para cada punto inicial {{{x0[k], k = 1 ... K}}} se genera de forma independiente una cadena MCMC para la distribución a posteriori 1. Reseteo del modelo llamando al método cycler::Initialize 1. Asignación del punto inicial con la sentencia {{{Real cycler::_.sampler::setStore(x0[k]);}}} 1. Creación del k-esimo gestor de la estimación BSR como una instancia local de la clase BysMcmc::@Estim local_estim 1. Generar la k-esima cadena local MCMC llamando al método {{{local_estim::Run}}} 1. Retorno de la cadena llamando a {{{ cycler::loadFullMcmc(?) }}} para evitar pisarla en el disco. 1. Chequeo de convergencia de las cadenas correspondientes a una misma variable. En principio el método es para comparar dos muestras. En principio, para extenderlo a más de dos muestras, habría que aplicarlo a todos los {{{K*(K-1)/2}}} pares posibles pero eso podría suponer demasiado trabajo, por lo podría ejecutarse sobre un conjunto representativo de pares elegidos aleatoriamente, que podrían ser entre {{{2*K}}} y {{{4*K}}} pares. 1. Burn-in individualizado de cada cadena 1. Mezclado de las partes no quemadas por concatenación simple. 1. Reporte de de resultados sobre la cadena conjunta == Generación en paralelo con BSR == En problemas grandes puede ser interesante utilizar una máquina para la generación de cada cadena. Como no es necesario un alto nivel de comunicación entre los distintos procesos involucrados no habría ningún problema en que ésta se desarrollara vía ficheros, pues lo único que habría que hacer sería almacenar el objeto {{{cycler}}} serializado en un fichero {{{cycler.oza}}} que pudieran leer todas las máquinas servidoras, las cuales guardarían las cadenas en ficheros {{{mcmc.k.bbm}}} Pero no sólo sería posible paralelizar las MCMC de la distribución a posteriori, si fuera necesario, también sería posible hacerlo en estos puntos * En el paso 2.2 de generación de candidatos sobre-dispersos a puntos iniciales que luego habría que mezclar en una sola matriz para la extracción del punto 2.3. * En el paso 5 de chequeo de la convergencia se pueden dividir los pares de cadenas seleccionadas para que el chequeo se hiciera en varias máquinas a la vez. Sólo habría que recolectar los resultados finales en un lugar preestablecido. === Monitorización de convergencia en paralelo === En el caso de generación en paralelo existiría la posibilidad de que el chequeo de la convergencia lo ejecutara un proceso aparte desde el ''master'' que fuera leyendo cada cierto tiempo las cadenas generadas hasta el momento por cada ''slave'' en sus respectivos archivos {{{mcmc.bbm}}}. Una vez leídas se procedería a chequear la convergencia y en caso positivo se ejecutaría un test de tipo Raftery sobre la cadena conjunta ya quemada, para establecer si la longitud total de la muestra es suficientemente significative y, si efectivamente es así, comunciar a los procesos de alguna forma que ya pueden parar de simular. Luego sólo restaría volver a leer las cadenas para incluir las simulaciones realizadas entre tanto para finalizar con los últimos dos puntos de mezclado y reporte.