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 7 (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.

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 e incluso sustituirlos por puntos en otras zonas infra-muestreadas.

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

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

Distribución del cardinal local

Así las cosas tenemos que el cardinal 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, pero lo que sí 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$

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

$r^n_{i,k}$

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

De esta forma tenemos la fórmula de aproximación

$\ln p_{i}=\mu_0+\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}+\xi_{i}$

en la que el error se postulará por comodidad normal, homocedástico e independiente

$\xi_{i}\sim N\left(0,\sigma\right)$

aunque cabría pensar que en los entornos de relieve más complicado el error puede ser mayor que en las zonas más suaves. Dicho de otra forma

$\ln p_{i}\sim N\left(\mu_{i},\sigma\right)$

con

$\mu_{i}=\mu_0+\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}$

Verosimilitud de los parámetros

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$ y $\sigma$ será la suma ponderada por las repeticiones del producto de la probabilidad anterior por la densidad del error de aproximación, luego la expresión de su logaritmo será

$L\left(\mu_0,\sigma\right)=\underset{i=1}{\overset{S'}{\sum}}s_{i}\left(\ln\left(P_{i}\left(\mu_0\right)\right)-\frac{S'}{2}\ln\left(2\pi\right)-S'\ln\left(\sigma\right)-\frac{1}{2\sigma^{2}}\underset{i=1}{\overset{S'}{\sum}}\xi_{i}^{2}\right)$

Como esta expresión depende de los errores de aproximación que no son contrastables habría que simularlos, lo cual podría resultar muy costoso, o bien aproximar la esperanza

$E\left[L\left(\mu_0,\sigma\right)\right]=\underset{i=1}{\overset{S'}{\sum}}s_{i}\left(E\left[\ln\left(P_{i}\right)\right]-\frac{S'}{2}\ln\left(2\pi\right)-S'\ln\left(\sigma\right)-\frac{1}{2}S'\right)$

Para ello se puede desarrollar para la función

$f\left(\xi;\alpha,\beta,\gamma\right)=\alpha\left(\ln\left(\gamma\right)+\xi\right)+\beta\ln\left(1-\exp\left(\ln\left(\gamma\right)+\xi\right)\right)$

la serie de Taylor de orden 2 centrada en 0

$f\left(\xi;\alpha,\beta,\gamma\right)=f\left(0\right)+\left(\alpha-\frac{\beta\gamma}{1-\gamma}\right)\xi-\frac{1}{2}\beta\frac{\gamma}{1-\gamma}\left(1+\frac{\gamma}{1-\gamma}\right)\xi^{2}+O\left(\xi^{3}\right)$

Si $\xi\sim N\left(0,\sigma\right)$ entonces la esperanza de esta función es aproximadamente

$E\left[f\left(\xi;\alpha,\beta,\gamma\right)\right]\approx f\left(0\right)-\frac{1}{2}\frac{\gamma}{1-\gamma}\left(1+\frac{\gamma}{1-\gamma}\right)\beta\sigma^{2}$

Como resulta que

$\ln\left(P_{i}\right)=\ln\left(\begin{array}{c}S\\h_{i}\end{array}\right)+f\left(\xi_{i};h_{i},S-h_{i},e^{\mu_{i}}\right)$

su esperanza se puede aproximar como

$E\left[\ln\left(P_{i}\right)\right]\approx\ln\left(\begin{array}{c}S\\h_{i}\end{array}\right)+E\left[f\left(\xi_i;h_{i},S-h_{i},e^{\mu_{i}}\right)\right]$

es decir

$E\left[\ln\left(P_{i}\right)\right]\approx\ln\left(\begin{array}{c}S\\h_{i}\end{array}\right)+h_{i}\mu_{i}+\left(S-h_{i}\right)\ln\left(1-e^{\mu_{i}}\right)-\frac{1}{2}\frac{e^{\mu_{i}}}{1-e^{\mu_{i}}}\left(1+\frac{e^{\mu_{i}}}{1-e^{\mu_{i}}}\right)\left(S-h_{i}\right)\sigma^{2}$

Para que esta expresión sea evaluable debe cumplirse

$1-e^{\mu_{i}}>0\forall i=1\ldots S'$

Así pues tendremos el problema de optimización bivariante

$\underset{\mu_{0},\sigma}{\max}\left\{ E\left[L\left(\mu_{0},\sigma\right)\right]\right\}$

Sujeto a

$\mu_0\leq-\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\}$

$\sigma>0$

Test de super-població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}\left(S-h_{i},h_{i+1}\right)$

Download in other formats:

Plain Text