3 | | 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. |
4 | | |
5 | | 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. |
6 | | |
7 | | 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 |
8 | | [[LatexEquation( x_{i} \wedge i=1 \dots S' )]] [[BR]][[BR]] |
9 | | y llamar |
10 | | [[LatexEquation( s_{i} \wedge i=1 \dots S' )]] [[BR]][[BR]] |
11 | | 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 |
12 | | [[LatexEquation( S=\underset{i=1}{\overset{S'}{\sum}}s_{i} )]] [[BR]][[BR]] |
13 | | |
14 | | Sean los [[LatexEquation( k )]] puntos muestrales vecinos de [[LatexEquation( x_{i} )]] en orden de proximidad al mismo |
15 | | [[LatexEquation( y_{i,j} \wedge i=1 \dots S \wedge j=1 \dots k )]] [[BR]][[BR]] |
16 | | |
17 | | Sea la distancia euclídea del punto [[LatexEquation( x_{i} )]] a su [[LatexEquation( k )]]-ésimo vecino más próximo[[BR]][[BR]] |
18 | | [[LatexEquation( r_{i,k}=\sqrt{\underset{d=1}{\overset{n}{\sum}}\left(y_{i,k,d}-x_{i,d}\right)^{2}} )]] [[BR]][[BR]] |
19 | | |
20 | | 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]] |
21 | | [[LatexEquation( h_{i} = k+s_i)]] [[BR]][[BR]] |
22 | | cantidad que se distribuye como una binomial [[BR]][[BR]] |
23 | | [[LatexEquation( \eta_{i} = B\left(S, p_i\rigtht))]] |
24 | | 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]] |
25 | | [[LatexEquation( p_{i}=\int_{\left\Vert y-x_{i}\right\Vert ^{2}\leq r_{i,k}}\pi\left(y\right)\mathrm{d}y )]] [[BR]] |
26 | | 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 |
27 | | [[LatexEquation( \ln\pi_{i}=\ln\pi\left(x_{i}\right)+\lambda_0 )]] [[BR]][[BR]] |
28 | | [[LatexEquation( \ln\pi_{i,j}=\ln\pi\left(y_{i,j}\right)+\lambda_0 )]] [[BR]] |
29 | | |
30 | | |
31 | | 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 |
32 | | [[LatexEquation( r^n_{i,k} )]] [[BR]] |
33 | | obteniendo la relación |
34 | | [[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]] |
35 | | en la que [[LatexEquation(\lambda_1)]] es una constante desconocida. Puesto que la probabilidad no puede ser mayor que 1 tenemos una cota de la constante desconocida: [[BR]] [[BR]] |
36 | | [[LatexEquation( \lambda_{1}<-\underset{i=S'}{\max}\left\{ \ln\left(s_{i}\pi_{i}+\underset{j=1}{\overset{k}{\prod}}\pi_{i,j}\right)+n\ln r_{i,k}\right\} )]] [[BR]] |
37 | | |
38 | | También es posible mejorar la aproximación de la integral por interpolación, concretamente mediante el [http://www.alglib.net/interpolation/inversedistanceweighting.php método de Sheppard de ponderación inversa a la distancia] que es muy eficiente pues no requiere de ninguna evaluación extra. |
39 | | |
40 | | La probabilidad de que el número de puntos que caen dentro de la hiperesfera sea exactamente [[LatexEquation(h)]] será por tanto [[BR]][[BR]] |
41 | | [[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]] |
42 | | y el logaritmo de dicha probabilidad del contraste será |
43 | | [[LatexEquation( \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) )]] [[BR]] |
44 | | expresión en la cual se puede sustituir la anterior aproximación de [[LatexEquation( p_{i} )]]. |
45 | | |
46 | | La probabilidad de que el número de puntos que caen dentro de la hiperesfera sea mayor o igual que [[LatexEquation(h)]] se calcula mediante la función [http://www.gnu.org/software/gsl/manual/html_node/Incomplete-Beta-Function.html beta incompleta] [[BR]] [[BR]] |
47 | | [[LatexEquation( \mathrm{Pr}\left[\eta_{i}\leq h_{i}\right]=I_{1-p}\left(S-h_{i},h_{i+1}\right) )]] [[BR]] |
48 | | |