id	summary	reporter	owner	description	type	status	priority	milestone	component	version	severity	resolution	keywords	cc
1049	[BysSampler] Post-procesado de cadenas de Markov	Víctor de Buen Remiro	Víctor de Buen Remiro	"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 [[LatexEquation( 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 [[LatexEquation( S' < S )]] puntos únicos 
  [[LatexEquation( x_{i} \wedge i=1 \dots S' )]] [[BR]][[BR]]
y llamar 
  [[LatexEquation( s_{i} \wedge i=1 \dots S' )]] [[BR]][[BR]]
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
  [[LatexEquation( S=\underset{i=1}{\overset{S'}{\sum}}s_{i} )]] [[BR]][[BR]]

Sean los [[LatexEquation( k )]] puntos muestrales vecinos de [[LatexEquation( x_{i} )]] en orden de proximidad al mismo
  [[LatexEquation( y_{i,j} \wedge i=1 \dots S \wedge j=1 \dots k )]] [[BR]][[BR]]

Sea la distancia euclídea del punto [[LatexEquation( x_{i} )]] a su [[LatexEquation( k )]]-ésimo vecino más próximo[[BR]][[BR]]
  [[LatexEquation( r_{i,k}=\sqrt{\underset{d=1}{\overset{n}{\sum}}\left(y_{i,k,d}-x_{i,d}\right)^{2}} )]] [[BR]][[BR]]

Así las cosas tenemos que el número total de puntos muestrales en la hiperesfera de radio [[LatexEquation( r_{i,k} )]] y centro [[LatexEquation( x_{i} )]] es [[BR]][[BR]]
  [[LatexEquation( h_{i} = k+s_i)]] [[BR]][[BR]]
cantidad que se distribuye como una binomial [[BR]][[BR]]
  [[LatexEquation( \eta_{i} = B\left(S, p_i\rigtht))]]
donde [[LatexEquation( p_i )]] es la probabilidad de la hiperesfera, es decir, la integral de la función de densidad en esa hiperesfera [[BR]][[BR]]
  [[LatexEquation( p_{i}=\int_{\left\Vert y-x_{i}\right\Vert ^{2}\leq r_{i,k}}\pi\left(y\right)\mathrm{d}y )]] [[BR]]
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 [[LatexEquation(\lambda_0)]] desconocida, evaluado en cada uno de los puntos muestrales, es decir, conocemos
  [[LatexEquation( \ln\pi_{i}=\ln\pi\left(x_{i}\right)+\lambda_0 )]] [[BR]][[BR]]
  [[LatexEquation( \ln\pi_{i,j}=\ln\pi\left(y_{i,j}\right)+\lambda_0 )]] [[BR]]

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
  [[LatexEquation( r^n_{i,k} )]] [[BR]]
obteniendo la relación
  [[LatexEquation( \ln p_{i}\approx\lambda_{1}+\ln\left(s_{i}\pi_{i}+\underset{j=1}{\overset{k}{\prod}}\pi_{i,j}\right)+n\ln r_{i,k} )]] [[BR]]
en la que [[LatexEquation(\lambda_1)]] es una constante desconocida.

La probabilidad de que el número de puntos que caen dentro de la hiperesfera sea exactamente [[LatexEquation(h)]] será por tanto [[BR]][[BR]]
  [[LatexEquation( 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}} )]] [[BR]]
y el logaritmo de dicha probabilidad del contraste será
  [[LatexEquation( \ln\left(P_{i}\right)=\ln\left(\begin{array}{c}S\\h_{i}\end{array}\right)+h_{i}p_{i}+\left(S-h_{i}\right)\ln\left(1-p_{i}\right) )]] [[BR]]
expresión en la cual se puede sustituir la anterior aproximación de [[LatexEquation( p_{i} )]]

Si la muestra fuera efectivamente perfectamente generada entonces sería lícito pensar que cierta transformación de Box-Cox de esta probabilidad debería tener una distribución normal independiente de cada punto[[BR]][[BR]]
  [[LatexEquation( T_{bc}\left(P_{i},\alpha,\beta\right)\sim N\left(\mu,\sigma^{2}\right) )]] [[BR]]
  [[LatexEquation( T_{bc}\left(P_{i},\alpha,\beta\right)=\begin{cases}\left(P_{i}+\alpha\right)^{\beta} & \forall\beta>0\\\ln\left(P_{i}+\alpha\right) & \forall\beta=0\end{cases} )]] [[BR]]

Es posible por lo tanto estimar los parámetros desconocidos  
  [[LatexEquation( \Theta=\left(\lambda_{1},\alpha,\beta,\mu,\sigma^{2}\right) )]]
que maximizan la verosimilitud de la normal así planteada.
  
Si suponemos que la cadena de partida no está perfectamente equilibrada pero sí ha recorrido el espacio convenientemente, entonces podemos buscar los outliers de las residuos en el modelo anterior, eliminando los puntos con residuos muy negativos pues corresponden a zonas superpobladas. 
"	task	new	highest	TOL Packages	Math		blocker