wiki:OfficialTolArchiveNetworkBysSamplerPostProccess

Context Navigation

close Warning: Can't synchronize with repository "(default)" (/var/svn/tolp does not appear to be a Subversion repository.). Look in the Trac log for more information.

Version 19 (modified by Víctor de Buen Remiro, 14 years ago) (diff)
--

Post-procesado de cadenas de Markov

Post-procesado de cadenas de Markov

Introducción

Los métodos tradicionales de post-procesado de cadenas de simulación basados en el burn-in y el thinning son demasiado arbitrarios para poder parametrizarlos de forma automática sin intervención del usuario. Además no solucionan uno de los principales problemas de los métodos accept-reject basados en paseos aleatorios, que es la alternacia de fases con exceso de repeticiones por usar un tamaño de paso demasiado grande con otras fases en las que se avanza poco en cada iteración por todo lo contrario.

Las cadenas simuladas con BysSampler cuentan con una ventaja adicional al conocerse el logaritmo de la verosimiltud de cada punto muestral, salvo una constante aditiva, pues esto permite contrastarla directamente con la masa local empírica, es decir el número de puntos que han sido generados en sus cercanías.

En una cadena perfectamente muestreada el número de puntos generados en torno a un punto dado debería ser proporcional a la verosimilitud media alrededor de dicho punto. Esto permite diseñar un criterio completamente objetivo para eliminar puntos de zonas supoerpobladas e incluso sustituirlos por puntos en otras zonas inframuestreadas.

Diseño de entornos locales solapados

El método propuesto será utilizar el algoritmo KNN que está disponible dentro del paquete TOL MatQuery, para encontrar los vecinos más próximos de cada punto de la muestra de tamaño $S$ . Como en los métodos de simulación tipo accept-reject hay por definición puntos repetidos, para que el algoritmo tenga sentido habría que tomar los $S' < S$ puntos únicos

$x_{i} \wedge i=1 \dots S'$

y llamar

$s_{i} \wedge i=1 \dots S'$

al número de veces que aparece cada uno en la muestra. Obviamente, la suma de los números de apariciones da el tamaño muestral

