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 1 (modified by Víctor de Buen Remiro, 14 years ago) (diff)

--

Post-procesado de cadenas de Markov

Los métodos tradicionales de post-procesado basados en el burn-in y el thinning son demasiado arbitrarios para poder parametrizarlos de forma automática sin intervención del usuario.

Las cadenas simuladas con BysSampler cuentan con una ventaja adicional al conocerse la log-likelihood de cada muestra, pues esto permite contrastarla directamente con la densidad local empírica de los puntos cercanos 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 sobre-muestreadas y sustituirlos por puntos en otras zonas infra-muestreadas.

Una posibilidad sería utilizar el algoritmo KNN para encontrar los vecinos más próximos de cada punto de la muestra de tamaño S' < S . Como en los métodos de simulación tipo accept-reject suele haber bastantes 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}

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}}

Así las cosas tenemos que el número total de puntos muestrales en la hiperesfera de radio r_{i,k} y centro x_{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 esa hiperesfera

p_{i}=\int_{\left\Vert y-x_{i}\right\Vert ^{2}\leq r_{i,k}}\pi\left(y\right)\mathrm{d}y

Esa integral sería algo muy costoso de evaluar, pero lo que sí conocemos sin coste adicional es el logaritmo de esa densidad, 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

Podemos pues aproximar dicha integral como el producto de la media de las densidades por el hipervolumen de la región hiperesférica, que será proporcional a

r^n_{i,k}

obteniendo la relación

\ln p_{i}\approx\ln\left(\tilde{p}_{i}\left(\lambda_{1}\right)\right)\approx\lambda_{1}+\ln\left(\pi_{i}+\underset{j=1}{\overset{k}{\prod}}\pi_{i,j}\right)+n\ln r_{i,k}

en la que \lambda_1 es una constante desconocida. Puesto que la probabilidad no puede ser mayor que 1 tenemos una cota de ella:

\lambda_{1}\leq\lambda_{2}=-\underset{i=S'}{\max}\left\{ \ln\left(\pi_{i}+\underset{j=1}{\overset{k}{\prod}}\pi_{i,j}\right)+n\ln r_{i,k}\right\}

También es posible mejorar la aproximación de la integral 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. La probabilidad de que el número de puntos que caen dentro de la hiperesfera sea exactamente h será por tanto

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 del contraste 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)

expresión en la cual se puede aproximar como.

\ln\left(P_{i}\right)\approx\ln\left(\tilde{P}_{i}\left(\lambda_{1}\right)\right)\approx\ln\left(\begin{array}{c}S\\h_{i}\end{array}\right)+h_{i}\ln\tilde{p}_{i}\left(\lambda_{1}\right)-\left(S-h_{i}\right)\tilde{p}_{i}\left(\lambda_{1}\right)

Buscaremos el valor de \lambda_1 que maximiza la verosimilitud aproximada

\underset{\lambda_{1}\leq\lambda_{2}}{\max}\left\{ \underset{i=1}{\overset{S'}{\sum}}s_{i}\ln\left(\tilde{P}_{i}\left(\lambda_{1}\right)\right)\right\}

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}\left(S-h_{i},h_{i+1}\right)