Inferencia estadística en procesos puntuales sobre grafos lineales
Abstract
Este trabajo constituye un recorrido que, desde el nivel del alumnado del Grado en Matemáticas, permite llegar a abordar problemas de inferencia no paramétrica sobre la función de intensidad de procesos puntuales definidos sobre grafos lineales. Para ello, hemos establecido una estructura en capítulos que ha de ser entendida como una consecución de peldaños de complejidad ascendente, en los que, de manera gradual, se van presentando los conceptos y herramientas necesarias hasta llegar al desarrollo del problema final (y objetivo de este trabajo). Comenzamos introduciendo una serie de conceptos esenciales de la Estadística y la Teoría de la Probabilidad que se necesitan y manejan a lo largo de todo el trabajo. En el Capítulo 1 nos centramos en un problema clásico y bien conocido que es la estimación de la función de densidad, focalizándonos en las técnicas no paramétricas, y en particular en los métodos núcleo. Este capítulo no ha de entenderse una mera introducción, ya que es de esencial importancia en el estudio de los procesos puntuales por la íntima relación existente entre las funciones de densidad y las funciones de intensidad. En el Capítulo 2, presentamos los procesos puntuales en el plano euclídeo y sus modelos esenciales, como los procesos de Poisson. Una vez que se manejan estos conceptos sobre el plano euclídeo, en el Capítulo 3 se incrementa un grado el nivel de complejidad presentándolos sobre grafos lineales, en donde ya no disfrutamos de ventajas como la existencia de la métrica euclídea o del concepto de gradiente. De nuevo, focalizamos nuestro estudio en los procesos de Poisson, estudiando en profundidad diversas técnicas de estimación no paramétrica de la función de intensidad, así como la cuestión clave de los selectores del parámetro ventana. Uno de los problemas más interesantes y estudiados en la literatura estadística es el de comparación de dos (o más) poblaciones. El Capítulo 4 se dedica íntegramente a la presentación de este problema en el marco de los procesos puntuales en grafos lineales. Además, se aportan soluciones innovadoras que consisten en tres test estadísticos que permiten concluir si, dados dos patrones de puntos sobre un mismo grafo lineal, provienen de procesos con funciones de intensidad proporcionales, es decir, con una densidad espacial común que se traduce en una misma estructura espacial (en términos de propiedades de primer orden). En el Capítulo 5 se lleva a cabo un exhaustivo estudio de simulación con varios escenarios en los que se comprueba la calidad y buen comportamiento de los métodos de contraste propuestos en el capítulo anterior. Para ello se presentan tanto resultados sobre el ajuste de cada contraste en distintos tipos de grafos lineales, así como de la potencia de los mismos. Para finalizar el trabajo, se presenta una aplicación de parte de los métodos descritos para el contraste de dos poblaciones sobre un conjunto de datos reales de accidentes de tráfico en la ciudad de Río de Janeiro (Brasil) entre 2019 y 2022. Además, en el apéndice de este trabajo puede encontrarse el código empleado, con garantías de reproducibilidad. Cabe decir que no se incluye la base de datos reales por cuestiones de confidencialidad
Full text
Traballo Fin de Grao Inferencia estadística en procesos puntuales sobre grafos lineales Ignacio González Pérez Julio, 2022 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA
GRAO DE MATEMÁTICAS Traballo Fin de Grao Inferencia estadística en procesos puntuales sobre grafos lineales Ignacio González Pérez Julio, 2022 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA
Trabajo propuesto Área de Coñecemento: Estatística e Investigación Operativa Título: Inferencia Estadística en procesos puntuales sobre grafos lineales Breve descrición do contido Este trabajo tiene por objetivo general realizar tareas de Inferencia Estadística en un contexto particular en el que los elementos de interés son eventos o sucesos que ocurren sobre un grafo lineal. La motivación del estudio de este tipo de elementos surge del hecho de que muchos eventos geolocalizados, como pueden ser los accidentes de trácos, colisiones con animales salvajes, eventos delictivos... tienen lugar sobre un soporte que puede ser caracterizado como un grafo lineal. Trabajar con este soporte lineal embebido en el espacio euclídeo bidimensional supone nuevos retos, asociados especialmente al uso de la métrica del camino más corto que sustituye a la habitual distancia euclídea. Este trabajo se estructurará como sigue: 1) Revisión de herramientas de Inferencia Estadística para espacios bidimensionales. 2) Introducción a los procesos puntuales sobre el plano euclídeo. 3) Introducción a la teoría de grafos lineales. 4) Estimación de la función de densidad en procesos puntuales sobre grafos lineales. 5) Ilustración con simulaciones o bases de datos reales. Recomendacións Outras observacións iii
Índice Resumen viii Introducción xiii 1. Estimación de la función de densidad en espacios euclídeos 1 1.1. Estimación no paramétrica de la densidad unidimensional . . . . . . . . . . . . . 1 1.2. Estimación no paramétrica de la densidad multidimensional . . . . . . . . . . . . 7 2. Procesos puntuales en el plano euclídeo 11 2.1. Procesospuntuales................................... 11 2.2. Funcióndeintensidad ................................. 14 3. Procesos puntuales en grafos lineales 19 3.1. Procesos puntuales en grafos lineales . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.2. Funcióndeintensidad ................................. 23 4. Comparación de funciones de intensidad 31 4.1. Test de Kolmogorov-Smirnov . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 4.2. TestdeCramervonMises............................... 34 4.3. Test no paramétrico basado en la función de riesgo relativo . . . . . . . . . . . . 35 4.4. Procedimiento de calibración . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 v
vi Índice 5. Estudio de simulación 41 5.1. Modelos bajo la hipótesis nula . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 5.2. Modelos bajo la hipótesis alternativa . . . . . . . . . . . . . . . . . . . . . . . . . 43 5.3. Resultados........................................ 45 6. Aplicación a datos reales 53 A. Código 59 A.1. Cálculo de los estadísticos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 A.2. Estimación de niveles críticos empleando el test de permutaciones . . . . . . . . . 69 A.3.Estudiodesimulación ................................. 74 Índice de notación 81 Bibliografía 85
xiv Introducción Denición 0.4. Sea X= (X1, . . . , Xm)0 un vector aleatorio m -dimensional. Este se dice absolutamente continuo si existe una función f:Rm→R , que llamaremos su función de densidad de probabilidad, tal que: P(X∈A) = ZA f(x)dx, para cualquier A⊂Rm perteneciente a la σ -álgebra relativa al espacio de probabilidad sobre el que está denido X . De las propiedades de la probabilidad se deduce que necesariamente f es una función no negativa y que RRmf(x)dx= 1 . Además, cualquier función que verique estas dos propiedades determina de forma unívoca una distribución de probabilidad en Rm ; es decir, las dos condiciones anteriores caracterizan las funciones de densidad de los vectores aleatorios absolutamente continuos. En este mismo contexto podemos denir la función de distribución de X y relacionarla con su función de densidad, como: Denición 0.5. Sea X un vector aleatorio m -dimensional. Denimos su función de distribución (conjunta) como la aplicación F:Rm→R tal que: F(x) = F(x1, . . . , xm) = P(X1≤x1, . . . , Xm≤xm) = Zxm −∞ . . . Zx1 −∞ f(y1, . . . , ym)dy1···dym, ∀x= (x1, . . . , xm)0∈Rm , lo que denotaremos usualmente por F(x) = P(X≤x) .
Capítulo 1 Estimación de la función de densidad en espacios euclídeos A la hora de estimar la función de densidad de un vector aleatorio X , el escenario usual suele ser contar con una muestra aleatoria simple de este vector; es decir, un conjunto de vectores aleatorios X1,...,Xn independientes entre sí e idénticamente distribuidos a X cuya realización nos proporcionará una muestra. En aquellos casos en los que se pueda suponer que la distribución de X pertenece a una familia paramétrica conocida, los parámetros desconocidos pueden estimarse a partir de la muestra empleando técnicas de momentos, máxima verosimilitud, etc. Este es el contexto de la estimación paramétrica de la función de densidad. Si carecemos de información acerca de la distribución de X , lo más natural es buscar un método que permita estimar la función de densidad dejando a los datos hablar por sí mismos; es decir, sin imponer ningún tipo de restricción a la forma de f . Entramos así en el contexto de la estimación no paramétrica de la función de densidad. 1.1. Estimación no paramétrica de la densidad unidimensional Supongamos inicialmente que nos encontramos en un caso unidimensional. La estimación más elemental de la densidad es la frecuencia relativa, que toma la forma de un histograma. Dividiendo la recta real en intervalos de igual longitud h 1 , que denominamos ventanas o bandas, la estimación de la función de densidad dada por el histograma en un punto x∈R es: b f(x;h) = n º de observaciones en la ventana que contiene a x nh . 1
2 1. Estimación de la función de densidad en espacios euclídeos (a) h= 2, x0=−3 . (b) h= 0.5, x0=−3 . (c) h= 1, x0=−3 . (d) h= 1, x0=−2.5 . Figura 1.1: distintos histogramas para una muestra de tamaño 100 de la distribución N(2,2.25) , donde h es el ancho de banda y x0 el extremo izquierdo de la primera ventana. La densidad real se representa de forma punteada. A la hora de construir un estimador como este, han de tomarse dos decisiones: el ancho de banda, h , y el punto inicial x0 . Ambas decisiones son capitales de cara a la forma nal del estimador. En la Figura 1.1 vemos como diferentes elecciones del ancho de banda o del punto inicial generan estimadores considerablemente distintos para la misma muestra. Este problema se soluciona, en primera instancia, introduciendo el concepto de histograma móvil. Fijado un ancho de banda h , para cada i∈ {1, . . . , n} se considera el histograma que tiene como única observación Xi y toma una ventana centrada en este punto. El estimador b f se obtiene promediando todos estos histogramas. De esta forma, podemos escribir: b f(x;h) = 1 nh n X i=1 1[−1/2,1/2] x−Xi h. Si notamos que 1[−1/2,1/2] es la función de densidad de una distribución uniforme en [−1/2,1/2] , podemos plantearnos sustituirla por otras funciones de densidad k que varíen más suavemente con la distancia, suavizando nuestro estimador. Surgen así los estimadores tipo núcleo de la función de densidad: b f(x;h) = 1 nh n X i=1 kx−Xi h=1 n n X i=1 kh(x−Xi), (1.1) donde hemos introducido la notación kh(·) = h−1k(·/h) . La función núcleo k ha de ser una densidad unimodal simétrica respecto al cero. Estas condiciones aseguran que el estimador b f sigue siendo una función de densidad. Uno de los criterios de referencia para determinar la calidad de un estimador es su error cuadrático medio, denotado por sus siglas en inglés, Mean Square Error : 1 En un contexto más general, estos intervalos no tienen que tener la misma longitud.
1.1. Estimación no paramétrica de la densidad unidimensional 3 Denición 1.1. Sea θ una cantidad desconocida y ˆ θ un estimador de θ . Se dene el error cuadrático medio de este estimador como: MSE ˆ θ=Eˆ θ−θ2=V ar ˆ θ+hEˆ θ−θi2 . Si f es la verdadera densidad de X , tenemos que, bajo muestreo aleatorio simple: Ehb f(x;h)i=1 n n X i=1 E[kh(x−Xi)] = E[kh(x−X)] = Zkh(x−y)f(y)dy =Zk(z)f(x−hz)dz. Antes de entrar en más detalle debemos explicitar las hipótesis con las que vamos a trabajar. Bajo estas podremos hacer cálculos asintóticos del MSE del estimador tipo núcleo, así como asegurar el buen comportamiento de este. Las hipótesis son: (Ui) El ancho de banda h≡hn se representa por una sucesión no aleatoria de valores que verican, omitiendo el subíndice n , que l´ım n→∞ h= 0 pero l´ım n→∞ nh =∞ . (Uii) El núcleo k es una densidad de probabilidad acotada, simétrica entorno al origen, y con momento µl(k) = Rxlk(x)dx nito para l= 0,1,2 . (Uiii) La densidad f es tal que f00 es continua y de cuadrado integrable. Antes de realizar estos cálculos asintóticos debemos introducir el concepto de o pequeña. Dadas (an)n∈N,(bn)n∈N dos sucesiones de números reales, diremos que (an)n∈N es una o pequeña de (bn)n∈N , y escribiremos an=o(bn) , si l´ım n→∞ |an/bn|= 0 . Notemos que bajo la hipótesis (Ui) cualquier función de h puede entenderse como una sucesión de números reales. La hipótesis (Uiii) permite hacer un desarrollo de Taylor de f(x−hz) centrado en x hasta orden 2, como puede consultarse en [10]. De esta forma: Ehb f(x;h)i=f(x)Zk(z)dz −hf0(x)Zzk(z)dz +1 2h2f00(x)Zz2k(z)dz +o(h2) =f(x) + 1 2h2f00(x)µ2(k) + o(h2), (1.2) en donde debemos notar que, bajo las hipótesis previas, µ1(k)=0 . Además, la hipótesis (Ui) garantiza la insesgadez asintótica del estimador; en efecto, ya que l´ım n→∞ Ehb f(x;h)i=f(x) . Analicemos ahora el término asociado a la varianza. Empleando la hipótesis de que tenemos muestreo aleatorio simple: V ar hb f(x;h)i=V ar "1 n n X i=1 kh(x−Xi)#=1 nV ar [kh(x−X)] =1 nZ[kh(x−y)−E[kh(x−X)]]2f(y)dy =1 nZk2 h(x−y)f(y)dy −2E[kh(x−X)] Zkh(x−y)f(y)dy
4 1. Estimación de la función de densidad en espacios euclídeos +E[kh(x−X)]2Zf(y)dy =1 nZk2 h(x−y)f(y)dy −1 nE[kh(x−X)]2 =1 nh Zk2(z)f(x−hz)dz −1 nE[b f(x;h)]2. Haciendo ahora un desarrollo de Taylor de f centrado en x , pero ahora a orden cero, y empleando el resultado de la ecuación (1.2), vemos que: V ar hb f(x;h)i=1 nh Zk2(z) [f(x) + o(1)] dz −1 n[f(x) + o(h)]2 =1 nhf(x)Zk2(z)dz +o[(nh)−1]−o(n−1) = (nh)−1R(k)f(x) + o[(nh)−1], (1.3) donde hemos denido R(g) = Rg2(x)dx . Debemos notar que la hipótesis (Ui) garantiza que la varianza del estimador converge a cero, lo que nos permite concluir que el estimador tipo núcleo es asintóticamente consistente. Teniendo entonces en cuenta los resultados obtenidos en las ecuaciones (1.2) y (1.3), concluimos que el error cuadrático medio del estimador tipo núcleo de f en x es: MSE hb f(x;h)i= (nh)−1R(k)f(x) + h4 4µ2(k)2f00(x)2+o(nh)−1+h4. (1.4) El principal problema que presenta la expresión obtenida en la ecuación (1.4) es que se trata de un indicador local de la calidad de la estimación. Para obtener un indicador global, integramos esta expresión para obtener el denominado error cuadrático medio integrado MISE: MISE hb f(·;h)i=Z MSE hb f(x;h)idx = (nh)−1R(k) + h4 4µ2(k)2R(f00) + o(nh)−1+h4. De quedarnos con los términos dominantes, obtenemos una estimación asintótica global de la calidad de nuestro estimador no paramétrico, el denominado error cuadrático medio integrado asintótico AMISE: AMISE hb f(·;h)i= (nh)−1R(k) + h4 4µ2(k)2R(f00), (1.5) donde el primer término es relativo a la varianza, mientras que el segundo es el relativo al sesgo cuadrático. Notemos que, dada la muestra, según h→0 el sesgo cuadrático decrece mientras que el término asociado a la varianza aumenta; es decir, según tomamos anchos de banda más pequeños vamos obteniendo estimadores con menor sesgo pero con mayor varianza. Esto es conoce como compensación varianza-sesgo, o en la literatura anglosajona como variance-bias trade-o . Este comportamiento se ejemplica en la Figura 1.2; en (a) tenemos un ancho de banda demasiado pequeño que genera una estimación con excesiva varianza pero poco sesgo, lo que nos
1.1. Estimación no paramétrica de la densidad unidimensional 5 (a) h= 0.2 . (b) h= 0.5 . (c) h= 1 . Figura 1.2: estimaciones tipo núcleo para una muestra empleada en la Figura 1.1. Se ha empleado en todas en núcleo de Epanechnikov con distintos anchos de banda. La densidad real se representa de forma punteada. lleva a decir que nuestro estimador está infrasuavizado. Según aumentamos h , vemos como el estimador se suaviza buscando un equilibrio sesgo-varianza tal como ocurre en (b). Finalmente, anchos de banda demasiado grandes generan estimadores con excesivo sesgo, denominados sobresuavizados, como el que vemos en (c). Esta compensación del error de estimación entre la varianza y el sesgo cuadrático motiva buscar el ancho de banda óptimo. Una opción es tomar h tal que se minimice el AMISE calculado en la ecuación (1.5). Derivando respecto a h , se concluye que el ancho de banda óptimo para este criterio es: h AMISE =R(k) nµ2(k)2R(f00)1/5 , (1.6) que si lo sustituimos de vuelta en la ecuación (1.5), llegamos a que: ínf h>0 AMISE hb f(·;h)i=5 4µ2(k)2R(k)4R(f00)1/5n−4/5. (1.7) De la ecuación (1.7) deducimos que el orden de consistencia (asintótico) de nuestro estimador no paramétrico de la función de densidad es n−4/5 , que es menor que la tasa usual de convergencia paramétrica de n−1 . Esto nos indica que, aunque mucho más versátiles, los estimadores tipo núcleo no convergen con tanta rapidez como los paramétricos. A la hora de elegir el ancho de banda con el que construir un estimador tipo núcleo, lo lógico sería emplear el obtenido en la ecuación (1.6). El problema reside en que este depende de cantidades teóricas que son desconocidas, pues hacen referencia a la densidad de X . Por ello, debemos determinar métodos que permitan seleccionar un h que proporcione una estimación de calidad. Una primera idea, propuesta en [27], es asumir la hipótesis de que nuestra muestra procede de una población normal. En este caso podemos calcular R(f00) = 3/(8√πσ5) , lo que permite
6 1. Estimación de la función de densidad en espacios euclídeos estimar el ancho de banda óptimo para el AMISE como: h Norm =8√πR(k) 3nµ2(k)21/5 σ. Ahora bien, en muchos de los casos ni siquiera se tienen indicios de normalidad; por ello, se han diseñado métodos más generales que permitan obtener un ancho de banda efectivo en ausencia de más hipótesis sobre la población que las que tenemos hasta el momento. Primeramente, notamos que el error cuadrático medio integrado puede descomponerse como: MISE hb f(·;h))i=EZb f(x;h)2dx−2EZb f(x;h)f(x)dx+Zf(x)2dx, y como Rf2(x)dx no depende de h , minimizar el MISE equivale a minimizar: MISE hb f(·;h)i−Zf(x)2dx =EZb f(x;h)2dx −2Zb f(x;h)f(x)dx. (1.8) A priori puede parecer que tenemos el mismo problema que con el ancho de banda óptimo en la ecuación (1.6), ya que el segundo sumando de la función objetivo anterior depende de f , que es desconocida. Ahora bien, si notamos que Rb f(x;h)f(x)dx =Ehb f(X, h)i , podemos tratar de estimar esta cantidad mediante n−1Pn i=1 b f(xi, h) . El problema con este estimador, como puede consultarse en [27], es que presenta un fuerte sesgo negativo al evaluar b f en puntos de la muestra. Para solucionar este problema, Bowman propone en [7] sustituir b f(xi, h) por la estimación no paramétrica de la densidad en xi empleando todos los datos de la muestra a excepción del propio xi , que denotaremos por b f−i(xi, h) . Surge así la función del ancho de banda: LSCV (h) = Zb f(x;h)2dx −2 n n X i=1 b f−i(Xi;h), que es un estimador insesgado de la función objetivo dada en la ecuación (1.8). En efecto, tenemos primeramente que: E"1 n n X i=1 b f−i(Xi;h)#=1 n n X i=1 Ehb f−i(Xi;h)i=1 n(n−1) n X i=1 n X j=1 j6=i E[kh(Xi−Xj)] . Dadas X y Z variables aleatorias independientes e idénticamente distribuidas con densidad f , la función de densidad de Y=X−Z es g(y) = Rf(x−y)f(x)dx . Así, bajo muestreo aleatorio simple: E"1 n n X i=1 b f−i(Xi;h)#=1 n(n−1) n X i=1 n X j=1 j6=iZkh(y)g(y)dy =ZZ kh(y)f(x)f(x−y)dxdy.
1.2. Estimación no paramétrica de la densidad multidimensional 7 Veamos que la esperanza del segundo sumando de la ecuación (1.8) coincide con este resultado. En efecto: EZb f(x;h)f(x)dx =Zf(x)Ehb f(x;h)idx =Zf(x)Zkh(x−y)f(y)dydx =ZZ kh(x−y)f(x)f(y)dxdy =ZZ kh(y)f(x)f(x−y)dxdy, tal y como queríamos ver. Así, LSCV es un estimador insesgado de la función objetivo dada en la ecuación (1.8). Por ello, eligiendo el ancho de banda h que minimice LSCV( h ) (que notemos es una función que podemos computar para cada h ) obtenemos un ancho de banda que, presumiblemente, conducirá a una estimación de f con un MISE cercano a su valor mínimo. El empleo de este tipo de estimadores, en los que se estiman cantidades evaluadas en una observación empleando todas las demás, recibe el nombre de validación cruzada (en inglés cross-validation). Es por ello que esta técnica para calcular un h óptimo se conozca como Least Squares Cross-Validation , ya que se emplean técnicas de validación cruzada para minimizar el error cuadrático. Técnicas más elaboradas de selección del ancho de banda construyen procesos en varias etapas, tras cada una de las cuales se obtiene un ancho de banda que proporciona una estimación de mejor calidad. Otros métodos consisten también en elegir otros funcionales que estimen el MISE, y tratar de minimizarlos. Estas técnicas de selección de h , junto con muchas otras, pueden consultarse en [27]. 1.2. Estimación no paramétrica de la densidad multidimensional En un escenario multivariante, X es ahora un vector aleatorio d -dimensional del cual poseemos una m.a.s. X1,...,Xn . A la hora de extender nuestro estimador no paramétrico de la función de densidad al caso multivariante, la principal diferencia se encuentra en que el ancho de banda h pasa a ser una matriz de ancho de banda H . Tomaremos H∈ F ={A∈ Md×d(R): A simétrica y denida positiva } . Esto es natural, ya que si recordamos que en el caso univariante h medía la dispersión asociada a los núcleos promediados, en el caso multivariante jugará el papel de una matriz de covarianzas. De esta forma, el estimador tipo núcleo d -dimensional de la función de densidad, como puede consultarse en [27], es: b fF(x;H) = 1 n n X i=1 KH(x−Xi), donde KH(x) = |H|−1/2K(H−1/2x), siendo la función núcleo K una densidad de probabilidad en Rd . Al ser H una matriz d×d simétrica, posee d(d+ 1)/2 entradas independientes. Con el objetivo de simplicar el problema de selección de esta matriz (en el cual profundizaremos más adelante) es costumbre imponer restricciones sobre la forma de H , como por ejemplo H∈ D ={ diag (h2 1, . . . , h2 d)∈ Md×d(R): h∈Rd} ,
8 1. Estimación de la función de densidad en espacios euclídeos (a) Diagrama de dispersión. (b) H= diag (1.72,192) . (c) H= diag (32,302) . (d) H= diag (1,142) . Figura 1.3: diagrama de dispersión (a) y curvas de nivel de estimaciones de la densidad bivariante. Se ha empleado un núcleo normal bivariante alineado con los ejes y varias matrices H∈ D . Datos extraídos de [26]. en cuyo caso el estimador pasa a tener la forma: b fD(x;h) = 1 n d Y j=1 hj −1n X i=1 Kx1−Xi1 h1 ,...,xd−Xid hd. Un ejemplo de estimación en el caso bivariante empleando una matriz H de esta clase puede verse en la Figura 1.3, donde podemos apreciar nuevamente el compensación varianza-sesgo según varían los valores de hj . Una mayor simplicación es tomar H∈ S ={h2Id:h > 0} . Este último caso presenta especial interés tanto teórico como práctico debido a la simplicación del estimador, el cual pasa a ser: b fS(x;h) = 1 nhd n X i=1 Kx−Xi h. (1.9) Para estudiar las propiedades de estos estimadores, debemos calcular el AMISE. Por simplicidad, nos restringiremos a estimadores de la forma dados en la ecuación (1.9), debido a que presentan la mejor relación coste-benecio tanto conceptual como computacional. Para poder obtener estimaciones asintóticas del error cuadrático de los estimadores anteriores, emplearemos la siguiente versión débil del teorema de Taylor en varias variables:
1.2. Estimación no paramétrica de la densidad multidimensional 9 Teorema 1.2. Sea g∈ C2(Rm,R) , y (αn)n∈N una sucesión de vectores de Rm convergente a 0 . Denotemos por ∇g(x) al vector gradiente de g en x , y por Hg(x) a la matriz Hessiana de g en x ; es decir, la matriz m×m tal que (Hg(x))ij =Dijg(x) . Además, dada A∈ Mm×m(R) denimos su traza como tr(A) = Pm i=1 Aii . Tenemos entonces que: g(x+αn) = g(x) + α0 n∇g(x) + 1 2tr Hg(x)αnα0 n+o(α0 nαn). Versiones más generales del teorema de Taylor en varias variables pueden consultarse en [24]. Además del resultado anterior, para obtener dichas estimaciones, así como para poder asegurar el buen comportamiento del estimador dado en la ecuación (1.9), asumiremos las siguientes hipótesis establecidas en [27]: (Mi) Las entradas de la matriz hessiana de f , Hf son continuas y de cuadrado integrable. (Mii) El ancho de banda h se representa por una sucesión no aleatoria de valores que verican, omitiendo el subíndice n , que l´ım n→∞ h= 0 pero l´ım n→∞ nhd=∞ . (Miii) La función núcleo K posee un soporte compacto, está acotada, y es tal que: RK(z)dz= 1 , RzK(z)dz=0 , y Rzz0K(z)dz=µ2(K)Id ; siendo µ2(K) = Rz2 iK(z)dz invariante con i . Para la esperanza tenemos entonces que: Ehb fS(x;h)i=1 hdEKx−X h=1 hdZKx−y hf(y)dy=ZK(z)f(x−hz)dz =ZK(z)f(x)−hz0∇f(x) + h2 2tr Hf(x)zz0+o(h2)dz =f(x)−h∇f(x)0ZzK(z)dz+h2 2ZK(z)tr Hf(x)zz0dz+o(h2) =f(x) + h2 2tr Hf(x)Zzz0K(z)dz+o(h2) =f(x) + h2 2µ2(K)tr [Hf(x)] + o(h2). El resultado anterior bajo la hipótesis (Mii) nos permite concluir, al igual que ocurría univariantemente, que el estimador dado en la ecuación (1.9) es asintóticamente insesgado. La varianza de este estimador es: V ar hb fS(x;h)i=1 nh2dV ar Kx−X h=1 nh2d"E"Kx−X h2#−EKx−X h2# =1 nh2dZKx−y h2 f(y)dy−1 nEhb f(x;h)i2 =1 nhdZK(z)2f(x−hz)dz−1 n[f(x) + o(h)]2
16 2. Procesos puntuales en el plano euclídeo (a) Patrón de puntos. (b) Estimación mediante b λ0 . (c) Estimación mediante b λD . Figura 2.2: localización de pinos negros Japoneses en una región de muestreo de un bosque (a), junto con dos estimaciones no paramétricas de su intensidad (b), (c). En ambas se ha empleado un núcleo Gaussiano estándar y el ancho de banda se ha tomado empleando la regla de Scott [25] isótropa. Datos extraídos del paquete [4]. Del estimador tipo núcleo dado en la ecuación (2.3) debemos destacar que se ha tomado como matriz de ancho de banda H=h2I2 , y que el núcleo L se entiende una función de densidad de probabilidad bivariante, isótropa 1 , unimodal y centrada en 0 . Notemos además que en el estimador b λ0 no dividimos entre N , a diferencia de en los estimadores estudiados en el Capítulo 1, ya que se está estimando una función de intensidad a través de su densidad asociada. El principal problema que presenta el estimador dado en la ecuación (2.3) es que decrece fuertemente cerca de la frontera de W . Si, por ejemplo, nuestro proceso puntual es la ubicación de árboles en una región de muestreo en un bosque, como en el ejemplo de la Figura 2.1, los árboles que se encuentran cerca de los límites de W , pero fuera de esta región, no se están considerando a la hora de estimar la intensidad. Esto se conoce como un efecto frontera. Este tipo de efectos causan que b λ0 presente un fuerte sesgo negativo cerca de la frontera de W . Tratando de corregir estos problemas se propone la corrección de Diggle [12]: b λD(y) = 1 h2 N X i=1 1 e(xi)Ly−xi h, donde e(z) = ZW L(z−w)dw. (2.4) Esta corrección trata de compensar los efectos de borde al ponderar por e(xi)−1 , que es una medida de cómo de cerca está xi de la frontera de W . Esto puede verse en la Figura 2.2, donde presentamos el patrón de puntos dado en la Figura 2.1 con dos estimaciones de la función de intensidad; en (b) se ha empleado el estimador dado en (2.3), mientras que en (c) se ha empleado 1 Tal que la distribución que representa posee por matriz de covarianzas Σ=σ2I2 .
2.2. Función de intensidad 17 la corrección de Diggle. Podemos ver en estas guras como b λ0 toma valores mucho menores que b λD cerca de la frontera de W , consecuencia de los efectos frontera. En determinados contextos en los que se estudia un proceso puntual se conocen los valores de una magnitud a lo largo de la región de observación. Estos observables se denominan covariables del proceso puntual, y acostumbran a describirse por una función Z:W→R . Idealmente, los valores de una covariable se conocen de forma precisa a lo largo de todo W ; en la práctica, suele conocerse la imagen de Z a lo largo de una malla de W sucientemente na. Si estudiando un proceso puntual en presencia de una covariable tenemos evidencias de inhomogeneidad, resulta natural pensar que la covariable pueda tener un efecto sobre la intensidad de dicho proceso puntual. Así, de forma alternativa a la estimación tipo núcleo, podemos tratar de buscar una función ρ tal que: λ(y) = ρ[Z(y)] . Como generalmente carecemos de ninguna información a cerca de la forma funcional de ρ , optamos por una estimación no paramétrica como la propuesta en [3]. La idea fundamental consiste en determinar la distancia entre dos puntos del plano no empleando la distancia euclidiana, sino la función que describe la covariable Z ; concretamente, empleando el área de la región comprendida entre dos curvas de nivel de esta función. Para ello, comenzamos deniendo la función de acumulación de Z como: C(z) = 1 |W|ZW 1{Z(w)<z}dw. Usualmente nos veremos en la necesidad de estimar C , al no conocer la forma precisa de Z a lo largo de todo W . Para su estimación, dividimos la región de observación en una na malla a lo largo de la cual conocemos Z , y estimamos la función de acumulación como: b C(z) = #{ puntos y de la malla en los que Z(y)≤z} #{ puntos de la malla }, donde #A denota el cardinal de un conjunto A . Una vez conocida (o estimada) C , estimamos ρ a través de un estimador tipo núcleo transformado: bρ(z) = 1 |W| N X i=1 kh[C(z)−C[Z(xi)]], lo que proporciona la estimación de la intensidad: b λ(y) = 1 |W| N X i=1 kh[C[Z(y)] −C[Z(xi)]], (2.5) donde k en este caso ha de ser una función núcleo apropiada para la estimación no paramétrica en el caso unidimensional, y kh(·) = h−1k(·/h) . Notemos que: C[Z(y)] −C[Z(xi)] ∝ZW 1{Z(w) entre Z(y) y Z(xi)}dw,
18 2. Procesos puntuales en el plano euclídeo que es el área encerrada entre las curvas de nivel de Z que pasan por y y xi . De esta forma, vemos como el estimador de la función de intensidad dado en la ecuación (2.5) no es más que un estimador tipo núcleo dónde la distancia entre dos puntos se mide a través del área de la región delimitada por las curvas de nivel de Z que pasan por dichos puntos, normalizada por el área de la región de observación. 2.2.2. Marcas Los puntos de una realización de un proceso puntual pueden llevar asociados múltiples atributos. Cualquier información asociada a cada punto de dicha realización se denomina marca. Esto nos lleva a hablar de patrones de puntos marcados: y ={(x1,m1),...,(xN,mN)}, donde x ={xi}N i=1 es un patrón de puntos en el plano y {mi}N i=1 es el conjunto de marcas de x . Esto lleva a denir un proceso puntual marcado como aquel mecanismo aleatorio cuya realización es un patrón de puntos marcado. Las marcas de un patrón de puntos pueden ser de múltiples tipos: desde valores categóricos hasta observaciones multivariantes. La gran diferencia entre una marca y una covariable de un proceso puntual es que las marcas son intrínsecas a este y sus valores se obtienen como parte la realización del proceso puntual, mientras que las covariables son extrínsecas; es decir, sus valores a lo largo de la región de observación no son aleatorios y no dependen de la realización del proceso puntual.
Capítulo 3 Procesos puntuales en grafos lineales En el capítulo anterior hemos estudiado los procesos puntuales en el plano euclídeo. Las realizaciones de estos procesos puntuales eran patrones de puntos en una región de observación W⊂R2 , que considerábamos conexa y compacta en la topología usual. Ahora bien, en determinadas ocasiones nos encontramos ante patrones de puntos que no se disponen a lo largo de un subconjunto 2-dimensional del plano, sino a lo largo de un conjunto 1-dimensional embebido en este: accidentes de tráco en una red de carreteras, defectos en una instalación eléctrica, escapes a lo largo de un sistema de cañerías, avistamientos de aves a lo largo de la línea de costa... Este tipo de eventos están constreñidos a ocurrir en espacios de una dimensión, lo que supone una diferencia sustancial con los patrones de puntos que hemos estudiado hasta el momento. Debemos entonces denir con precisión los espacios en los que vamos a estudiar este tipo de patrones de puntos. Introducimos para ello el concepto de grafo lineal [8]: Denición 3.1. Un grafo lineal G es una dupla G= (V,A) donde V={v1,...,vnv} ⊂ R2 es el conjunto de nodos o vértices del grafo, y A={l1, . . . , lnl} es el conjunto de aristas del grafo. Estas aristas son segmentos que empiezan y acaban en nodos del grafo: ∀i∈ {1, . . . , nl} ∃ vi1,vi2∈ V, tales que vi16=vi2 y li={(1 −t)vi1+tvi2∈R2:t∈[0,1]}. Asumiremos que dos aristas únicamente pueden intersecarse en nodos del grafo. Denimos además el subconjunto del plano representado por este grafo como la unión de sus aristas: LG=Snl i=1 li⊂R2 , que supondremos conexo. Debemos comentar que existe cierta dualidad en la representación de un grafo lineal. Primeramente, es frecuente abusar del lenguaje y denominar grafo lineal, o simplemente grafo, tanto a la dupla G como al subconjunto del plano LG . Intuitivamente, podría parecer que la descripción apropiada de un grafo sería LG , pues es lo que representamos grácamente y sobre la que se puede 19
20 3. Procesos puntuales en grafos lineales efectuar cálculo diferencial e integral. Indistintamente, la representación en dupla nodos-aristas de un grafo lineal es apropiada tanto desde un punto de vista algebraico como computacional. Empleando grafos lineales pueden representarse una amplia variedad de espacios. Retomando los ejemplos previos: en un mapa las calles o carreteras se representan por aristas y los cruces por nodos; en una instalación eléctrica los cables se identican con las aristas y las clemas con los nodos; el litoral de una determinada región puede aproximarse por segmentos consecutivos, unidos pos nodos. Esta versatilidad hace de los grafos lineales la herramienta apropiada para describir los espacios 1-dimensionales contenidos en el plano sobre los que se observan los patrones de puntos que venimos comentando. Ahora bien, los grafos lineales presentan diferencias sustanciales con los subconjuntos del plano empleados en el capítulo anterior. Primeramente, no todos los puntos de LG presentan las mismas propiedades locales. Por ejemplo, el número de aristas conectadas en cada nodo de G puede depender del nodo en cuestión. La segunda gran diferencia es la métrica empleada en el grafo. En el plano hemos empleado la distancia euclídea, es decir, la distancia entre dos puntos es la longitud del segmento que los une. El problema surge en que la gran mayoría de los grafos no son convexos, careciendo de sentido emplear la métrica euclídea para medir distancias sobre ellos. La forma apropiada de medir distancias en un grafo lineal es empleando la métrica del camino más corto. Para denir este concepto de forma precisa necesitamos el concepto de camino [8]: Denición 3.2. Sea G= (V,A) un grafo lineal. Dados v,w∈ V , un camino en G entre ellos es un conjunto de nodos {v1,...,vs}⊂V tal que v1=v , vs=w y para todo i∈ {1, . . . , s−1}∃li∈ A conectando vi y vi+1 . Denimos la longitud de este camino como Ps−1 i=1 kvi+1 −vik. Un ciclo en G será un camino en G que empieza y termina en el mismo nodo. Al igual que con el concepto de grafo lineal, existe una dualidad en la forma de representar un camino en un grafo: bien como un conjunto ordenado de nodos o como el conjunto de puntos del plano que forman una poligonal entre los nodos inicial y nal. Podemos denir entonces la métrica del camino más corto en un grafo lineal como: Denición 3.3. Sea G= (V,A) un grafo lineal. Dados v,w∈ V , denimos la distancia del camino más corto (o simplemente la distancia) entre v y w como la menor de las longitudes de los caminos en G entre v y w . La denotaremos por dG(v,w) . Notemos que dG está bien denida, ya que el conjunto de caminos entre dos nodos es nito al serlo V . El cálculo del camino más corto entre dos nodos de un grafo lineal es un problema de programación. El algoritmo más extendido para el cálculo de la distancia del camino más corto entre dos nodos de un grafo lineal es el denominado algoritmo de Dijkstra, cuya formulación puede consultarse en [8]. Una demostración de que este algoritmo en efecto encuentra el camino
3.1. Procesos puntuales en grafos lineales 21 más corto entre dos nodos de un grafo lineal en un tiempo polinomial puede encontrarse en [9]. Ahora bien, en múltiples ocasiones vamos a necesitar calcular la distancia (más corta) entre dos puntos p,q∈ LG que no sean necesariamente nodos de G . Para ello, lo que haremos será añadir p y q al conjunto de nodos de G , modicando también su conjunto de aristas, obteniendo un nuevo grafo e G=e V,e A tal que e V=VS{p,q} y LG=Le G . Hecho esto, calculamos la distancia entre p y q como de G(p,q) empleando el algoritmo de Dijkstra. El procedimiento para añadir un nuevo nodo a un grafo lineal puede verse en detalle en [22]. 3.1. Procesos puntuales en grafos lineales Al examinar patrones de puntos en un grafo lineal, surgen las mismas preguntas que cuando el patrón de puntos se disponía en una región del plano: ¾se distribuyen los puntos uniformemente en el grafo?, ¾depende la densidad de puntos de alguna variable explicativa?, ¾depende la densidad de puntos cerca de un nodo del número de aristas que se conectan en él?. Ya sabemos que estas preguntas no aluden al patrón de puntos en sí, sino al proceso responsable de su generación. Debemos por ello introducir el objeto estadístico que nos permita estudiar patrones de puntos en un grafo. Estos mecanismos serán los procesos puntuales en grafos lineales: Denición 3.4. Un proceso puntual en un grafo lineal G= (V,A) es un mecanismo aleatorio que genera localizaciones distribuidas en LG⊂R2 , y cuya realización es un patrón de puntos. Denotaremos los procesos puntuales en un grafo G por letras mayúsculas P ; mientras que denotaremos los patrones de puntos en este grafo por p ={p1,...,pN}⊂LG⊂R2 . Empleamos por defecto la letra p para denotar los elementos de un patrón de puntos en un grafo p para enfatizar que, aunque sigan siendo puntos del plano, están constreñidos a encontrarse en un subconjunto 1-dimensional del plano como es LG . Es por este mismo motivo que los procesos puntuales en un grafo los denotaremos por defecto por P en vez de por X . En la Figura 3.1 encontramos ejemplos de patrones de puntos en grafos. En (a) se muestra la ubicación de telas de araña a lo largo de las juntas de una pared de ladrillo. El material del que están compuestos los ladrillos y el mortero hacen que las telas de araña solo puedan encontrarse en las juntas, por lo que este conjunto de datos ha de entenderse como un patrón de puntos en un grafo lineal. En (b) encontramos espinas dendríticas en un fragmento de la red dendrítica de una neurona. Este es un ejemplo donde se emplea un grafo lineal para modelar un espacio unidimensional, similar al ejemplo de avistamiento de aves en el litoral. Al igual que en los procesos puntuales en el plano, trabajaremos con procesos puntuales en
22 3. Procesos puntuales en grafos lineales (a) (b) (c) Figura 3.1: ejemplos de patrones puntuales en grafos lineales extraídos del paquete [4]. grafos lineales nitos, que son aquellos tales que cualquier realización p de P es un patrón de puntos nito en el grafo; y tales que el número de puntos observado en cualquier subconjunto L0 G⊂ LG , lo que denotamos por N( P ∩L0 G) , es una variable aleatoria bien denida. Para que los razonamientos que haremos en secciones posteriores sean correctos, todas las subregiones 1 L0 G⊂ LG que consideremos serán conexas y compactas como subconjuntos de R2 , de tal forma que pueden identicarse con un grafo lineal G0 tal que L0 G=LG0 . Además, diremos que dos subregiones de LG son disjuntas si a lo sumo se intersecan en un conjunto de medida nula. Al igual que cuando estudiamos procesos puntuales en el plano, debemos empezar deniendo los procesos puntuales en grafos lineales en los que existe el mayor grado de aleatoriedad. Al igual que en los procesos puntuales en el plano, estos procesos se denominan de aleatoriedad espacial completa (CSR), y se construyen en base a las hipótesis siguientes: Homogeneidad: los puntos no tienen preferencia por ninguna subregión de LG . Independencia: dadas dos regiones disjuntas, las observaciones en una de ellas carecen de inuencia sobre las de la otra. Orden: hay una probabilidad despreciable de que una región sucientemente pequeña contenga más de un punto. De forma más precisa: l´ım |L0 G|→0 1 |L0 G|PN( P ∩L0 G)≥2= 0,∀L0 G⊂ LG, donde |L0 G| denota la longitud total del grafo L0 G (la suma de las longitudes de sus aristas). Al igual que en procesos puntuales en el plano, la hipótesis de homogeneidad implica que el número esperado de puntos observados en una subregión LG ha de ser proporcional a su longitud: ∃λ∈R+ tal que EN( P ∩L0 G)=λ|L0 G|,∀L0 G⊂ LG. 1 Que inicialmente no tienen que ser subgrafos en el sentido estricto de la teoría de grafos lineales.
3.2. Función de intensidad 23 Como ya sabemos, λ recibe el nombre de intensidad del proceso puntual. Por su parte, la hipótesis de independencia garantiza que dadas dos subregiones L0 G,L00 G⊂ LG disjuntas, las variables aleatorias N( P ∩ L0 G) y N( P ∩ L00 G) son independientes. Un razonamiento análogo al realizado en el capítulo anterior permite demostrar que si P es un proceso puntual CSR en un grafo G con intensidad λ , entonces la variable aleatoria N( P ∩ L0 G) sigue una distribución de Poisson de parámetro µ=λ|L0 G| , para todo L0 G⊂ LG . Por ello, a los procesos puntuales CSR en un grafo se les denomina también procesos puntales de Poisson homogéneos en un grafo lineal. La primera desviación de los procesos puntuales CSR en un grafo la encontramos si prescindimos de la hipótesis de homogeneidad. Surgen así los procesos puntuales de Poisson inhomogéneos en un grafo. Estos se caracterizan porque la intensidad es ahora una función de la posición en el grafo λ(p) . Dada una subregión L0 G⊂ LG , podemos realizar un razonamiento análogo al de la Sección 2.1.2, dividiéndola en regiones de longitud arbitrariamente pequeña, lo que permite concluir que si P es un proceso puntual de Poisson inhomogéneo en un grafo G , entonces: EN( P ∩L0 G)=ZL0 G λ(p)dp,∀L0 G⊂ LG, (3.1) donde por dp denotamos el diferencial de línea en L0 G . Recordando que las hipótesis de independencia y orden siguen vigentes en este tipo de procesos puntuales, y la reproductividad de la distribución de Poisson, tenemos que N( P ∩L0 G) sigue una distribución de Poisson de parámetro µ=RL0 Gλ(p)dp , de acuerdo con la ecuación (3.1), para todo L0 G⊂ LG . Al igual que ocurría en los procesos puntuales en el plano, la realización de un proceso puntual en un grafo lineal puede no consistir únicamente en las ubicaciones espaciales de los puntos en el grafo, sino de más información acompañando a estas. Como ya sabemos, esta información adicional vinculada a los eventos se conoce como una marca del proceso puntual. Un ejemplo de una realización de un proceso puntual marcado en un grafo lineal lo encontramos en la Figura 3.1, (c), donde presentamos eventos criminales en el mapa de calles de un barrio de Chicago. El tipo de crimen (asalto, robo, allanamiento, robo de vehículo...) constituye una marca al proceso puntual. Debemos recordar que las marcas son intrínsecas al proceso puntual, y por ello sus valores se obtienen como parte de la realización de este. Por su parte, aquellas magnitudes cuyos valores se conocen a lo largo del grafo, pero que son extrínsecas al proceso puntual, se conocen como covariables del proceso puntual en el grafo. Al igual que ocurría en el plano, una covariable a un proceso puntual P en un grafo lineal G se representa mediante una función Z:LG→R . 3.2. Función de intensidad Como ya sabemos, todo proceso puntual queda caracterizado en términos de primer orden por su función de intensidad, y los procesos puntuales en grafos lineales no son una excepción. Dado
24 3. Procesos puntuales en grafos lineales P un proceso puntual en un grafo lineal G , se dene su función de intensidad λ:LG→[0 ,∞) como una función de la posición en el grafo vericando la ecuación (3.1). A la vista de esta ecuación, podemos interpretar la función de intensidad de un proceso puntual en un grafo lineal como el número esperado de puntos por unidad de longitud. Supongamos que poseemos un patrón de puntos p en un grafo lineal G , y que deseamos emplear la información que este nos brinda para estimar la intensidad del proceso puntual P en G que lo genera. Si tenemos evidencias de que P es un proceso homogéneo, su intensidad no varía a lo largo de LG . Recordando que la intensidad de un proceso puntual en un grafo lineal se interpreta como el número esperado de puntos por unidad de longitud, en [3] se propone el estimador: b λ=N( p ) |LG|, (3.2) que, al igual que ocurría en el plano, es un estimador insesgado de la intensidad de P . En el caso de que no tengamos evidencias signicativas a cerca de la homogeneidad de P , debemos asumir que su intensidad varía a lo largo de LG . Al carecer de ninguna información a cerca de la forma funcional de λ , la estrategia usual empleada será la estimación no paramétrica. Una primera idea podría ser no tener en cuenta el hecho de que nuestro patrón de puntos solo puede observarse sobre LG . Al estar este conjunto embebido en el plano euclídeo, podríamos tomar las coordenadas espaciales de los puntos de p , entender este patrón de puntos como un patrón de puntos en una región W⊂R2 que contenga a LG , y emplear las herramientas de estimación discutidas en el Capítulo 2. Este tipo de técnicas lleva a resultados incorrectos: si por ejemplo P posee intensidad constante, el procedimiento de estimación que hemos descrito tenderá a sobreestimar la intensidad en aquellas zonas donde el grafo sea más denso, y a subestimarla donde la concentración de aristas sea menor. Para ejemplicar este problema, podemos comparar las estimaciones de la intensidad (asumiendo homogeneidad) empleando el patrón de puntos dado en la Figura 3.1, (b). Empleando el estimador dado en la ecuación (3.2) obtenemos b λ= 0.06 µ m −1 ; mientras que si obviamos la estructura de grafo lineal sobre la que se encuentran nuestros puntos y estimamos la intensidad empleando el estimador dado en la ecuación (2.2) concluímos que b λ= 0.00255 µ m −2 , que vemos es una subestimación considerable. Este problema no afecta solo a la función de intensidad. Si empleamos únicamente las coordenadas espaciales de un patrón de puntos en un grafo para estudiar la correlación mediante las distancias entre pares de puntos, comúnmente se concluirá que existe un fenómeno de clusterización a cortas distancias (debido a que los puntos se disponen sobre las aristas del grafo) y homogeneidad a largas distancias (debido a la distancia entre aristas del grafo). Estos fenómenos se han observado en datos reales, como puede verse [28] en más detalle.
3.2. Función de intensidad 25 Otra estrategia más elaborada consiste en dar cabida a la restricción espacial impuesta por el grafo sustituyendo la distancia euclídea por la distancia del camino más corto en los estimadores tipo núcleo de la función de intensidad. Para ello, podríamos plantearnos sustituir las funciones núcleo L(y−xi) empleadas en las estimaciones dadas en las ecuaciones (2.3) y (2.4), por una función núcleo apropiada para la estimación en una dimensión con argumento la distancia del camino más corto en el grafo: k[dG(p,pi)] , con p∈ LG y pi∈ p . Indistintamente, esta tampoco es una estrategia fructuosa. Como puede consultarse en [22], la función núcleo inducida en el grafo k[dG(·,pi)] no es siquiera una función de densidad en LG , al no estar normalizada. La razón fundamental de ello es que, a diferencia de lo que ocurre en R , en un grafo el número de puntos a una distancia de uno dado no tiene que ser siempre dos, debido a la no homogeneidad de los grafos que ya hemos comentado. 3.2.1. Estimadores tipo núcleo equitativos Una posible solución al problema de la función núcleo es emplear estimadores de la intensidad que no solamente empleen la distancia del camino más corto, sino que tengan en cuenta la estructura de LG como subconjunto 1-dimensional del plano. Dado un punto p∈ LG y un ancho de banda h denimos Lp={q∈ LG:dG(p,q)≤h} . En [22] se propone considerar como estimador de la función de intensidad: b λ(p) = N X i=1 Gpi,h(p), (3.3) donde, dado q∈ LG , se dene la función núcleo en LG centrada en q para el ancho de banda h como una función Gq,h :LG→R vericando que: Gq,h(p)≥0∀p∈ LG, Gq,h(p) = 0 ∀p/∈ Lq y ZLq Gq,h(p)dp= 1. (3.4) Para que el estimador propuesto en la ecuación (3.3) sea un estimador insesgado de la función de intensidad, en [22] se propone construir las funciones núcleo Gq,h en términos de la denominada función núcleo básica w . Estas son funciones de la distancia más corta desde q a través de LG , y sobre las que asumiremos las siguientes hipótesis: (Ri)RLGw[dG(q,p)]dq= 1 para todo p∈ LG . (Rii)w es una función no negativa, continua y no creciente respecto de la distancia del camino más corto entre dos puntos de LG . (Riii) Fijado el ancho de banda h , w[dG(q,p)] = 0 siempre que dG(q,p)≥h .
32 4. Comparación de funciones de intensidad Ahora bien, en el caso de que los procesos puntuales ante los que nos encontremos sean procesos puntuales denidos en un grafo lineal, aún no se han descrito este tipo de técnicas para comparación de funciones de intensidad. En este capítulo trataremos entonces de desarrollar procedimientos, innovadores en el campo, que permitan comparar las funciones de intensidad de dos procesos puntuales denidos en un mismo grafo lineal. Sean P 1 y P 2 dos procesos puntuales en un grafo lineal G con funciones de intensidad λ1 y λ2 , respectivamente. De estos procesos puntuales poseemos sendas realizaciones: p 1 y p 2 , patrones de puntos en LG . Denotamos Ni= # p i , con i= 1,2 . Nuestro objetivo es entonces estudiar si estos dos procesos puntuales poseen la misma estructura espacial. Notemos que esto no se traduce en que tengan la misma función de intensidad, ya que estas podrían estar en distintas escalas; es decir, el número esperado de puntos podría ser diferente, aunque su distribución espacial fuera la misma. Lo que vamos a contrastar es entonces que estos procesos puntuales tengan la misma función de densidad relativa, ya que notemos que estas están normalizadas en LG , lo que es equivalente a que sus funciones de intensidad sean proporcionales. Así, la hipótesis nula es: H0:∃η > 0 tal que λ1(p) = ηλ2(p)∀p∈ LG. (4.1) Antes de proceder con la construcción de los distintos estadísticos para contrastar la hipótesis nula dada en (4.1), debemos hacer un inciso. Estrictamente hablando, si queremos formular la hipótesis nula de igualdad de funciones de densidad para dos procesos puntuales en un grafo cualquiera, debemos hacerlo a través de la proporcionalidad de las funciones de intensidad condicionadas, λc(p|q) , que representa la intensidad del proceso puntual en p∈ LG condicionada a observar un punto en q∈ LG . Ahora bien, bajo la hipótesis de que nuestros procesos sean procesos puntuales de Poisson, la hipótesis de independencia nos dice que λc(p|q) = λ(p) para todo p,q∈ LG . Por ello, si asumimos que trabajamos con procesos de Poisson, podemos formular la hipótesis nula que queremos contrastar tal y como se detalla en (4.1). 4.1. Test de Kolmogorov-Smirnov La hipótesis nula que queremos contrastar establece una relación de proporcionalidad punto a punto entre funciones de intensidad. Traduciremos esta igualdad a una serie de subconjuntos L0 G⊂ LG , y estudiaremos en qué grado de verica. Integrando la igualdad dada en (4.1) en L0 G , tenemos que: EN( P 1∩L0 G)=ηEN( P 2∩L0 G),∀L0 G⊂ LG, (4.2) donde únicamente permitimos subconjuntos de LG que pertenezcan a su σ -álgebra de Borel (la σ -álgebra en LG generada por todos los abiertos en la topología inducida por la métrica del camino más corto). La ecuación (4.2) motiva una estimación global de η como bη=N1/N2 . Si
4.1. Test de Kolmogorov-Smirnov 33 denotamos Ni(L0 G) = #( p i∩L0 G) para todo L0 G⊂ LG e i= 1,2 , inspirándonos en [30] podemos denir una medida de discrepancia como: D(L0 G) = N1(L0 G)−bηN2(L0 G)=N1 N1(L0 G) N1−N2(L0 G) N2,∀L0 G⊂ LG. (4.3) Una propuesta de estadístico para determinar la discrepancia de la hipótesis nula sería el supremo de todas las posibles discrepancias dadas en la ecuación (4.3), con L0 G⊂ LG . El problema es que la σ -álgebra de Borel de LG es una familia de subconjuntos demasiado grande. Por ello, para procesos puntuales en el plano euclídeo, en [30] se encuentra suciente tomar el supremo de las discrepancias dadas en (4.3) entre un π -sistema de W que genere la σ -álgebra de Borel. Denición 4.1. Sea Φ un conjunto, diremos que P ⊂ P (Φ) es un π -sistema de Φ si P 6=∅ y A∩B∈ P para todo A, B ∈ P . Denimos también la σ -álgebra generada por P como la menor σ -álgebra de Φ que lo contiene. Finalmente, diremos que un π -sistema en Φ genera una σ -álgebra en Φ si esta está contenida en la σ -álgebra generada por el π -sistema. Así, inspirándonos en el estadístico dado en [30], tomamos P un π -sistema en LG que genere su σ -álgebra de Borel, y proponemos como estadístico para contrastar la hipótesis nula (4.1): TKS =1 ξrN1N2 N1+N2 sup L0 G∈P N1(L0 G) N1−N2(L0 G) N2, (4.4) donde ξ es una constante normalizadora. En [30], la constante análoga a ξ garantiza la convergencia del estadístico a un puente Browniano cuando se estudian procesos puntuales en el plano. Sin tener ningún resultado teórico acerca de la convergencia del estadístico propuesto en la ecuación (4.4), optamos por generalizar la estimación de ξ dada en [30] como sigue: sea {Lγ}Γ γ=1 una partición de LG tal que en todo elemento de la partición se observa al menos un punto de alguno de los dos patrones observados, y denimos, para todo γ∈ {1,...,Γ} e i= 1,2 : c Ni(Lγ) = Ni·N1(Lγ) + N2(Lγ) N1+N2 , que no es más que una estimación del número de puntos observados en cada elemento de la partición para cada patrón de puntos. Usando estas cantidades estimamos entonces ξ como: b ξ= 1 Γ−1 Γ X γ=1 "[N1(Lγ)−c N1(Lγ)]2 c N1(Lγ)+[N2(Lγ)−c N2(Lγ)]2 c N2(Lγ)# 1/2 , (4.5) que representa una medida de discrepancia entre el número de puntos observados en cada patrón para cada elemento de la partición, y los valores estimados empleando los dos patrones de forma conjunta. En efecto, vemos como cada uno de los sumandos de (4.5) es de la forma de un estadístico χ2 . Notemos que (4.5) está bien denido ya que c Ni(Lγ)>0 para todo i= 1,2 y
34 4. Comparación de funciones de intensidad γ∈ {1,...,Γ} , ya que por construcción hemos garantizado que en todos los elementos de la partición se observa al menos un punto de alguno de los dos patrones observados. A la hora de calcular el estadístico dado en la ecuación (4.4), han de elegirse tanto el π - sistema, P , como la partición del grafo. Una partición natural del grafo es su conjunto de aristas A . Ahora bien, este podría no ser una partición válida, ya que podría darse que en alguna arista no se observara ningún punto en ninguno de los dos patrones. Esto no es un gran problema, ya que podemos tomar cada una de las aristas en las que no se observe ningún punto y unirlas, como subconjunto de LG , a una arista en la que si se observe algún punto, obteniendo así una partición de LG válida. Notemos que esto es equivalente a restringir la suma dada en la ecuación (4.5) a aquellas aristas en las que se observe algún punto. La elección del π -sistema, P , no es una cuestión tan directa. La opción que proponemos tiene como principal objetivo simplicar la computación del estadístico propuesto en la ecuación (4.4). Proponemos tomar P como el conjunto de bolas en LG , centradas en uno de sus puntos, respecto a la métrica del camino más corto en LG , es decir: P={BLG(q, r): r > 0} donde BLG(q, r) = {p∈ LG:dG(p,q)< r}. Para denir P es necesario elegir el punto q a partir del cual se construye el π -sistema empleado en el estadístico. La elección de este punto depende del tipo de grafo ante el que nos encontremos. Como por defecto trabajaremos con grafos conexos, podemos distinguir si estos grafos no son o son árboles, es decir, si poseen o no ciclos. Un ejemplo de grafo con ciclos es el dado en la Figura 3.1 (a). En este tipo de grafos una elección razonable es un punto centrado en el grafo. En caso de que nuestro grafo sea un árbol, como por ejemplo el que vemos en la Figura 3.1 (b), resulta natural tomar como punto base para la construcción de P el nodo raíz. 4.2. Test de Cramer von Mises Para el procedimiento que describiremos en esta sección, debemos recordar que la hipótesis nula formulada en (4.1) es equivalente a la igualdad de las funciones de densidad relativas. Por ello, bajo la hipótesis nula, la distancia (en un determinado espacio de funciones) entre estas funciones de densidad ha de ser cero. Así, una medida de discrepancia respecto de la hipótesis nula puede ser la distancia en L 2 entre las estimaciones de las funciones de densidad. Recordemos que si b λi(p) es una estimación de la función de intensidad de P i como las estudiadas en el Capítulo 3, con i= 1,2 , entonces b λi(p)/Ni es una estimación de su función de densidad.
4.3. Test no paramétrico basado en la función de riesgo relativo 35 Proponemos entonces como estadístico de contraste para la hipótesis nula dada en (4.1): TCvM =ZLG"b λ1(p) N1−b λ2(p) N2#2 dp. (4.6) Recordemos como, bajo la hipótesis nula, las funciones de densidad asociadas a ambos procesos puntuales son iguales. Por ello, mayores valores del estadístico reejan mayor discrepancia con la hipótesis nula. A la hora de computar el estadístico propuesto en la ecuación (4.6), surge la cuestión a cerca de cómo de similares han de ser las estimaciones de la intensidad de cada uno de los procesos puntuales. Idealmente, nos gustaría que estas se construyeran de la forma más parecida posible, con el objetivo de que las discrepancias detectadas por el estadístico se deban a diferencias en la estructura espacial de los procesos puntuales estudiados, y no a diferencias causadas por los distintos métodos de estimación. Por ejemplo, si empleamos estimaciones no paramétricas como las estudiadas en el Capítulo 3, parece conveniente emplear en ambos casos el mismo tipo de núcleos equitativos, o en ambos casos la estimación basada en la ecuación del calor. Una vez decididos los estimadores a emplear, el cómputo de este estadístico TCvM no necesita de la elección de otros elementos, como sí ocurría en el test de Kolmogorov-Smirnov presentado en la Sección 4.1. Si bien es cierto que esta ventaja del aligeramiento en cuanto al proceso se ve difuminada por un mayor coste computacional, tal y como explicaremos en el Capítulo 5. 4.3. Test no paramétrico basado en la función de riesgo relativo Otra forma alternativa de reescribir la hipótesis nula dada en (4.1) es emplear la función de riesgo relativo (cociente entre intensidades): H0:r(p) = λ1(p) λ2(p) es constante a lo largo de LG. (4.7) Así, trataremos de contrastar si el cociente de intensidades varía a lo largo del grafo. Ahora bien, notemos que debido al posible distinto número esperado de puntos observados en la realización de cada uno de los procesos puntuales, no conocemos el valor de la constante a la que iguala la función de riesgo relativo. Por ello, resulta más apropiado contrastar que el cociente de las funciones de densidad de cada proceso puntual en el grafo sea constate e igual a 1, o equivalentemente que su logaritmo sea nulo en todo el grafo. Por ello, estimamos el logaritmo de la función de riesgo relativo como: bρ(p) = ln "N2 N1b λ1(p) b λ2(p)#.
36 4. Comparación de funciones de intensidad Así, podemos formular un test de efecto de p sobre bρ(p) : si reindexamos los puntos observados como p 1S p 2={pj}N j=1 , con N=N1+N2 , podemos entender los valores de {bρ(pj)}N j=1 como una muestra de una variable respuesta frente a la variable explicativa dada por la posición en el grafo, con valores asociados {pj}N j=1 . Tenemos entonces un problema de regresión en el grafo, que nos permite reformular el test con hipótesis nula dada en (4.7) como el test de efecto siguiente: H0:Ebρ(pj)|pj=µ frente a Ha:Ebρ(pj)|pj=m(pj), (4.8) siendo m una función a lo largo del grafo, usualmente denominada función de regresión. Necesitamos entonces herramientas de regresión en grafos lineales. Este es otro de los temas que no han sido abordados previamente en la literatura existente. Por ello, vamos a desarrollar un método no paramétrico de estimación de la función de regresión en un grafo lineal. Si en un contexto de regresión tenemos n valores de la variable respuesta {yj}n j=1 y n valores de una variable explicativa {xj}n j=1 , una estimación no paramétrica de la función de regresión es la dada por el estimador de Nadaraya-Watson que, tal y como se detalla en [27], podemos escribir como: bm(x) = Pn j=1 Sh(x−xj)yj Pn j=1 Sh(x−xj), (4.9) donde S es una función núcleo, h es el parámetro de ventana, y Sh(·) = h−1S(·/h) . Las condiciones sobre S pueden consultarse en [27], aunque los núcleos habitualmente empleados en la estimación no paramétrica de la función de densidad siguen siendo válidos en este contexto. Ahora bien, cuando nos encontramos ante un proceso puntual en un grafo lineal, la estimación de la función de densidad de este proceso ha de hacerse empleando las técnicas estudiadas en el Capítulo 3. Si empleamos núcleos equitativos, podemos escribir una estimación de la función de densidad de un proceso puntual en un grafo mediante una estimación de su intensidad como: 1 Nb λ(p) = 1 N N X j=1 Gpj,h(p), (4.10) siendo {pj}N j=1 un patrón de puntos en el grafo, y Gpj,h el núcleo equitativo centrado en pj con ancho de banda h . Comparando entonces los estimadores de la función de densidad en R y en un grafo lineal, e inspirándonos en el estimador de Nadaraya-Watson dado en la ecuación (4.9), proponemos un estimador no paramétrico de la función de regresión en un grafo lineal como: bm(p) = PN j=1 Gpj,h(p)yj PN j=1 Gpj,h(p), (4.11) siendo {yj=bρ(pj)}N j=1 el conjunto de observaciones de la variable respuesta, asociados a una serie de localizaciones en el grafo {pj}N j=1 ⊂ LG , y donde Gpj,h sigue siendo un núcleo equitativo centrado en pj con ancho de banda h .
4.3. Test no paramétrico basado en la función de riesgo relativo 37 En el contexto del problema de regresión ante el que nos encontramos, tenemos N=N1+N2 , los valores observados de la variable respuesta son {bρ(pj)}N j=1 asociados a las localizaciones en el grafo {pj}N j=1 = p 1S p 2⊂ LG . De acuerdo con (4.8), el modelo bajo la hipótesis nula puede estimarse mediante la media muestral: bµ=1 N N X j=1 bρ(pj). Por otra parte, la estimación de la función de regresión bajo la hipótesis alternativa emplea el estimador no paramétrico dado en (4.11): bm(pi) = PN j=1 Gpj,h(pi)bρ(pj) PN k=1 Gpk,h(pi)= N X j=1 Si,j bρ(pj), donde hemos denido la matriz de suavizado S∈ MN×N(R) como: Si,j =Gpj,h(pi) PN k=1 Gpk,h(pi), i, j ∈ {1, . . . , N}. Como el contraste de efecto propuesto en (4.8) puede entenderse como un contraste entre modelos anidados (dónde el modelo bajo la hipótesis nula es un caso particular del modelo bajo la hipótesis alternativa), emplearemos un test F. Denimos los residuos del modelo bajo la hipótesis nula como: RSS 0= N X j=1 bρ(pj)−bµ2, y bajo la hipótesis alternativa como: RSS a= N X j=1 bρ(pj)−bm(pj)2. Como bajo la hipótesis nula estimamos un único parámetro, el número de grados de libertad de los residuos bajo la hipótesis nula es df0=N−1 . Para denir los grados de libertad de los residuos bajo la hipótesis alternativa seguimos la estrategia propuesta en [6], donde, en analogía con el modelo lineal general, se propone tomar dfa=tr(IN−S) . Entonces el estadístico de contraste para el test de efecto viene dado por: TNP =( RSS 0− RSS a)/(df0−dfa) RSS a/dfa . (4.12) Al igual que ocurría en el test de Cramer von Mises presentado en la Sección 4.2, para computar el valor del estadístico propuesto en la ecuación (4.12) debemos elegir cómo vamos a llevar a cabo las distintas estimaciones no paramétricas requeridas. Reparamos primeramente en
38 4. Comparación de funciones de intensidad el cálculo de {bρ(pj)}N j=1 ; debemos procurar de nuevo que las estimaciones de la intensidad de cada uno de los procesos puntuales sean lo más similares en lo que al proceso de estimación se reere, por ello, es recomendable emplear el mismo tipo de núcleos equitativos, o la estimación basada en la ecuación del calor, en ambas estimaciones. Además, para que los valores de {bρ(pj)}N j=1 estén bien denidos, necesitamos que los anchos de banda sean sucientemente grandes como para que la estimación de las funciones de intensidad de los procesos puntuales en estudio no se anulen en ninguno de los puntos de los dos patrones observados. Esto hace que no podamos tomar anchos de banda demasiado pequeños, lo que aumentará el coste computacional del cálculo de los valores {bρ(pj)}N j=1 , ya el coste del cálculo de los núcleos equitativos aumenta con h . Nótese que únicamente vamos a necesitar los valores de las estimaciones de la intensidad en los puntos de los patrones observados. Por ello, puede ser recomendable determinar estas cantidades empleando validación cruzada; es decir: bρ(pj) = ln N2 N1−1 b λ1,−j(pj) b λ2(pj) si pj∈ p 1 ln N2−1 N1 b λ1(pj) b λ2,−j(pj) si pj∈ p 2 , donde b λi,−j es la estimación de la función de intensidad de P i empleando como patrón de puntos p i\{pj} , con j∈ {1, . . . , N} e i= 1,2 . A la hora de computar la estimación no paramétrica de la función de regresión en el grafo, debemos elegir nuevamente qué tipo de núcleos equitativos vamos a usar, así como qué ancho de banda, ya que este no tiene porqué coincidir con ninguno de los empleados anteriormente en el cálculo de los valores {bρ(pj)}N j=1 . 4.4. Procedimiento de calibración A lo largo de este capítulo hemos descrito tres procedimientos para construir estadísticos de contraste sobre la hipótesis nula dada en (4.1). Ahora bien, para armar si a un determinado nivel de signicación tenemos evidencias signicativas a favor o en contra de la hipótesis nula, el valor del estadístico de contraste no es suciente, necesitamos conocer su nivel crítico [16]: Denición 4.2. Sea el contraste de hipótesis H0:θ∈Θ0 frente a Ha:θ∈Θa , para el cual se posee un estadístico conveniente T(x1, . . . , xs) , siendo (x1, . . . , xs) una muestra. Se denomina nivel crítico o p-valor de una muestra (y1, . . . , yl) para el contraste en cuestión mediante el estadístico T a: p(y1, . . . , yl) = sup θ∈Θ0 P[T(x1, . . . , xs)≥T(y1, . . . , yl)|θ],
4.4. Procedimiento de calibración 39 donde P(A|B) denota la probabilidad del suceso A condicionada al suceso B . Notemos que carecemos de resultados teóricos acerca de la distribución de nuestros estadísticos bajo la hipótesis nula que nos permitieran calcular p-valores de manera sencilla a partir de dicha distribución. Por ello, una posible solución es emplear el test de permutaciones adaptado al problema de las dos muestras. Considerados dos patrones de puntos que entendemos como realizaciones de sendos procesos puntuales en un grafo, se calcula uno de los tres estadísticos de contraste que hemos propuesto. Para estimar su nivel crítico, se consideran todas las posibles permutaciones de los puntos de los patrones observados, tomando dos subconjuntos de {pj}N j=1 de cardinales N1 y N2 . Para cada una de estas permutaciones se calcula el valor del estadístico de contraste. Hecho esto, el p-valor se estima como la fracción de estos estadísticos que son mayores que el de partida (el asociado a los patrones de puntos sin permutar). En el caso de que ninguno de estos estadísticos sea mayor que el de partida (lo que conduciría a una estimación nula del nivel crítico) consideraremos como estimación del p-valor la mitad del inverso del número de permutaciones realizadas. En la práctica este procedimiento presenta el problema de que el número total de permutaciones es computacionalmente inabordable a partir de tamaños muestrales no muy grandes. Por ello, una práctica estándar es considerar un subconjunto de np permutaciones, escogidas al azar, de entre todas las posibles. La elección de np no es cuestión menor, ya que debe escogerse sucientemente grande como para dar resultados ables, al tiempo que ha de ser computacionalmente realizable. Como se discute en [19], tomar np= 5000 resulta suciente en la mayoría de escenarios, llegando en casos excepcionales a tomar np= 10000 permutaciones. Por ello, de partida jamos np= 7500 .
40 4. Comparación de funciones de intensidad
Capítulo 5 Estudio de simulación En el Capítulo 4 hemos diseñado procedimientos para contrastar si dos patrones de puntos, observados en un mismo grafo, están generados por procesos puntuales con funciones de intensidad proporcionales. En el presente capítulo presentamos el exhaustivo estudio de simulación llevado a cabo para analizar el comportamiento de nuestras propuestas. Al tratarse de procedimientos de contraste, el estudio de su comportamiento radica fundamentalmente en valorar dos cuestiones: el nivel (probabilidad de rechazar cuando la hipótesis nula es cierta) y la potencia (la probabilidad de rechazar cuando la hipótesis nula es falsa). Para evaluar estas cantidades emplearemos técnicas de Monte Carlo: jado un modelo, generaremos M= 1000 realizaciones (réplicas Monte Carlo) de cada uno de los dos procesos puntuales considerados, y determinaremos la conclusión de cada uno de los procedimientos de contraste empleando el test de permutaciones introducido en la Sección 4.4. (a) Spiders . (b) Dendrite . Figura 5.1: representación de los grafos Spiders (a) y Dendrite (b) que consideraremos en el estudio de simulación. 41
48 5. Estudio de simulación 5.3.2. Test de Cramer von Mises En primer lugar notemos la distinta naturaleza entre el estadístico de Cramer von Mises y el de Kolmogorov-Smirnov. El test de Cramer von Mises no calcula distancias entre pares de puntos, sino que mide distancias entre estimaciones de la función de intensidad. En un primer paso, necesitamos decidir qué estimación vamos a emplear. Recordando las consideraciones hechas en la Sección 4.2, hemos optado por emplear en todos los casos la estimación basada en la ecuación del calor, introducida en la Sección 3.2.2, escogiendo en cada caso el ancho de banda mediante la regla de Scott modicada, denida en la ecuación (3.6). Estas elecciones tienen como principal objetivo reducir el coste computacional del cálculo de este estadístico. Aunque el contraste de Cramer von Mises no está condicionado a que el grafo sea o no un árbol, por consistencia con el estudio realizado para el test de Kolmogorov-Smirnov estudiaremos el nivel y la potencia del mismo tanto en Spiders como en Dendrite . A pesar de las elecciones hechas en lo referente al método de estimación y selección de los correspondientes parámetros, el número elevado de nodos de los grafos propuestos hace computacionalmente inabordable en tiempos factibles el cálculo sobre los mismos de este estadístico. (a) Subgrafo de Spiders . (b) Subgrafo de Dendrite . Figura 5.6: representación de los subgrafos de Spiders (a) y de Dendrite (b) empleados en el estudio del nivel y la potencia para el estadístico de Cramer von Mises. Con el objetivo de reducir este coste computacional, pero manteniendo la esencia de los modelos descritos anteriormente optamos por considerar subgrafos de Spiders y de Dendrite que tengan un menor número de nodos, pero sigan representando la estructura geométrica de los mismos. Los subgrafos que vamos a emplear pueden verse en la Figura 5.6, y se han obtenido limitándonos a un subconjunto (en el plano euclídeo) de los grafos originales. Notemos que el subgrafo de Spiders sigue teniendo ciclos, y que el de Dendrite sigue siendo un árbol. Esta acción no fue suciente, por ello, para reducir todavía más el coste computacional reduciremos el número
5.3. Resultados 49 (a) Subgrafo de Spiders . (b) Subgrafo de Dendrite . Figura 5.7: representación de la función de densidad asociada a la intensidad inhomogénea, λInh , denida en la ecuación (5.2), en los subgrafos de Spiders (a) y Dendrite (b). de permutaciones consideradas a np= 5000 , elección que según [19] no pone en compromiso el proceso de calibrado. Los modelos bajo la hipótesis nula que consideraremos para el estudio del nivel serán los que hemos detallado en la Sección 5.1: el modelo con intensidad homogénea denido en la ecuación (5.1), y en la ecuación (5.2) para el caso inhomogéneo. Las funciones de densidad relativas a estas funciones de intensidad pueden verse representadas en los subgrafos de Spiders y de Dendrite en la Figura 5.7, dónde vemos que el comportamiento espacial es el mismo que el observado en la Figura 5.2 para los grafos originales. En la Tabla 5.5 y Tabla 5.6 mostramos los resultados del calibrado del nivel del estadístico Cramer von Mises para el subgrafo de Spiders y de Dendrite respectivamente. Podemos observar como a partir de tamaños muestrales esperados de orden m= 100 el contraste está bien calibrado, si bien es cierto que para el menor nivel de signicación, α= 0.01 , presenta ciertas limitaciones, que se ven subsanadas al aumentar el tamaño muestral. Para el estudio de la potencia del test de Cramer von Mises consideraremos los modelos bajo la hipótesis alternativa que hemos introducido en la Sección 5.2. De nuevo emplearemos los subgrafos de Spiders y Dendrite que hemos introducido en la Figura 5.6. En este caso los subconjuntos L+ G en los que se añaden en media a puntos homogéneamente vienen dados por: L+ Subgrafo Spiders ={p= (p1, p2)∈ L Subgrafo Spiders :p∈[735.25,984.375] ×[735.25,956.25]}, y L+ Subgrafo Dendrite ={p= (p1, p2)∈ L Subgrafo Dendrite :p2>295.9625}.
50 5. Estudio de simulación Modelo homogéneo Modelo inhomogéneo α= 0.1α= 0.05 α= 0.01 α= 0.1α= 0.05 α= 0.01 m= 20 0.077 0.035 0.001 0.073 0.026 0.002 m= 50 0.076 0.037 0.005 0.087 0.038 0.004 m= 100 0.080 0.046 0.002 0.097 0.05 0.011 m= 200 0.107 0.044 0.003 0.098 0.051 0.006 m= 500 0.090 0.043 0.013 0.092 0.047 0.008 m= 750 0.101 0.044 0.009 0.081 0.038 0.008 Tabla 5.5: proporción de rechazos bajo la hipótesis nula a distintos niveles de signicación, α∈ {0.1,0.05,0.01} , para el test de Cramer von Mises en el subgrafo de Spiders , con modelos homogéneo e inhomogéneo y tamaños muestrales en media m= 20,50,100,200,500 y 750 . Modelo homogéneo Modelo inhomogéneo α= 0.1α= 0.05 α= 0.01 α= 0.1α= 0.05 α= 0.01 m= 20 0.086 0.042 0.005 0.089 0.045 0.007 m= 50 0.091 0.040 0.006 0.088 0.048 0.003 m= 100 0.108 0.048 0.003 0.084 0.041 0.011 m= 200 0.099 0.039 0.013 0.100 0.048 0.010 m= 500 0.106 0.040 0.005 0.102 0.037 0.003 m= 750 0.089 0.048 0.008 0.098 0.048 0.008 Tabla 5.6: proporción de rechazos bajo la hipótesis nula a distintos niveles de signicación, α∈ {0.1,0.05,0.01} , para el test de Cramer von Mises en el subgrafo de Dendrite , con modelos homogéneo e inhomogéneo y tamaños muestrales en media m= 20,50,100,200,500 y 750 . Estas elecciones tienen como objetivo que |L+ G|/|LG| sea aproximadamente el mismo en todos los casos. Atendiendo a estas deniciones, las funciones de intensidad bajo la hipótesis alternativa pueden verse representadas en la Figura 5.8 (caso homogéneo) y en la Figura 5.9 (escenario inhomogéneo). En la Tabla 5.7 y en la Tabla 5.8 se encuentran los resultados del estudio de la potencia del contraste de Cramer von Mises, en los escenarios homogéneo e inhomogéneo respectivamente. Lo primero que debemos notar es que en este caso hemos tomado los mismos valores de a tanto para el subgrafo de Spiders como para el de Dendrite , ya que se alcanzan valores de potencia muy altos en ambos soportes para los mismos valores de a .
5.3. Resultados 51 (a) Subgrafo de Spiders . (b) Subgrafo de Dendrite . Figura 5.8: representación de la función de intensidad homogénea perturbada de acuerdo a la ecuación (5.3) en los subgrafos de Spiders (a) y Dendrite (b) para el caso m= 100 y a= 50 . (a) Subgrafo de Spiders . (b) Subgrafo de Dendrite . Figura 5.9: representación de la función de intensidad inhomogénea perturbada de acuerdo a la ecuación (5.3) en los subgrafos de Spiders (a) y Dendrite (b) para el caso m= 300 . En (a) hemos empleado a= 100 y en (b) a= 10 . Subgrafo de Spiders Subgrafo de Dendrite a= 3m/4a=m/2a=m/4a= 3m/4a=m/2a=m/4 m= 50 0.873 0.557 0.168 0.977 0.843 0.377 m= 100 0.999 0.927 0.378 1 0.986 0.669 m= 200 1 1 0.790 1 1 0.938 m= 500 1 1 0.994 1 1 1 Tabla 5.7: proporción de rechazos bajo la hipótesis alternativa al nivel de signicación α= 0.05 para el test de Cramer von Mises en los subgrafos Spiders y Dendrite , en el caso homogéneo ( λ1=λHom ), para diferentes tamaños muestrales esperados del primer patrón de puntos ( m= 50,100,200 y 500 ) y diferentes valores del número de puntos añadidos de forma esperada homogéneamente en L+ G en el segundo patrón de puntos a .
52 5. Estudio de simulación Subgrafo de Spiders Subgrafo de Dendrite a= 3m/4a=m/2a=m/4a= 3m/4a=m/2a=m/4 m= 50 0.905 0.612 0.184 0.985 0.848 0.394 m= 100 1 0.926 0.414 1 0.992 0.680 m= 200 1 0.999 0.738 1 1 0.994 m= 500 1 1 0.999 1 1 1 Tabla 5.8: proporción de rechazos bajo la hipótesis alternativa al nivel de signicación α= 0.05 para el test de Cramer von Mises en los subgrafos Spiders y Dendrite , en el caso inhomogéneo ( λ1=λInh ), para diferentes tamaños muestrales esperados del primer patrón de puntos ( m= 50,100,200 y 500 ) y diferentes valores del número de puntos añadidos de forma esperada homogéneamente en L+ G en el segundo patrón de puntos a . En ambas tablas vemos como para cada valor de m la potencia aumenta con a . Esto era de esperar, ya que nos alejamos de la hipótesis nula. Nótese además como, jado a en función de m , según aumentamos también m aumenta la potencia, ya que tenemos más información para discernir si las funciones de intensidad de los procesos puntuales que generan los patrones observados son o no proporcionales. Estos comportamientos esperables se observan en los dos grafos que hemos estudiado. Por todo ello, podemos concluir que los resultados del estudio de la potencia para el test de Cramer von mises son satisfactorios.
Capítulo 6 Aplicación a datos reales Los accidentes de carretera se han convertido, en el último medio siglo, en un serio problema de seguridad ciudadana a lo largo de todo el mundo: Brasil [1], India [18], Irán [15], Korea [29], etc. El estudio de la distribución de accidentes de tráco a lo largo de una red de carreteras, analizando los puntos de mayor acumulación de accidentes (puntos negros), o cómo su distribución se ve afectada por factores externos, resulta una labor esencial de cara a mejorar la seguridad vial, y en esencia, salvar vidas. En concreto, el ayuntamiento de Río de Janeiro ha puesto en marcha un plan dedicado a mejorar la seguridad en la circulación por las carreteras de la ciudad. A raíz de esto, hemos conseguido un conjunto de datos reales sobre los que vamos a ilustrar las técnicas desarrolladas en el Capítulo 4 1 . La base de datos cuenta con un total de 270908 entradas, cada una de las cuales corresponde con un accidente de tráco en una carretera de Río de Janeiro entre el 21 de marzo de 2019 y el 4 de mayo de 2022. Estos eventos son reportados por conductores/transeúntes a través de la plataforma Waze 2 . Para cada evento se conocen: sus coordenadas geográcas, la fecha y hora en la que el accidente fue reportado, el tipo de vía en el que tuvo lugar, su dirección postal, e indicadores de la calidad de la medida asociados a posteriores noticaciones de que el accidente se ha reportado de forma correcta. Para poder estudiar la distribución de los accidentes de tráco a lo largo de la red de carreteras de Río de Janeiro, necesitamos un grafo lineal que la describa. Para que este grafo sea computacionalmente manejable, optamos por considerar únicamente las carreteras principales de la ciudad, obviando vías secundarias. Así, para la construcción del grafo emplearemos como referencia el mapa de carreteras de Río de Janeiro que puede verse en la Figura 6.1 (a). Toma1 Agradecemos al profesor Rodrigo S. Targino de la Fundación Getulio Vargas por habernos proporcionado esta base de datos. 2 La página principal de esta plataforma puede verse aquí. 53
54 6. Aplicación a datos reales (a) Mapa de las principales carreteras de Río, extraído de [17]. (b) Grafo lineal representando la red de carreteras principales de Río de Janeiro. Figura 6.1: mapa de las principales carreteras de Río de Janeiro (a), y grafo lineal representativo, construido a partir del anterior (b). (a) Conjunto de datos original. (b) Accidentes ocurridos en vías principales. Figura 6.2: representación sobre el grafo lineal de las coordenadas geográcas de los accidentes de carretera recogidos en la base de datos (a), y de aquellos que tuvieron lugar en vías principales (b). mos como conjunto de nodos los cruces de las carreteras de este mapa, y consideramos luego las aristas que se identican con las carreteras que unen dichos nodos. Obtenemos así el grafo lineal que vemos en la Figura 6.1 (b), que posee 534 nodos y 806 aristas. Notemos que este no deja de ser una aproximación, ya que además de considerar solo las vías principales, aproxima tramos curvos por segmentos rectos. Ahora bien, a pesar de que nosotros hemos considerado únicamente para la construcción de nuestro grafo lineal las carreteras principales de Río de Janeiro, la base de datos de la que disponemos contiene también accidentes en vías secundarias, como puede apreciarse en la Figura 6.2 (a). Como nuestra base de datos incluye una variable categórica que nos indica el tipo de vía en el que se reportó el accidente, de ahora en adelante nos restringiremos a aquellos accidentes que ocurran sobre las carreteras principales: 241804 que eventos pueden verse representados en la Figura 6.2 (b).
55 (a) Patrón de puntos. (b) Estimación de la densidad. Figura 6.3: representación del patrón de puntos de los accidentes ocurridos en vías principales de lunes a viernes (a), junto con una estimación de su función de densidad asociada (b), obtenida empleando la ecuación del calor escogiendo el ancho de banda con la regla de Scott modicada. (a) Patrón de puntos. (b) Estimación de la densidad. Figura 6.4: representación del patrón de puntos de los accidentes ocurridos en vías principales los sábados y domingos (a), junto con una estimación de su función de densidad (b), obtenida empleando la ecuación del calor escogiendo el ancho de banda con la regla de Scott modicada. El primer problema de interés es comparar la distribución espacial de los accidentes que ocurren durante días laborables con los que ocurren en n de semana. Para ello, extraemos los accidentes de lunes a viernes (189996 eventos) y los ocurridos a lo largo del n de semana (51681 eventos). Estos dos patrones, junto con estimaciones no paramétricas de su función de densidad, pueden verse en la Figura 6.3 y en la Figura 6.4. Para este par de patrones de puntos calcularemos el nivel crítico asociado al contraste tanto de Kolmogorov-Smirnov como de Cramer von Mises, empleando np= 10000 permutaciones. Para el cómputo del estadístico de Kolmogorov-Smirnov se ha tomado como punto base del π -sistema en el grafo lineal de las carreteras principales de Río de Janeiro el punto que hemos representado en la Figura 6.7. En ambos casos la estimación del nivel crítico ha sido de 5·10−5 ; es decir, tenemos evidencias signicativas en contra de la hipótesis nula. Podemos por ello concluir que existen evidencias signicativas a favor de que la estructura espacial de los accidentes de tráco en las carreteras principales de Río de Janeiro cambia de los días de semana a los nes de semana.
56 6. Aplicación a datos reales (a) Patrón de puntos. (b) Estimación de la densidad Figura 6.5: representación del patrón de puntos de los accidentes ocurridos en vías principales de 10 a 13 horas (a), junto con una estimación de su función de densidad asociada (b), obtenida empleando la ecuación del calor escogiendo el ancho de banda con la regla de Scott modicada. (a) Patrón de puntos. (b) Estimación de la densidad Figura 6.6: representación del patrón de puntos de los accidentes ocurridos en vías principales de 20 a 23 horas (a), junto con una estimación de su función de densidad (b), obtenida empleando la ecuación del calor escogiendo el ancho de banda con la regla de Scott modicada. El segundo problema que consideraremos será comparar la distribución espacial de los accidentes que ocurren en los dos tramos de hora punta. Extraemos para ello los accidentes que tienen lugar de 10 a 13 horas (45311 eventos) y de 20 a 23 horas (57159 eventos). Estos patrones de puntos se encuentran representados en la Figura 6.5 y en la Figura 6.6, respectivamente, junto con estimaciones no paramétricas de sus funciones de densidad asociadas. Nuevamente, para estos dos patrones de puntos estimamos el nivel crítico para los contrastes de Kolmogorov-Smirnov (empleando de nuevo como punto base el representado en la Figura 6.7) y de Cramer von Mises, empleando np= 10000 permutaciones. De nuevo, ambas estimaciones arrojaron un p-valor de 5·10−5 . Teniendo evidencias en contra de la hipótesis nula, podemos concluir que existen evidencias signicativas a favor de que la estructura espacial de los accidentes de tráco en las carreteras principales de Río de Janeiro cambia de la franja horaria de 10 a 13 a la franja horaria de 20 a 23.
57 Figura 6.7: representación del punto base del π -sistema del grafo de las carreteras principales de Río de Janeiro empleado para el contraste de Kolmogorov-Smirnov. Este punto se ha calculado empleando el código descrito para tal efecto en la Sección A.1.1. Los resultados obtenidos pueden resultar sorprendentes, ya que a primera vista no se observan grandes diferencias en las estimaciones de la densidad de las Figura 6.3 y la Figura 6.4, y tampoco entre las de la Figura 6.5 y la Figura 6.6. Ahora bien, las conclusiones obtenidas parecen tener sentido si pensamos, por ejemplo, que de lunes a viernes los accidentes de tráco ocurrirán en mayor medida en las zonas de mayor actividad laboral, mientras que los nes de semana los eventos se concentrarán en mayor medida en las zonas de ocio, cuyo emplazamiento es diferente en la ciudad de Río de Janeiro. Un comportamiento similar es de esperar también de la estructura espacial de los accidentes en la red de carreteras cuando se comparan las dos horas punta. Es importante notar que debido a los tamaños muestrales tan elevados que se manejan en este análisis, cualquier mínima variación entre los patrones puede ser detectada como signicativa y derivar en un rechazo de la hipótesis nula. Sería interesante, y será objeto de trabajo futuro, resolver problemas en rangos de tiempo más especícos que den lugar a comparaciones más razonables.
64 A. Código 1 library(spatstat) 2 NPtest=function(pp1,pp2,h,Heat=F,cont=F,ker="epa",CV=F){ 3 #Recordemos que el test no paramétrico calculará estimaciones de la función de riesgo relativo en los puntos de ambos patrones observados, y luego planteará un test de efecto en la regresión de estos valores estimados frente a las posiciones de los puntos en el grafo, para lo que se realizará una regresión no paramétyrica empleando una generalización a grafos lineales del estimador de Nadaraya-Watson. 4 #pp1 y pp2 son los patrones de puntos en el mismo grafo. 5 #h es el ancho de banda para la estimación no paramétrica de la función de intensidad y para la construcción de la matriz de suavizado. Se admite que se pase un vector con tres entradas, las dos primeras con el ancho de banda para la estimación de la densidad de cada proceso, y la tercera para el cómputo de la matriz de suavizado 6 #Heat es una variable binaria que nos indica si la estimación de la función de intensidad y los núcleos empleando en la regresión no paramétrica se hace (Heat=T) o no (Heat=F) empleando la estimación basada en la ecuación del calor. En caso de que Heat=F, se consideran núcleos equitativos. 7 #cont es una variable binaria. En el caso de que Heat=F, se emplearán núcleos equitativos discontinuos si cont=F y continuos si cont=T 8 #ker indica la función núcleo básica con la que se construyen los núcleos equitativos en el caso de que Heat=F. Por defecto se toma el núcleo de Epanechnikov. 9 #CV es una variable binaria que nos indica si la estimación de la función de riesgo relativo en los puntos observados se hace empleando validación cruzada (CV=T) o no (CV=F) 10 #notemos como por defecto la estimación de las funciones de intensidad y el cálculo de la matriz de suavizado se hace empleando núcleos equitativos discontinuos tomando como función núcleo básica el núcleo de Epanechnikov, y la estimación de la función de riesgo relativo se hace sin emplear validación cruzada. 11 #si no se especifica ancho de banda, este se toma para cada proceso puntual empleando la regla de Scott modificada, y para el caso de la matriz de suavizado se consideran los dos patrones superpuestos
A.1. Cálculo de los estadísticos 65 12 if (missing(h)){ 13 h=numeric(3) 14 h[1]=sqrt(sum(diag(cov(coords(as.ppp(pp1)))))) 15 h[1]=h[1]*(3*npoints(pp1))^(-1/5) 16 h[2]=sqrt(sum(diag(cov(coords(as.ppp(pp2)))))) 17 h[2]=h[2]*(3*npoints(pp2))^(-1/5) 18 #teniendo en cuenta que h[3] lo vamos a usar para estimar la intensidad empleando tanto puntos de 19 #ambos patrones, parece una elección natural usar el ancho de banda dado por la regla se scott superponiendo ambos patrones de puntos: 20 pp3=superimpose(pp1,pp2) 21 h[3]=sqrt(sum(diag(cov(coords(as.ppp(pp3)))))) 22 h[3]=h[3]*(3*npoints(pp3))^(-1/5) 23 } 24 # si se especifica un único ancho de banda, se entiende el mismo para todos los casos 25 if (length(h)==1){h=c(1,1,1)*h} 26 N1=npoints(pp1) 27 N2=npoints(pp2) 28 #recordemos que las estimaciones de las funciones de intensidad de los procesos puntuales no pueden anularse en ninguno de los puntos observados. Como no hemos impuesto ninguna condición sobre el ancho de banda para garantizar esto, lo que haremos será imponer un treshold a las estimaciones de la intensidad, tal que estas nunca sean menores que el valor de la estimación homogénea de la intensidad del proceso en cuestión. Estimamos la intensidad de forma homogénea como: 29 I1=N1/volume(domain(pp1)) 30 I2=N2/volume(domain(pp2)) 31 #calculamos entonces las estimaciones de la intensidad de cada proceso puntual en los puntos. Tenemos que distinguir en función de los valores de Heat y de CV 32 if (Heat==F){ 33 #estimamos las funciones de intensidad empleando la ecuación del calor. Para simplificar su evaluación, las tomamos como objetos linfun
66 A. Código 34 lambda1=as.linfun(densityEqualSplit(pp1, sigma=h[1], continuous = cont, kernel=ker, verbose=F, savehistory = F)) 35 lambda2=as.linfun(densityEqualSplit(pp2, sigma=h[2], continuous = cont, kernel=ker, verbose=F, savehistory = F)) 36 if (CV==F){ 37 #de no requerir validación cruzada evaluamos directamente 38 l1=c(lambda1(coords(pp1)),lambda1(coords(pp2))) 39 l2=c(lambda2(coords(pp1)),lambda2(coords(pp2))) 40 #en aquellos puntos en los que se estime una intensidad menor que el estimador homogéneo, sustituimos la estimacuón por la homogénea 41 l1[l1<I1]=I1 42 l2[l2<I2]=I2 43 #y estimamos entonces el logaritmo de la función de riesgo relativo 44 rho=log((N2/N1)*(l1/l2)) 45 }else{ 46 #de querer emplear validación cruzada, las estimaciones que tenemos nos valen para el otro patrón de puntos, pero para el respectivo tenemos que calcularlas empleando validación cruzada: 47 l11=densityEqualSplit(pp1, sigma=h[1], continuous = cont, kernel=ker, verbose=F, savehistory = F, at="points", leaveoneout = T) 48 attr(l11,"sigma")=NULL 49 l22=densityEqualSplit(pp2, sigma=h[2], continuous = cont, kernel=ker, verbose=F, savehistory = F, at="points", leaveoneout = T) 50 attr(l22,"sigma")=NULL 51 l12=lambda1(coords(pp2)) 52 l21=lambda2(coords(pp1)) 53 #nuevamente imponemos el treshold de la intensidad homogénea 54 l11[l11<I1]=I1 55 l12[l12<I1]=I1 56 l21[l21<I2]=I2 57 l22[l22<I2]=I2
A.1. Cálculo de los estadísticos 67 58 #y estimamos 59 rho=c(log((N2/(N1-1))*(l11/l21)),\n log(((N2-1)/N1)*(l12/l22))) 60 } 61 } else{ 62 #si Heat=F, todo es análogo cambiando la forma en la que se estima la intensidad: 63 lambda1=as.linfun(densityHeat.lpp(pp1,sigma=h[1],verbose=F)) 64 lambda2=as.linfun(densityHeat.lpp(pp2,sigma=h[2],verbose=F)) 65 if (CV==F){ 66 l1=c(lambda1(coords(pp1)),lambda1(coords(pp2))) 67 l2=c(lambda2(coords(pp1)),lambda2(coords(pp2))) 68 l1[l1<I1]=I1 69 l2[l2<I2]=I2 70 rho=log((N2/N1)*(l1/l2)) 71 }else{ 72 l11=densityHeat.lpp(pp1, sigma=h[1], verbose=F, at="points", leaveoneout = T) 73 attr(l11,"sigma")=NULL 74 l12=lambda1(coords(pp2)) 75 l21=lambda2(coords(pp1)) 76 l22=densityHeat.lpp(pp2, sigma=h[2], verbose=F, at="points", leaveoneout = T) 77 attr(l22,"sigma")=NULL 78 l11[l11<I1]=I1 79 l12[l12<I1]=I1 80 l21[l21<I2]=I2 81 l22[l22<I2]=I2 82 rho=c(log((N2/(N1-1))*(l11/l21)), 83 log(((N2-1)/N1)*(l12/l22))) 84 } 85 } 86 #una vez estimada la función de riegso relativo, llega el momento de efectuar la regresión. Bajo la hipótesis nula del test de efecto 87 mu=mean(rho) 88 #bajo la hipótesis alternativa del test de efecto, estimamos la función de regresión a través de la matriz de suavizado S 89 S=matrix(0,nrow=N1+N2,ncol=N1+N2)
68 A. Código 90 #ahora calculamos cada elemento, que es S_ij=G_{p_j}(p_i). Tenemos que distinguir como queremos calcular estos núcleos. 91 if (Heat==F){ 92 #empezamos con los puntos de pp1 93 for (j in 1:N1){ 94 #calculamos G_{p_j} como una linfun 95 Gj=as.linfun(densityEqualSplit(pp1[j],sigma=h[3], 96 continuous = cont,kernel = ker,verbose = F,savehistory = F)) 97 #y la evaluamos en los puntos de los patrones observados 98 S[,j]=c(Gj(coords(pp1)),Gj(coords(pp2))) 99 } 100 #y ahora hacemos lo mismo basando en los puntos de pp2 101 for (j in 1:N2){ 102 #calculamos G_{p_j} como una linfun 103 Gj=as.linfun(densityEqualSplit(pp2[j],sigma=h[3], 104 continuous = cont,kernel = ker,verbose = F,savehistory = F)) 105 S[,j+N1]=c(Gj(coords(pp1)),Gj(coords(pp2))) 106 } 107 }else{ 108 #prodecemos de fora análoga pero estimando mediante la ecuación del calor 109 for (j in 1:N1){ 110 Gj=as.linfun(densityHeat.lpp(pp1[j], sigma=h[3], verbose = F)) 111 S[,j]=c(Gj(coords(pp1)),Gj(coords(pp2))) 112 } 113 for (j in 1:N2){ 114 Gj=as.linfun(densityHeat.lpp(pp2[j], sigma=h[3], verbose = F)) 115 S[,j+N1]=c(Gj(coords(pp1)),Gj(coords(pp2))) 116 } 117 } 118 #ahora que ya tenemos la matriz S, debemos dividir cada elemento de S entre la suma de los elementos de su fila 119 S=sweep(S,1,apply(S,1,sum),"/")
A.2. Estimación de niveles críticos empleando el test de permutaciones 69 120 #una vez ya tenemos S podemos calcular los residuos cuadráticos de cada modelo de regresión bajo la hipótesis nula y alternativa. 121 #bajo la hipótesis nula 122 RSS0=sum((rho-mu)^2) 123 #bajo la hipótesis alternativa 124 RSSa=sum((rho-S%*%rho)^2) 125 #los grados de libertad bajo la hipótesis nula son 126 df0=N1+N2-1 127 #ya que se estima un único parámetro. Bajo la hipótesis alternativa, en analogía con el modelo lineal general 128 dfa=N1+N2-sum(diag(S)) 129 #Con todo esto ya podemos calcular el estadístico 130 Fhat=((RSS0-RSSa)/(df0-dfa))/(RSSa/dfa) 131 return(Fhat) 132 } 133 #EJEMPLOS 134 #empleando la ecuación del calor 135 NPtest(spiders[seq(1,48,2)], spiders[seq(2,48,2)], Heat = T) 136 #empleando la ecuación del calor, y estimando el log-riesgo relativo mediante validación cruzada 137 NPtest(spiders[seq(1,48,2)], spiders[seq(2,48,2)], Heat = T, CV = T) A.2. Estimación de niveles críticos empleando el test de permutaciones Habiendo especicado ya el código que permite calcular los tres estadísticos de contraste que hemos propuesto en el Capítulo 4, debemos especicar el código empleado para nuestro estudio de simulación. Primeramente debemos introducir las funciones que ejecutan el test de permutaciones para estimar el nivel crítico asociado al valor de un estadístico de contraste. A.2.1. Test de Kolmogorov-Smirnov Para el estadístico del test de Kolmogorov-Smirnov hemos empleado la siguiente función: 1 library(spatstat) 2 PermutationsKS=function(pp1,pp2,q,nboots){ 3 #pp1 y pp2 son los patrones de puntos observados
70 A. Código 4 #q es el punto base del pi-sistema que se desea emplear en el cálculo del estadístico de Kolmogorov-Smirnov 5 #nboots es el número de permutaciones que se van a considerar para estimar el estadístico. 6 t =KStest(pp1,pp2,q) 7 na = npoints(pp1) 8 nb = npoints(pp2) 9 n = nb + na 10 #queremos ahora combinar los dos patrones de puntos en uno solo 11 comb = superimpose(pp1,pp2) 12 #en el caso de que nboots no se haya pasado como entero 13 nboots = as.integer(nboots) 14 reps = bigger = 0L 15 boot_t=numeric(nboots) 16 for(idx in 1:nboots){ 17 #tomamos los subconjutnos de índices de forma aleatoria 18 e = sample.int(n, na, T) 19 f = sample.int(n, nb, T) 20 #calculamos el estadístico para las muestras permutadas 21 boot_t[idx] = KStest(comb[e], comb[f], q) 22 #si el valor del estadístico es mayor que el obtenido en las muestras originales, acumulamos 23 if (boot_t[idx] >= t){bigger = 1L + bigger } 24 } 25 #lo que vamos a devolver es el estadístico y el p-valor, que calculamos como la fracción de estadísticos calculados a partir de muestras permutadas mayores que el de partida 26 out = c(t, bigger/nboots) 27 #en el caso de que nuestro p-valor estimado sea cero, admitimos que este condicide con la mitad de la resolución 28 if (out[2] == 0) {out[2] = 1/(2 * nboots)} 29 #damos alrgo de formato 30 details = c(na, n - na, nboots) 31 names(details) = c("n1", "n2", "n.boots") 32 attributes(out) = list(details = details) 33 names(out) = c("Test Stat", "P-Value") 34 #finalmente le damos formato a la salida 35 out2=list("Test Stat"=out[1],"Pvalue"=out[2])
A.2. Estimación de niveles críticos empleando el test de permutaciones 71 36 return(out2) 37 } 38 #EJEMPLO 39 #EJEMPLO: 40 PermutationsKS(spiders[seq(1,48,2)], spiders[seq(2,48,2)], lpp(c(562.5,562.5), domain(spiders)), nboots = 5000L) A.2.2. Test de Cramer von Mises Para el test de Cramer von Mises y el test no paramétrico basado en la función de riesgo relativo las funciones que efectúan el test de permutaciones son análogas a la anterior, ya que únicamente hay que cambiar los parámetros de entrada y la llamada a las funciones que efectúan el cálculo del estadístico. Para el test de Cramer von Mises tenemos que: 1 library(spatstat) 2 PermutationsCvM=function(pp1,pp2,h,contin=F,kernel="epa",Heateq=F,nboots){ 3 #pp1 y pp2 son los patrones de puntos observados 4 #q es el punto base del pi-sistema que se desea emplear en el cálculo del estadístico de Kolmogorov-Smirnov 5 #nboots es el número de permutaciones que se van a considerar para estimar el estadístico. 6 t =CvMtest(pp1,pp2,h,cont=contin,ker=kernel,Heat=Heateq) 7 na = npoints(pp1) 8 nb = npoints(pp2) 9 n = nb + na 10 #queremos ahora combinar los dos patrones de puntos en uno solo 11 comb = superimpose(pp1,pp2) 12 #en el caso de que nboots no se haya pasado como entero 13 nboots = as.integer(nboots) 14 reps = bigger = 0L 15 boot_t=numeric(nboots) 16 for(idx in 1:nboots){ 17 #tomamos los subconjutnos de índices de forma aleatoria 18 e = sample.int(n, na, T) 19 f = sample.int(n, nb, T) 20 #calculamos el estadístico para las muestras permutadas 21 boot_t[idx] = CvMtest(comb[e], comb[f], h, cont=contin, ker=kernel,Heat=Heateq)
72 A. Código 22 #si el valor del estadístico es mayor que el obtenido en las muestras originales, acumulamos 23 if (boot_t[idx] >= t){bigger = 1L + bigger } 24 } 25 #lo que vamos a devolver es el estadístico y el p-valor, que calculamos como la fracción de estadísticos calculados a partir de muestras permutadas mayores que el de partida 26 out = c(t, bigger/nboots) 27 #en el caso de que nuestro p-valor estimado sea cero, admitimos que este condicide con la mitad de la resolución 28 if (out[2] == 0) {out[2] = 1/(2 * nboots)} 29 #damos alrgo de formato 30 details = c(na, n - na, nboots) 31 names(details) = c("n1", "n2", "n.boots") 32 attributes(out) = list(details = details) 33 names(out) = c("Test Stat", "P-Value") 34 #finalmente le damos formato a la salida 35 out2=list("Test Stat"=out[1],"Pvalue"=out[2]) 36 return(out2) 37 } 38 #EJEMPLO 39 PermutationsCvM(spiders[seq(1,48,2)], spiders[seq(2,48,2)], Heateq = T , nboots = 5000L) A.2.3. Test no paramétrico basado en la función de riesgo relativo Finalmente para el test no paramétrico basado en la función de riesgo relativo tenemos: 1 library(spatstat) 2 PermutationsNP=function(pp1,pp2,h,con=F,k="epa",Heq=F,cross=F,nboots){ 3 #pp1 y pp2 son los patrones de puntos observados 4 #q es el punto base del pi-sistema que se desea emplear en el cálculo del estadístico de Kolmogorov-Smirnov 5 #nboots es el número de permutaciones que se van a considerar para estimar el estadístico. 6 t =NPtest(pp1,pp2,h,cont=con,ker=k,Heat=Heq,CV=cross) 7 na = npoints(pp1) 8 nb = npoints(pp2)
A.2. Estimación de niveles críticos empleando el test de permutaciones 73 9 n = nb + na 10 #queremos ahora combinar los dos patrones de puntos en uno solo 11 comb = superimpose(pp1,pp2) 12 #en el caso de que nboots no se haya pasado como entero 13 nboots = as.integer(nboots) 14 reps = bigger = 0L 15 boot_t=numeric(nboots) 16 for(idx in 1:nboots){ 17 #tomamos los subconjutnos de índices de forma aleatoria 18 e = sample.int(n, na, T) 19 f = sample.int(n, nb, T) 20 #calculamos el estadístico para las muestras permutadas 21 boot_t[idx] = NPtest(comb[e], comb[f], h, cont=con, ker=k,Heat=Heq,CV=cross) 22 #si el valor del estadístico es mayor que el obtenido en las muestras originales, acumulamos 23 if (boot_t[idx] >= t){bigger = 1L + bigger } 24 } 25 #lo que vamos a devolver es el estadístico y el p-valor, que calculamos como la fracción de estadísticos calculados a partir de muestras permutadas mayores que el de partida 26 out = c(t, bigger/nboots) 27 #en el caso de que nuestro p-valor estimado sea cero, admitimos que este condicide con la mitad de la resolución 28 if (out[2] == 0) {out[2] = 1/(2 * nboots)} 29 #damos alrgo de formato 30 details = c(na, n - na, nboots) 31 names(details) = c("n1", "n2", "n.boots") 32 attributes(out) = list(details = details) 33 names(out) = c("Test Stat", "P-Value") 34 #finalmente le damos formato a la salida 35 out2=list("Test Stat"=out[1],"Pvalue"=out[2]) 36 return(out2) 37 } 38 #EJEMPLO 39 PermutationsNP(spiders[seq(1,48,2)], spiders[seq(2,48,2)], Heq = T , nboots = 5000L)
80 A. Código
Índice de notación 0 vector cero 1{·} función indicadora |A| área o longitud del conjunto A A conjunto de aristas de un grafo lineal AMISE error cuadrático medio integrado asintótico βmσ -álgebra de Borel en Rm #A cardinal del conjunto A CSR aleatoriedad espacial completa dG distancia del camino más corto en el grafo lineal G G grafo lineal E esperanza F función de distribución f función de densidad b f estimación de la función de densidad Gp,h función núcleo equitativo centrada en p con ancho de banda h ∇g vector gradiente de una función g que tome valores escalares deg(·) grado de un nodo de un grafo lineal h ancho de banda h AMISE ancho de banda óptimo respecto al AMISE h Norm ancho de banda óptimo respecto al AMISE para poblaciones normales H matriz de ancho de banda Hg matriz Hessiana de una función g que tome valores escalares Ha hipótesis alternativa 81
82 Índice de notación H0 hipótesis nula Id matriz identidad d -dimensional K función núcleo multidimensional k función núcleo unidimensional L densidad de probabilidad bivariante, isótropa, unimodal y centrada en 0 LG conjunto de puntos del plano del grafo lineal G l arista de un grafo lineal λ función de intensidad λc función de intensidad condicionada b λ estimación de la función de intensidad ln logaritmo en base e M número de réplicas Monte Carlo m marca de un proceso puntual MISE error cuadrático medio integrado MSE error cuadrático medio N(·) medida de contar Ni(·) medida de contar restringida al i -ésimo patrón de puntos n tamaño muestral determinista np número de permutaciones empleado en la estimación del nivel crítico N(µ,Σ) distribución normal de vector de medias µ y matriz de covarianzas Σ ||·|| norma euclídea Ω espacio muestral ω función núcleo básica P proceso puntual en un grafo lineal P probabilidad Pπ -sistema p patrón de puntos en un grafo lineal p punto de un grafo lineal p nivel crítico o p-valor P( A ) conjunto de partes del conjunto A Rm espacio euclídeo m -dimensional R+ conjunto de los números reales positivos
Índice de notación 83 RSS suma de residuos cuadráticos S matriz de suavizado S función núcleo apropiada para la regresión no paramétrica T estadístico TKS estadístico de Kolmogorov-Smirnov TCvM estadístico de Cramer von Mises TNP estadístico no paramétrico basado en la función de riesgo relativo V conjunto de nodos de un grafo lineal V ar varianza W región de observación de un proceso puntual en el plano X proceso puntual en el plano euclídeo X vector aleatorio X variable aleatoria x patrón de puntos en el plano euclídeo x punto del espacio euclídeo multidimensional x número real y patrón de puntos marcado en el plano Z(·) covariable a un proceso puntual
84 Índice de notación
Bibliografía [1] Bacchieri, G., y Barros, A. J. (2011). Trac accidents in Brazil from 1998 to 2010: many changes and few eects. Revista de Saúde Pública , 45, 949-963. [2] Baddeley, A., Nair, G., Rakshit, S., McSwiggan, G., y Davies, T. M. (2021). Analysing point patterns on networks A review. Spatial Statistics , 42. [3] Baddeley, A., Rubak, E. y Turner, R. (2015). Spatial Point Patterns: Methodology and Applications with R . Chapman & Hall. [4] Baddeley, A. y Turner, R. (2005). Spatstat: An R Package for Analyzing Spatial Point Patterns. Journal of Statistical Software , 12(6), 1-42. [5] Borrajo, M. I., González-Manteiga, W. y Martínez-Miranda, M. D. (2020). Testing for signicant dierences between two spatial patterns using covariates. Spatial Statistics , 40. [6] Bowman, A. W. y Azzalini, A. (1997). Applied smoothing techniques for data analysis: the kernel approach with S-Plus illustrations. Oxford Science Publications. [7] Bowman, A. W. (1984). An alternative method of cross-validation for the smoothing of density estimates. Biometrika , 71(2), 353-360. [8] Dantzig, G. B. y Thapa, M. N. (1997). Linear programming 1: introduction. Springer. [9] Dantzig, G. B. y Thapa, M. N. (2003). Linear programming 2: Theory and Extensions. Springer. [10] de Burgos, J. (2007). Cálculo Innitesimal de una variable (Segunda edición). McGrawHill. [11] DeGroot, M.H. (1988). Probabilidad y estadística (Segunda edición). Addison-Wesley Iberoamericana 85
86 Bibliografía [12] Diggle, P.J. (1985). A kernel method for smoothing point process data. Journal of the Royal Statistical Society (Series C) , 34, 138-147. [13] Diggle, P.J. (2013). Statistical Analysis of Spatial and Spatio-Temporal Point Patterns . Chapman & Hall. [14] Fuentes-Santos, I., González-Manteiga, W. y Mateu, J. (2021). Testing similarity between rst-order intensities of spatial point processes. A comparative study. Communications in Statistics-Simulation and Computation , 1-21. [15] Ghadirzadeh, M. R., Shojaei, A., Khademi, A., Khodadoost, M., Kandi, M., Alaeddini, F., y Moradi, S. (2015). Status and trend of deaths due to trac accidents from 2001 to 2010 in Iran. Iranian Journal of Epidemiology, 11(2), 13-22. [16] Gómez-Villegas, M. A. (2005). Inferencia estadística. Ediciones Díaz de Santos. [17] Google Maps ( https://www.google.es/maps/@-22.9140693,-43.5860658,11z?hl=es ). Consultado el 05/06/2022 [18] Gopalakrishnan, S. (2012). A public health perspective of road trac accidents. Journal of family medicine and primary care, 1(2), 144. [19] Marozzi, M. (2004). Some remarks about the number of permutations one should consider to perform a permutation test. Statistica , 64(1), 193-201. [20] McSwiggan, G., Baddeley, A. y Nair, G. (2017). Kernel density estimation on a linear network. Scandinavian Journal of Statistics , 44(2), 324-345. [21] Okabe, A., Satoh, T. y Sugihara, K. (2009). A kernel density estimation method for networks, its computational method and a GIS-based tool. International Journal of Geographical Information Science , 23(1), 7-32. [22] Okabe, A. y Sugihara, K. (2012). Spatial analysis along networks: statistical and computational methods. John Wiley & Sons. [23] Rakshit, S., Davies, T., Moradi, M. M., McSwiggan, G., Nair, G., Mateu, J. y Baddeley, A. (2019). Fast kernel smoothing of point patterns on a large network using two-dimensional convolution. International Statistical Review , 87(3), 531-556. [24] Rodríguez, G. (2003). Diferenciación de Funciones de Varias Variables Reales . Servizo de Publicacións e Intercambio Cientíco da USC. [25] Scott, D.W. (1992). Multivariate Density Estimation: Theory, Practice, and Visualization . John Wiley & Sons.
Bibliografía 87 [26] Venables, W. N. y Ripley, B. D. (2002) Modern Applied Statistics with S. (Cuarta edición). Springer. [27] Wand, M.P. y Jones, M.C. (1995). Kernel Smoothing . Chapman & Hall. [28] Yamada, I. y Thill, J. C. (2004). Comparison of planar and network K-functions in trac accident analysis. Journal of Transport Geography , 12(2), 149-158. [29] Yang, B. M., y Kim, J. (2003). Road trac accidents and policy interventions in Korea. Injury control and safety promotion, 10(1-2), 89-94. [30] Zhang, T. y Zhuang, R. (2017). Testing proportionality between the rst-order intensity functions of spatial point processes. Journal of Multivariate Analysis , 155, 78-82.