$S=\underset{i=1}{\overset{S'}{\sum}}s_{i}$

No está claro por el momento cuál podría ser el criterio para seleccionar el número $k$ de vecinos en cada entorno, pero debe ser en cualquier caso un número bastante pequeño en relación con $S'$ pues se trata de observar el comportamiento a nivel local. También ha de ser bastante pequeño en términos absolutos porque la complejidad del algoritmo KNN crece cuadráticamente con el tamaño del vecindario. Tampoco puede ser excesivamente pequeño porque entonces quedarían amplias zonas del espacio sin recubrir por no estar lo suficientemente cerca de ninguno de los puntos muestreados. Lo ideal sería tomar un sistema de entornos que recubriera de forma conexa y compacta la muestra, aunque eso no parece facil de comprobar a primera vista. Tampoco es en principio obligatorio tomar el mismo número de puntos en cada entorno aunque por el momento supondremos que es así. Sean los $k$ puntos muestrales vecinos de $x_{i}$ en orden de proximidad al mismo

$y_{i,j} \wedge i=1 \dots S \wedge j=1 \dots k$

Sea la distancia euclídea del punto $x_{i}$ a su $k$ -ésimo vecino más próximo

$r_{i,k}=\sqrt{\underset{d=1}{\overset{n}{\sum}}\left(y_{i,k,d}-x_{i,d}\right)^{2}}$

En cada punto muestral definiremos pues el entorno local como la hiperesfera de radio $r_{i,k}$ y centro $x_{i}$

$\Omega_{i}=\left\{ y\in\mathbb{R}^{n}\left|\left\Vert y-x_{i}\right\Vert ^{2}\leq r_{i,k}\right.\right\}$

Se trata de ver qué entornos están superpoblados y cuáles están inframuestreados, para lo cual deberemos estudiar la distribución de probabilidad del tamaño muestral incluido dentro de cada uno.

Distribución del tamaño muestral local

Así las cosas tenemos que el tamaño muestral local, es decir, el número total de puntos muestrales en $\Omega_{i}$ , es

$h_{i} = k+s_i$

cantidad que se distribuye como una binomial

$\eta_{i} = B\left(S, p_i\rigtht)$

donde $p_i$ es la probabilidad de la hiperesfera, es decir, la integral de la función de densidad en cada entorno local

$p_{i}=\int_{\Omega_{i}}\pi\left(y\right)\mathrm{d}y$

Aproximación de la probabilidad del entorno local

La anterior integral sería algo muy costoso de evaluar así que hay que aproximarla por el método de Montecarlo, como el producto de la media de las verosimilitudes en un conjunto de puntos del entorno por el hipervolumen del mismo

$p_{i}=\frac{1}{N}\underset{j=1}{\overset{N}{\sum}}\pi\left(z_{i,j}\right)V\left(\Omega_{i}\right)\wedge z_{i,j}\in\Omega_{i}\forall j=1\ldots N\wedge i=1\ldots S'$

El error en este tipo de aproximaciones decrece proporcionalmente a la raíz del número de puntos en el que evaluamos la verosimilitud pero sólo tenemos $h_{i}$ puntos interiores. Por otra parte tampoco conocemos la verosimilitud sino una función porporcional a la misma. Es decir, lo único conocemos sin coste adicional es el logaritmo de la verosimilitud, salvo una constante $\lambda_0$ desconocida, evaluado en cada uno de los puntos muestrales, es decir, conocemos

$\ln\pi_{i}=\ln\pi\left(x_{i}\right)+\lambda_0$

$\ln\pi_{i,j}=\ln\pi\left(y_{i,j}\right)+\lambda_0$

Como evaluar la función demasiadas veces podría ser muy costoso la única forma de mejorar la aproximación de la integral es por interpolación, concretamente mediante el método de Sheppard de ponderación inversa a la distancia que es muy eficiente pues no requiere de ninguna evaluación extra. Esto es especialemente recomendable si existen grandes diferencias en las verosimilitudes de los distintos puntos del vecindario. La interpolación será mejor realizarla en términos logarítmicos pues eso suavizará la función.

Para ello generaremos $N$ puntos con distribución uniforme en la hiperesfera

$\left\Vert z_{i,j}-x_{i}\right\Vert ^{2}\leq r_{i,k} \wedge j=1 \dots N$

y calculamos la aproximación del logaritmo de la verosimilitud en cada uno de ellos mediante la fórmula de ponderación de Sheppard

$\ln\tilde{\pi}\left(z\right)=\frac{\underset{j=0}{\overset{k}{\sum}}w_{j}\left(z\right)\ln\pi_{i,j}}{\underset{j=0}{\overset{k}{\sum}}w_{j}\left(z\right)}\wedge y_{i,0}=x_{i}\wedge\pi_{i,0}=\pi_{i}\wedge w_{j}\left(z\right)=\left\Vert z-y_{i,j}\right\Vert ^{-2}$

Si el número $k+1$ de puntos básicos de la interpolación $y_{i,j} \wedge j=0 \dots k$ , es menor o igual que la dimensión del espacio $n$ la anterior aproximación no es lo suficientemente flexible, pero eso es facil de evitar añadiendo más puntos básicos, simplemente tomando sus vecinos, los vecinos de sus vecinos y así sucesivamente hasta que haya al menos $n$ puntos básicos distintos. Gracias al algoritmo KNN esto no supondrá apenas ningún sobrecoste. Como el volumen de la hiperesfera es proporcional a

$r^n_{i,k}$

nos queda la fórmula de aproximación

$\ln p_{i}=\mu_{0}+\nu_{i}$

$\nu_{i}=\ln\left(\frac{1}{N}\underset{j=1}{\overset{N}{\sum}}\tilde{\pi}_{i}\left(z_{i,j}\right)\right)+n\ln r_{i,k}$

Puesto que una probabilidad ha de ser menor que 1 su logaritmo es siempre negativo, por lo que tenemos una cota para la constante

$\mu_{0}<-\underset{i=S'}{\max}\left\{ \nu_{i}\right\}$

Necesitamos un criterio razonable para establecer cuál debe ser el valor de la constante $\mu_{0}$ para poder continuar con los cálculos, y dado que conocemos la forma de la distribución podemos encontrar su valor máximo-verosímil.

Verosimilitud del parámetro

La probabilidad de que el número de puntos que caen dentro de la hiperesfera sea exactamente $h$ para la binomial definida anteriormente es

$P_i = \mathrm{Pr}\left[\eta_{i}=h_{i}\right]=\left(\begin{array}{c}S\\h_{i}\end{array}\right)p_{i}^{h_{i}}\left(1-p_{i}\right)^{S-h_{i}}$

y el logaritmo de dicha probabilidad será

$\ln\left(P_{i}\right)=\ln\left(\begin{array}{c}S\\h_{i}\end{array}\right)+h_{i}\ln p_{i}+\left(S-h_{i}\right)\ln\left(1-p_{i}\right)$

La verosimilitud de $\mu_0$ dada la muestra observada, bajo la hipótesis de independencia entre los distintos entornos, será el productorio de las probabilidades del número de puntos efectivamente encontrados en cada uno, teniendo en cuenta que los puntos repetidos deben multiplicarse tantas veces como aparecen. La expresión de su logaritmo será por tanto

$L\left(\mu_{0}\right)=\underset{i=1}{\overset{S'}{\sum}}s_{i}\ln\left(P_{i}\right)$

En realidad los entornos cercanos no pueden ser independientes entre sí, pues de hecho comparten puntos, pero en primera instancia daremos por buena la hipótesis de independencia, simplemente por comodidad y porque no está claro que sea demasiado importante el efecto de la dependencia.

Así pues tendremos el problema de optimización univariante

$\underset{\mu_{0}}{\max}\left\{ L\left(\mu_{0}\right)\right\}$

Sujeto a

$\mu_{0}<-\underset{i=S'}{\max}\left\{ \nu_{i}\right\}$

Resolviendo este problema obtenemos el valor de $\mu_{0}$ que nos permite completar el resto de cálculos.

Función de distribución

La probabilidad de que el número de puntos que caen dentro de la hiperesfera sea mayor o igual que $h$ se calcula mediante la función beta incompleta

$\mathrm{Pr}\left[\eta_{i}\leq h_{i}\right]=I_{1-p_i}\left(S-h_{i},h_{i+1}\right)$

Filtrado de la superpoblación

En los entornos en los que la probabilidad de exceso de muestra,

$\mathrm{Pr}\left[\eta_{i}\leq h_{i}\right]$

sea muy alta habrá que eliminar puntos, empezando por los repetidos, hasta que se entre dentro de un margen razonable ...

Estrategia de recolonización

En los entornos en los que la probabilidad de defecto de muestra,

$\mathrm{Pr}\left[\eta_{i}\geq h_{i}\right]$

sea muy alta habrá que añadir más mediante un mecanismo que asegure la convergencia. Una posibilidad sería continuar el mismo método utilizado en la generación de la muestra analizada comenzando por los puntos centrales de los entornos más despoblados hasta compensar la masa faltante. Pero dada la información acumulada sería quizás más razonable utilizar un generador de candidatos con media en los puntos centrales en lugar de usar un paseo aleatorio. Incluso se podría usar el método de ensayo múltiple generalizado usando como precandidatos los mismos puntos generados anteriormente para la aproximación de la integral.

Download in other formats:

Plain Text