| 4 | |
| 5 | = Package BysSampler = |
| 6 | |
| 7 | == Introducción == |
| 8 | |
| 9 | 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. |
| 10 | |
| 11 | 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: |
| 12 | |
| 13 | * '''Escalares''': Utilizan un método escalar en cada dimensión y el método de Gibbs para vincularlas conjuntamente: |
| 14 | * '''ARMS''': Adaptive Rejection Metropolis Sampling |
| 15 | * '''SLICE''': [http://en.wikipedia.org/wiki/Slice_sampling Slice sampling] |
| 16 | * '''Vectoriales''': Simulan directamente la variable aleatoria vectorial como un sólo bloque: |
| 17 | * '''MH''': [http://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm Metropolis-Hastings] |
| 18 | * '''MTM''': [http://en.wikipedia.org/wiki/Multiple-try_Metropolis Multiple Try Metropolis] |
| 19 | * '''GMTM''': [http://jmlr.csail.mit.edu/proceedings/papers/v9/pandolfi10a/pandolfi10a.pdf Generalized Multiple Try Metropolis-Hastings] |
| 20 | * '''CGMC''': [http://www.google.es/url?sa=t&source=web&cd=2&sqi=2&ved=0CCEQFjAB&url=http%3A%2F%2Fwww.people.fas.harvard.edu%2F~junliu%2FTechRept%2F98folder%2Fmtm.ps.gz&ei=szXQTIr9HIzqOZPlgdAE&usg=AFQjCNEYiDDLHk2M5Zmn5Wr23aKqS2_7LA&sig2=MpGM_qLG9by78g_2mq-VbA Conjugate gradient Monte Carlo] |
| 21 | * '''RRMC''': [http://books.google.es/books?id=Dk_ou-gqnHQC&pg=PA134&lpg=PA134&dq=random+ray+monte+carlo+Liu&source=bl&ots=UfoWWBBWuh&sig=K3q23CbKBkf1CWYmx2QLcCIqPDg&hl=es&ei=mkTRTIXJOZHQjAfC7b2oDA&sa=X&oi=book_result&ct=result&resnum=1&ved=0CBgQ6AEwAA#v=onepage&q=random%20ray%20monte%20carlo%20Liu&f=false Random Ray Monte Carlo] |
| 22 | |
| 23 | 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. |
| 24 | |
| 25 | 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. |
| 26 | |
| 27 | 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. |
| 28 | |
| 29 | 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. |
| 30 | |
| 31 | == Jerarquía de clases con base en {{{Class @Generic}}} == |
| 32 | |
| 33 | La clase {{{@Generic}}} es la clase base de la que heredan todos los simuladores. |
| 34 | |
| 35 | 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 |
| 36 | {{{ |
| 37 | #!cpp |
| 38 | //Target log-density |
| 39 | Real pi_logDens(VMatrix z); |
| 40 | //Internal drawer |
| 41 | VMatrix _draw(VMatrix x); |
| 42 | }}} |
| 43 | Por supuesto el usuario puede heredar sus propias clases que implementen algoritmos no suministrados por el sistema. |
| 44 | |
| 45 | Opcionalmente, el usuario podrá implementar el método |
| 46 | {{{ |
| 47 | #!cpp |
| 48 | //Will be called after each accepted or rejected simulation |
| 49 | Real doAfterDraw(Real void) { True }; |
| 50 | }}} |
| 51 | 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 |
| 52 | |
| 53 | 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: |
| 54 | {{{ |
| 55 | #!cpp |
| 56 | //Number of variables to be simulated |
| 57 | Real get.n(Real void); |
| 58 | //Returns current drawn values |
| 59 | VMatrix get.x(Real void); |
| 60 | //Sets current drawn values |
| 61 | Real set.x(VMatrix x); |
| 62 | }}} |
| 63 | |
| 64 | === {{{Class @AcceptReject }}} === |
| 65 | 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. |
| 66 | |
| 67 | ==== {{{Class @RandWalk }}} ==== |
| 68 | |
| 69 | 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). |
| 70 | |
| 71 | === {{{Class @MetHas }}} === |
| 72 | La clase {{{@MetHas : @AcceptReject}}} implementa el más simple de los métodos de simulación vectorial: el método [http://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm Metropolis-Hastings], que puede ser muy eficiente cuando no hay muchas variables y estas no están muy correlacionadas. |
| 73 | |
| 74 | ==== {{{Class @MetHas.RandWalk }}} ==== |
| 75 | 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) |
| 76 | |
| 77 | |
| 78 | === {{{Class @MulTryMet}}} === |
| 79 | La clase {{{@MulTryMet: @AcceptReject}}} implementa el método [http://en.wikipedia.org/wiki/Multiple-try_Metropolis 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) [[LatexEquation( \lambda \left( x,y \right) $$)]] . 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 [http://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm 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. |
| 80 | |
| 81 | ==== {{{Class @MulTryMet.RandWalk}}} ==== |
| 82 | La clase {{{@MulTryMet.RandWalk : @MulTryMet, @RandWalk}}} implementa la modalidad más simple de [http://en.wikipedia.org/wiki/Multiple-try_Metropolis 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 |
| 83 | |
| 84 | ==== {{{Class @MulTryMet.RandRay}}} ==== |
| 85 | 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 [http://oai.cwi.nl/oai/asset/1755/1755A.pdf Hit and Run] pero evitando la optimización unidimensional que puede resultar demasiado costosa. |
| 86 | |
| 87 | === {{{Class @Arima}}} === |
| 88 | 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. |
| 89 | |
| 90 | ==== {{{@Arima.MetHas.RandWalk}}} ==== |
| 91 | Simula un modelo ARIMA con el método Random Walk Metropolis-Hastings (RWMH) |
| 92 | |
| 93 | ==== {{{@Arima.MulTryMet.RandWalk}}} ==== |
| 94 | Simula un modelo ARIMA con el método Random Walk Multiple Try Metropolis (RWMTM) |
| 95 | |
| 96 | ==== {{{@Arima.MulTryMet.RandRay}}} ==== |
| 97 | Simula un modelo ARIMA con el método Random Ray Multiple Try Metropolis (RRMTM) |