Package BysSampler
Introducción
El paquete BysSampler consiste en una jerarquía de clases para el muestreo bayesiano de cadenas de Markov de una distribución de probabilidad vectorial arbitraria, partiendo únicamente del logaritmo de la función de densidad correspondiente.
Existe una batería de métodos disponibles en la literatura capaces en teoría de dar respuesta a este tipo de problemas en determinadas circunstancias:
- Escalares: Utilizan un método escalar en cada dimensión y el método de Gibbs para vincularlas conjuntamente:
- ARMS: Adaptive Rejection Metropolis Sampling
- SLICE: Slice sampling
- Vectoriales: Simulan directamente la variable aleatoria vectorial como un sólo bloque:
Cada uno de estos algoritmos requiere algún otro tipo de información auxiliar particular de cada caso, pero se intentará que los valores por defecto sean adecuados en la mayoría de las ocasiones, de forma que el usuario no necesite saber demasiado acerca de los detalles de cada uno, salvo en casos en los que la eficiencia o la complejidad del problema lo requieran. El objetivo final sería la implementación de un algoritmo que de forma heurística seleccionara el método más adecuado a un problema en base a sus dimensiones y a un análisis inicial de la función de densidad conjunta.
Los métodos escalares suelen ser más lentos en ejecución y además, si la región de densidad no nula es acotada, requieren la implementación de un método que proporcione los límites de cada variable condicionados a los valores actuales del resto de variables, lo cual puede ser a veces más costoso que la propia simulación.
Los métodos vectoriales en cambio sólo necesitan partir de un vector con densidad no nula. A cambio, la mayoría de ellos precisan de una distribución de probabilidad que genere candidatos plausibles que serán luego aceptados o rechazados según un criterio determinado por cada algoritmo. El sistema proporcionará un generador universal de tipo paseo aleatorio con control de tamaño de paso automático que será suficiente en la mayor parte de las ocasiones. Algunos algoritmos pueden requerir además de alguna otra función que cumpla determinadas condiciones pero para todos ellos se dará al menos una forma por defecto transparente para el usuario.
Cuando el modelo es muy complejo y está compuesto por bloques bien diferenciados lo mejor es aplicar el algoritmo más adecuado a cada uno de ellos y unir los resultados mediante el método de Gibbs.
Jerarquía de clases con base en Class @Generic
La clase @Generic es la clase base de la que heredan todos los simuladores.
Los métodos virtuales puros fundamentales son los que definirán el logaritmo de la función de densidad definir por el usuario para cada problema y el método de muestreo definido por el sistema para cada algoritmo
//Target log-density Real pi_logDens(VMatrix z); //Internal drawer VMatrix _draw(VMatrix x);
Por supuesto el usuario puede heredar sus propias clases que implementen algoritmos no suministrados por el sistema.
Opcionalmente, el usuario podrá implementar el método
//Will be called after each accepted or rejected simulation
Real doAfterDraw(Real void) { True };
que será llamado inmediatamente después de cada llamada al método _draw, lo cual podrá ser utilizado para tareas arbitrarias pero está especialmente indicado para montar una simulación de Gibbs por bloques, pues permite condicionar el resto de bloques al nuevo estado del bloque recién simulado
También tiene una serie de métodos virtuales puros de acceso al vector de variables de forma que se independice la metodología de simulación del mecanismo de almacenamiento de la cadena:
//Number of variables to be simulated Real get.n(Real void); //Returns current drawn values VMatrix get.x(Real void); //Sets current drawn values Real set.x(VMatrix x);
Class @AcceptReject
La clase @AcceptReject : @Generic sirve de base a todos los métodos en los que se genera un candidato generado según una función de transición con densidad Q(x,y) que luego puede ser aceptado o rechazado según cierto criterio.
Class @RandWalk
La clase @RandWalk : @AcceptReject genera los candidatos como un paseo aleatorio en el que en cada paso se da un salto desde el último punto aceptado a otro de su entorno mediante una normal multidimensional Q(x,y) : y ~ N(x,COV).
Class @MetHas
La clase @MetHas : @AcceptReject implementa el más simple de los métodos de simulación vectorial: el método Metropolis-Hastings, que puede ser muy eficiente cuando no hay muchas variables y estas no están muy correlacionadas.
Class @MetHas.RandWalk
Si se hereda de @RandWalk y @MetHas al mismo tiempo se obtiene la modalidad más usual: el método Random Walk Metropolis-Hastings (RWMH)
Class @MulTryMet
La clase @MulTryMet: @AcceptReject implementa el método Multiple Try Metropolis (MTM) basado en generar varios pre-candidatos cada vez y seleccionar uno de ellos de forma también aleatoria pero premiando a los más prometedores según un criterio predefinido mediante una función de pesos producto de la densidad objetivo, densidad de transición de los precandidatos y una función de transición positiva y simétrica llamada lambda(x,y) 
. Con una buena selección del número k de precandidatos y de la función lambda(x,y) se puede mejorar bastante la convergencia y disminuir la autocorrelación de la cadena con respecto al método Metropolis-Hastings, Evidentemente esto sucede a costa de un mayor numero de evaluaciones de la densidad por lo que no es trivial determinar el valor óptimo de k.
Class @MulTryMet.RandWalk
La clase @MulTryMet.RandWalk : @MulTryMet, @RandWalk implementa la modalidad más simple de Multiple Try Metropolis: el método Random Walk Multiple Try Metropolis (RWMTM). Por tratarse de una función de transición simétrica se puede tomar
Class @MulTryMet.RandRay
La clase @MulTryMet.RandRay: @MulTryMet implementa la modalidad Random Ray Multiple Try Metropolis (RRMTM) en la que todos los precandidatos se generan en una misma dirección seleccionada aleatoriamente, de forma simular al método Hit and Run pero evitando la optimización unidimensional que puede resultar demasiado costosa.
Class @Arima
Esta clase implementa los métodos necesarios para manejar un modelo en el que las variables son los coeficientes de los polinomios estacionarios de un modelo ARIMA que hay que contrastar con una serie temporal de ruido ARIMA normal. Los valores iniciales se simulan en cada paso de forma condicionada a los polinomios simulados pero no se almacenan como cadena de Markov. En una próxima versión podría añadirse esta capacidad de forma opcional.
@Arima.MetHas.RandWalk
Simula un modelo ARIMA con el método Random Walk Metropolis-Hastings (RWMH)
@Arima.MulTryMet.RandWalk
Simula un modelo ARIMA con el método Random Walk Multiple Try Metropolis (RWMTM)
@Arima.MulTryMet.RandRay
Simula un modelo ARIMA con el método Random Ray Multiple Try Metropolis (RRMTM)
