[[PageOutline]] = Generación de candidatos factibles en un politopo = == Descripción del problema == Tenemos que generar una cadena de Markov de muestras que se distribuyan con cierta función de la que se conoce su log-densidad salvo una constante, [[LatexEquation( \ln \pi \left( x \right) )]] a la par que cumplen un conjunto de restricciones de desigualdad lineal [[LatexEquation( A x \le a \wedge x\in\mathbb{R}^{n} )]] [[LatexEquation( x\in\mathbb{R}^{n} )]] [[LatexEquation( a\in\mathbb{R}^{r} )]] [[LatexEquation( A\in\mathbb{R}^{r\times n} )]] Lo único que deben cumplir las restricciones es que el politopo generado por ellas, o sea, el conjunto de puntos que las cumplen, tenga medida no nula en [[LatexEquation( \mathbb{R}^{n} )]], para que tengan sentido los algoritmos que se van a plantear. La matriz [[LatexEquation(A)]] puede ser singular o tener cualquier forma con tal cumpla esa condición. En particular puede haber restricciones no activas, es decir, que se pueden quitar sin que cambie en absoluto la forma del politopo generado. Para obtener un elemento de la cadena dado el anterior se utilizará algún método similar al Metropolis-Hastings, como podría ser el Multiple Try Metropolis. Todos ellos necesitan un método de generación aleatoria de candidatos o precandidatos que luego serán aceptados o no. Nos caben dos posibilidades a la hora de diseñar el generador de candidatos: 1. utilizar un generador de candidatos que cumpla las restricciones por construcción. 1. usar un generador no restringido y luego rechazar los candidatos no factibles. La ventaja de la segunda es que el generador puede ser simétrico y evitarse el cálculo de la verosimilitud del candidato, pero el problema es que puede ser que tarde mucho en encontrar uno factible si el punto actual está demasiado cerca de la frontera, lo cual será muy habitual si el punto de máxima verosimilitud se encuentra fuera del politopo definido por las anteriores inecuaciones. En la siguiente figura se observa un caso en el que el punto máximo-verosímil sin restricciones (en verde) se encuentra fuera del politopo (en gris), por lo que el punto máximo-verosimil restringido (en rojo) se encuentra en la frontera del politopo. Las elipses verdes son líneas isoprobables, es decir con la misma densidad por lo que ese punto se encuentra en la intesección con el politopo de la mayor elipse externa al mismo. La cadena tenderá a estar cerca de ese punto y por tanto cerca de la frontera. Es evidente que cualquier entorno de generación simétrica de candidatos (en rosa), tendrá como máximo la mitad de área en la zona factible, lo cual es perfectamente asumible. [[Image(source:/tolp/OfficialTolArchiveNetwork/BysSampler/doc/image/RandWalk.InPolytope.chart.1026074815.png)]] Es claro que no es preciso que el punto se encuentre en la frontera para que sus entornos tengan incluso menos de la mitad de sus puntos dentro del politopo, como ocurre con el entorno amarillo de la figura anterior. La pregunta es ¿qué pasaría si no estuviera en la frontera por inclumplir sólo una de las restricciones sino que incumpliera varias al mismo tiempo? Pues que la probabilidad de generar un punto factible decrecería exponencialmente, dependiendo también de los ángulos formados por los hiperplanos que definen cada restricción. Por ejemplo, si todos los ángulos son de 90 grados cada restricción reduce a la mitad el espacio factible, así que con [[LatexEquation(k)]] restricciones el ratio de aceptación se reduciría a [[LatexEquation(2^{-k})]] Y eso ya no es asumible a nada que se tengan 3 ó mas restricciones incumplidas. En dos dimensiones sólo se pueden cruzar dos restricciones activas al mismo tiempo pero es fácil de extrapolar lo que ocurriría en espacios de dimensiones más altas: [[Image(source:/tolp/OfficialTolArchiveNetwork/BysSampler/doc/image/RandWalk.InPolytope.chart.726488582.png)]] En la práctica se observa que cuando hay muchas restricciones la cadena tiende a colapsar en un punto de la frontera, lo cual impide que haya una convergencia real en un intervalo de tiempo razonable. Si se utiliza un generador de tipo Gibbs le pasa exactamente lo mismo, pero la consecuencia es que el método directamente fracasa y no es capaz de simular el siguiente punto. Se hace necesario por lo tanto disponer de generadores de candidatos que por construcción estén siempre dentro del politopo, es decir generadores de candidatos factibles. Tal generador no podrá ser simétrico por lo que será necesario no sólo poder generar muestras sino también calcular su densidad de una forma eficiente. == Definiciones previas == Dado un punto [[LatexEquation(x)]] interior al politopo [[LatexEquation(A x \le a)]] la distancia al [[LatexEquation(i)]]-ésimo hiperplano viene dada por la fórmula [[LatexEquation( \left\langle x,\mathcal{H}\left(A_i,a_i\right)\right\rangle = \frac{\left|\left(\overset{n}{\underset{j=1}{\sum}}A_{ij}x_{j}\right)-a_{i}\right|}{\sqrt{\overset{n}{\underset{j=1}{\sum}}A_{ij}^{2}}} = \frac{\left|A_{i}x-a_{i}\right|}{\sqrt{A_{i}^{T}A_{i}}} \forall i=1\ldots r)]] Llamaremos distancia de un punto interior a la frontera del politopo a la menor de las distancias a cada uno de sus hiperplanos, la cual por ser el punto interior ha de ser forzosamente no negativa [[LatexEquation( \left\langle x,\mathcal{P}\left(A,a\right)\right\rangle = \underset{i=1\ldots r}{min}\left\{ \left\langle x,\mathcal{H}\left(A_i,a_i\right)\right\rangle \} \ge 0 )]] Obsérvese que los hiperplanos inactivos nunca darán la distancia mínima para un punto factible por lo que efectivamente no influyen en absoluto en los cálculos, salvo que aumentan el tiempo de cálculo de forma espúrea. Sin embargo no es en absoluto trivial determinar si una restricción es inactiva así que será responsabilidad del analista no introducir demasiadas. Por este motivo no es posible calcular la distancia de un punto exterior al politopo si antes no se excluyen las restricciones inactivas. == Búsqueda de un punto inicial estrictamente interior == Para encontrar un punto estrictamente interior hay que encontrar un punto que cumpla las restricciones artificialmente aumentadas [[LatexEquation(A x \le a-\delta \wedge \delta>0 )]] que producirán un politopo reducido estrictamente interior al original. Es posible configurar el vector [[LatexEquation(\delta )]] para que la distancia mínima a la frontera sea mayor que una cota mínima [[LatexEquation(\rho_0>0 )]] predefinida por el usuario [[LatexEquation(\delta_i = \rho_0 \sqrt{A_{i}^{T}A_{i}} )]] De esta forma si un punto está en el [[LatexEquation(i)]]-ésimo hiperplano del politopo reducido [[LatexEquation(A_i x = a_i-\delta_i )]] entonces la distancia al politopo original será [[LatexEquation(\frac{\delta_i}{\sqrt{A_{i}^{T}A_{i}}} = \rho_0 )]] El mejor punto inicial en estas condiciones sería el resultado de optimizar la función de log-densidad [[LatexEquation( \min \ln \pi \left( x \right) )]] sujeto a [[LatexEquation(A x \le a-\delta )]] [[Image(source:/tolp/OfficialTolArchiveNetwork/BysSampler/doc/image/RandWalk.InPolytope.chart.1051120617.png)]] Si el proceso de optimización resulta muy lento hay que recordar que tampoco es necesario partir del óptimo, pues cualquier punto con log-densidad significativa podría servir, aunque eso a su vez inducirá un periodo de convergencia mayor de la cadena de Markov. Habría que combinar unas condiciones de parada razonables como pudiera ser un tiempo máximo y un valor mínimo para la log-densidad. == Paseo aleatorio hiperesférico == Se propone en primer lugar utilizar un generador con distribución uniforme en una hiperesfera centrada en el último punto generado y que esté incluida estrictamente en el politopo lo cual será cierto si el radio es menor que la distancia del punto a la frontera: [[LatexEquation( \begin{array}{c} y=x + h u\\ A y \leq a\\ \left\Vert u\right\Vert = 1 \\ 0 < h \le \rho < \left\langle x,\mathcal{P}\left(A,a\right)\right\rangle \end{array} )]] Al igual que se hace con los generadores simétricos habituales es necesario fijar un tamaño de paso máximo [[LatexEquation(s)]] para que se generen puntos no demasiado lejanos del actual, de forma que el ratio de aceptación sea razonable. De esta forma se podría definir el radio de la hiperesfera como [[LatexEquation(\rho=\rho\left(x; s, \lambda\right)=min\left\{ s,\left\langle x, \lambda \mathcal{P}\left(A,a\right)\right\rangle \right\} \wedge 0<\lambda<1)]] Para poder muestrear puntos cercanos a la frontera sería conveniente que [[LatexEquation(\lambda)]] fuera cercano a 1, pero no tanto como para que los candidatos se quedaran demasiado cerca de ésta, obligando a que el radio se estrechara cada vez más. [[Image(source:/tolp/OfficialTolArchiveNetwork/BysSampler/doc/image/RandWalk.InPolytope.chart.726488582_b.png)]] Obsérvese como este método tendrá siempre dificultades al acercarse demasiado a la frontera, y especialmente si se acerca a un vértice demasiado agudo. === Función de densidad === La densidad del generador será porporcional al inverso del volumen de la hiperesfera y su logaritmo será, salvo una constante [[LatexEquation( \ln Q\left(x,y\right) = n \ln \rho\left(x; s, \lambda\right))]] La densidad depende por tanto del radio, el cual depende de la posición del punto de partida respecto a las fronteras del politopo y no es por tanto un generador simétrico necesariamente, auqnue sí puede serlo eventualmente si la distancia del punto a la frontera es menor que el tamaño de paso actual. En realidad podría usarse cualquier otra densidad truncada en la hiperesfera, como por ejemplo una multinormal estandarizada, pero eso puede complicar demasiado el cálculo de la función de densidad, la cual habrá que calcular dos veces por cada candidato o precandidato generado. === Función generatriz === Para generar un candidato [[LatexEquation(y)]] a partir del actual [[LatexEquation(x)]] con esta distribución se procederá como sigue 1. Si el último punto generado es estrictamente interior al politopo y su distancia a la frontera es mayor que cierta cota positiva prefijada [[LatexEquation( \rho_0 > 0)]] entonces usamos dicho punto como origen de candidatos 1. Primero se simulará un vector multinormal estandarizado [[BR]] [[BR]] [[LatexEquation(w\sim N\left(0,I\right)\in\mathbb{R}^{n} )]] [[BR]] [[BR]] 1. Dividiendo ese vector por su norma euclídea se obtendrá una dirección unitaria uniforme en la frontera de la hiperesfera de radio 1 centrada en el origen [[BR]] [[BR]] [[LatexEquation( u=\frac{w}{\left\Vert w\right\Vert _{2}} )]] [[BR]] [[BR]] 1. Luego se multiplica ese vector por un radio [[LatexEquation(h)]] con densidad directamente proporcional a la hipersuperficie de dimensión n-1 de la hiperesfera correspondiente[[BR]] [[BR]] [[LatexEquation( \begin{array}{c} F\left(h\right)=\intop_{0}^{h}c\cdot t^{n-1}\mathbf{d}t=\left.\frac{c}{n}t^{n}\right]_{0}^{h}=\frac{c}{n}h^{n}\\ F\left(\rho\right)=\frac{c}{n}\rho^{n}=1\Rightarrow c=n\rho^{-n}\\ F\left(h\right)=\left(\frac{h}{\rho}\right)^{n}\end{array} )]] [[BR]] [[BR]] Para ello se genera una uniforme unitaria y se calcula la inversa de la función de distribución en ella [[BR]] [[BR]] [[LatexEquation( \gamma \sim \left[0,1\rigt] )]] [[BR]] [[BR]] [[LatexEquation( h = F^{-1}\left(\gamma\right) = \rho \gamma^{\frac{1}{n}})]] [[BR]] [[BR]] 1. Si el punto está en la frontera o a una distancia menor de [[LatexEquation( \rho_0 )]] entonces se retomará el último punto interior que cumplía dicha condición como origen de los candidatos, y se ejecuta desde el principio el proceso, pero luego se compararán las probabilidades de los candidatos con las del último punto generado, para asegurar así la convergencia, según el principio básico del [http://en.wikipedia.org/wiki/Detailed_balance balance-detallado]. Podríamos decir que se trata de un paseo aleatorio con retroceso, pues en ocasiones el origen de candidatos puede retroceder y no coincidir con el último pnto generado. Esto es perfectamente lícito pues no hay ninguna regla que obligue a que la distribución de los candidatos sea un paseo aleatorio. De hecho podría utilizarse perfectamente un generador de candidatos con origen fijo, aunque sólo funcionaría en la práctica si se conoce una buena aproximación de la función de densidad que sea fácil de generar y calcular. Aún así, como el comportamiento mayoritario será el de un paseo aleatorio mantendremos la nomenclatura con todas las reservas enunciadas. == Paseo aleatorio radial asimétrico == A continuación de define otra forma de generador factible similar a la anterior pero que podría adaptarse mejor cuando el punto se acerca demasiado a la frontera, e incluso es válida si está en la propia frontera. Además se puede utilizar con restricciones de desigualdad no necesariamente lineales [[LatexEquation( g_i \left( x \right) \le 0 \forall i=1 \dots r )]] Se trataría de tomar una dirección unitaria aleatoria siguiendo los mismos dos primeros pasos del punto anterior, y luego encontrar un intervalo en esa dirección que incluya al punto actual pero no necesariamente centrado en él, para generar el candidato uniforme en dicho intervalo [[LatexEquation( y = x + h u )]] [[LatexEquation(u \wedge \left\Vert u\right\Vert = 1)]] [[LatexEquation(h \sim U\left[h^{-},h^{+}\right] \wedge h^{-} \le 0 \le h^{+} )]] [[Image(source:/tolp/OfficialTolArchiveNetwork/BysSampler/doc/image/RandWalk.InPolytope.chart.726488582_c.png)]] === Función generatriz === Los pasos para muestrear un candidato serían los siguientes 1. Primero se simulará un vector multinormal estandarizado [[BR]] [[BR]] [[LatexEquation(w\sim N\left(0,I\right)\in\mathbb{R}^{n} )]] [[BR]] [[BR]] 1. Dividiendo ese vector por su norma euclídea se obtendrá una dirección unitaria uniforme en la frontera de la hiperesfera de radio 1 centrada en el origen [[BR]] [[BR]] [[LatexEquation( u=\frac{w}{\left\Vert w\right\Vert _{2}} )]] [[BR]] [[BR]] 1. Si los puntos [[LatexEquation( x - s u)]] y [[LatexEquation( x + s u)]] son ambos factibles entonces se toma [[LatexEquation( h^{-} = -s)]] y [[LatexEquation( h^{+} = +s )]] 1. Si [[LatexEquation( x - s u)]] no es factible habría que buscar mediante un método univariante, como el de la bisectriz o el de Fibonacci, el mínimo valor de [[LatexEquation( h^{-} )]] para el que el punto es factible, y que debe estar obligatoriamente en el intervalo [[LatexEquation(\left[-s,0\right] )]] 1. Análogamente, si [[LatexEquation( x + s u)]] no es factible habría que buscar el máximo valor de [[LatexEquation( h^{-} )]] para el que el punto es factible, y que debe estar obligatoriamente en el intervalo [[LatexEquation(\left[0,s\right] )]] 1. Si el intervalo es propio, es decir, si [[LatexEquation( h^{-} < h^{+})]], entonces se puede generar la uniforme anterior. 1. Si el intervalo no es propio entonces el punto generado es obligatoriamente el mismo en el que se estaba. Para evitar que se quede demasiado tiempo en el mismo punto, en la siguiente iteración se tomará el último punto con intervalo propio como origen de los candidatos y se repite desde el principio generando una nueva dirección unitaria aleatoria, pero se compararán las probabilidades de los candidatos con las del último punto generado, para asegurar así la convergencia. Para evitar intervalos demasiado pequeños se puede relajar la condición del cuarto paso a que su longitud sea mayor que cierta cota mínima predefinida [[LatexEquation( h^{+} - h^{-} > \h_0 > 0)]], === Función de densidad === La función de densidad en este caso es simplemente proporcional al inverso de la longitud del intervalo, luego su logaritmo salvo una constante será [[LatexEquation( \ln Q\left(x,y\right) = -\ln\left(h^{+}-h^{-}\right) )]] Los pasos del cálculo del intervalo serán ahora los siguientes: 1. Primero se toma el vector de dirección que une los puntos [[BR]] [[BR]] [[LatexEquation(w = y-x )]] [[BR]] [[BR]] 1. Dividiendo ese vector por su norma euclídea se obtendrá la dirección unitaria de la hiperesfera de radio 1 centrada en el origen [[BR]] [[BR]] [[LatexEquation( u=\frac{w}{\left\Vert w\right\Vert _{2}} )]] [[BR]] [[BR]] 1. Se repiten los pasos 3, 4 y 5 realizados en la generación para calcular el intervalo Está claro que la longitud del intervalo depende del punto de origen por lo que el generador no es simétrico. == Paseo aleatorio mixto == El paseo aleatorio hiperesférico tiene la ventaja de ser muy rápido pues sólo se necesita la distancia a la frontera, pero tiene el inconveniente que no es capaz de acercarse demasiado a la frontera y menos aún a los vértices donde a veces podría concentrarse la verosimilitud. El paseo aleatorio radial asimétrico tiene la ventaja de que puede acercarse mucho más a la frontera, de hecho puede simular puntos fronterizos, pero resulta mucho más caro de evaluar, pues puede precisar resolver hasta dos problemas de optimización univariante para cada generación y para cada evaluación de la log-densidad. Si dos métodos de generar cadenas de Markov cumplen la propiedad del [http://en.wikipedia.org/wiki/Detailed_balance balance-detallado] entonces cualquier mezcla de ellas también la cumple, por lo que es perfectamente lícito alternar ambos métodos según convenga. En este caso podría utilizarse el método hiperesférico cuando el punto se encuentre a una distancia de la frontera mayor o igual que [[LatexEquation( \rho_0 )]] y el método radial asimétrico en otro caso. == Bibliografía == Hay poquísima literatura acerca de simulación con restricciones: [http://www.mims.manchester.ac.uk/research/probability-statistics/research-reports/psrr19-2006.pdf Optimal Scaling for Random Walk Metropolis on Spherically Constrained Target Densities]