Herramientas estadísticas para el tratamiento de datos de estrellas binarias de la misión astrométrica GAIA

Traballo Fin de Grao Herramientas estadísticas para el tratamiento de datos de estrellas binarias de la misión astrométrica GAIA Ángel López Oriona 2018/2019 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA GRAO DE MATEMÁTICAS Traballo Fin de Grao Herramientas estadísticas para el tratamiento de datos de estrellas binarias de la misión astrométrica GAIA Ángel López Oriona Febrero, 2019 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA iii iv Trabajo propuesto Área de Conocimiento: Astronomía y Astrofísica - Estadística e Investigación Operativa. Título: Herramientas estadísticas para el tratamiento de datos de estrellas binarias de la misión astrométrica GAIA. Breve descripción del contenido La Agencia Espacial Europea lanzó en el año 2013 la misión astrométrica GAIA, que tiene por objetivo trazar el mapa tridimensional más preciso de nuestra galaxia, la Vía Láctea, con una muestra de hasta mil millones de estrellas. Entre ellas se encuentran las estrellas binarias, nuestro objeto de estudio, cuyos datos no serán publicados hasta 2021. En base al modelo simulado del Universo, elaborado durante la preparación de la propia misión GAIA, este trabajo pretende estudiar posibles técnicas estadísticas para analizar las relaciones entre los parámetros, tanto físicos como dinámicos, de las estrellas dobles, aplicables a futuras observaciones reales: regresiones lineales y no lineales, inferencia en modelos paramétricos de distribución y regresión, y tests de bondad de ajuste de los modelos considerados. Recomendaciones Haber cursado las materias de Fundamentos de Astronomía y Modelos de Regresión y Análisis Multivariante. Otras observaciones El trabajo consta de dos campos relativamente distintos, intentaremos mostrar la relación entre ellos de la mejor manera posible. Índice general Resumen viii Introducción xi 0.1. El Problema de los Dos Cuerpos . . . . . . . . . . . . . . . . . . . . . . . . xi 0.1.1. Preliminares................................ xi 0.1.2. ElProblema................................ xii 0.2. Coordenadas astronómicas y movimientos propios . . . . . . . . . . . . . . . xiii 0.2.1. Coordenadas ecuatoriales absolutas y coordenadas galácticas . . . . xiv 0.2.2. Movimientos propios de las estrellas . . . . . . . . . . . . . . . . . . xvi 0.3. ParalajeEstelar..................................xvii 0.4. EstrellasDobles..................................xix 0.4.1. Elementos orbitales un Sistema Binario . . . . . . . . . . . . . . . . xxii 0.5. Magnitud de una estrella . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxiii 1. La Misión Gaia 1 1.1. Generalidades................................... 1 1.2. Instrumentos y proceso de medición . . . . . . . . . . . . . . . . . . . . . . 3 1.2.1. Constitución del satélite. Scanning Law . . . . . . . . . . . . . . . . 3 1.2.2. Medición de errores astrométricos . . . . . . . . . . . . . . . . . . . . 4 1.3. Catálogos de Gaia y lenguaje de consultas . . . . . . . . . . . . . . . . . . . 5 1.3.1. El Gaia Archive. Catálogos . . . . . . . . . . . . . . . . . . . . . . . 5 1.3.2. LenguajeADQL ............................. 7 1.3.3. Cross-match................................ 7 1.4. Códigos de R y consultas ADQL . . . . . . . . . . . . . . . . . . . . . . . . 8 2. Completitud y contrastes 9 2.1. Introducción.................................... 9 2.2. Toma y filtrado de datos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 v vi ÍNDICE GENERAL 2.3. Análisis de completitud de GDR1 y GDR2 . . . . . . . . . . . . . . . . . . . 15 2.3.1. Análisis de completitud de GDR1 . . . . . . . . . . . . . . . . . . . . 15 2.3.2. Análisis de completitud de GDR2 . . . . . . . . . . . . . . . . . . . . 19 2.4. Contrastes de distribuciones en GDR2 . . . . . . . . . . . . . . . . . . . . . 21 2.4.1. Conceptos teóricos para los contrastes . . . . . . . . . . . . . . . . . 21 2.4.2. Datos para los contrastes de distribuciones . . . . . . . . . . . . . . . 27 2.4.3. Contraste de distribución para los paralajes . . . . . . . . . . . . . . 29 2.4.4. Contraste de distribución para los movimientos propios . . . . . . . . 36 3. Inferencia de distancias estelares 43 3.1. Generalidades sobre el tratamiento de los paralajes . . . . . . . . . . . . . . 43 3.2. El problema de la estimación de la distancia . . . . . . . . . . . . . . . . . . 45 3.3. El problema de los paralajes negativos en GDR2 . . . . . . . . . . . . . . . 53 3.4. Inferencia bayesiana de distancias estelares . . . . . . . . . . . . . . . . . . . 56 3.4.1. Planteamiento del problema . . . . . . . . . . . . . . . . . . . . . . . 56 3.4.2. PD Uniforme Impropia y Uniforme Propia . . . . . . . . . . . . . . . 58 3.4.3. PD de decaimiento exponencial de la densidad de volumen estelar . . 62 3.4.4. Estimaciones de distancias para estrellas dobles de GDR2C . . . . . 67 Bibliografía 71 Apéndice A 73 Apéndice B 75 xiv INTRODUCCIÓN Láctea. 0.2.1. Coordenadas ecuatoriales absolutas y coordenadas galácticas El objetivo de las coordenadas en Astronomía es el de identificar cada punto del firmamento. Un concepto clásico en este campo es el de esfera celeste (o bóveda celeste), definida como una esfera de radio arbitrario (normalmente la tomamos de radio unidad) concéntrica con la Tierra (que suponemos aquí una esfera, llamada esfera terrestre) en la cual están proyectados todos los astros. Gira en sentido contrario al de la Tierra (es lo que llamamos movimiento diurno aparente) respecto a un eje llamado eje del mundo, prolongación del eje terrestre. Sus polos norte y sur simplemente son la extensión a la misma de los polos norte y sur de la Tierra. Constituye un modelo muy intuitivo para poder situar cualquier astro. Sobre ella se definen las llamadas coordenadas astronómicas, de las que hay varios tipos, según qué planos se tomen como referencia. Todas las coordenadas astronómicas se definen tomando como base las coordenadas polares esféricas, que, como es bien conocido, son tres, una distancia y dos ángulos. Sin embargo, el hecho de que al mirar al cielo todos los astros parezcan situados a la misma distancia nos permite simplificar la situación, despreciando la coordenada que marca la distancia. Es por eso que cualquiera de las coordenadas astronómicas queda definida por dos ángulos. Nos vamos a centrar ahora en las llamadas coordenadas ecuatoriales absolutas, pues son uno de los dos tipos de coordenadas con los que trataremos en los análisis llevados a cabo en las secciones venideras. Estas coordenadas tienen la ventaja respecto a otras de no depender del lugar de la Tierra en el que está situado el observador ni del movimiento diurno. Para poder definir estas coordenadas necesitamos introducir brevemente una serie de conceptos. El ecuador celeste es una extensión del ecuador de la esfera terrestre (plano perpendicular al eje de rotación de la Tierra pasando por su centro); simplemente consiste en extender este plano a la esfera celeste. Como bien hemos visto en el Problema de los Dos Cuerpos, la Tierra describe una órbita plana (elíptica) alrededor del Sol, que se denomina eclíptica (línea en la que se producen los eclipses). El plano que contiene a dicha órbita se conoce como plano de la eclíptica. La intersección de dicho plano con el ecuador celeste se conoce como línea de los nodos. Además, estos dos planos forman un ángulo llamado oblicuidad de la eclíptica, que se denota por , y cuyo valor actual (varía a causa de los movimientos de precesión y nutación del eje de rotación terrestre) es aproximadamente 23o 260. Tanto el ecuador como el plano de la eclíptica como su intersección son independientes de cualquier lugar de la superficie terrestre en el que esté situado un observador. Como consecuencia del movimiento orbital de la Tierra podemos considerar un movimiento (el 0.2. COORDENADAS ASTRONÓMICAS Y MOVIMIENTOS PROPIOS xv que percibimos desde la Tierra) del Sol alrededor de la Tierra, cuya trayectoria en la esfera celeste será precisamente la eclíptica. La línea de los nodos interseca a la esfera celeste en dos puntos, los llamados punto aries ypunto libra. El punto aries es áquel por el que el Sol pasa el 21 de Marzo, en el llamado equinoccio de primavera (para habitantes del hemisferio norte). Este punto (también independiente del observador) va a ser clave en la definición de las coordenadas ecuatoriales absolutas. El sistema de referencia que usamos para definir las mismas va a ser el siguiente conjunto de vectores ortonormales con origen el centro de la Tierra: −→ v1, vector unitario en el dirección del punto aries, −→ v3, vector unitario en la dirección del polo norte y −→ v2=−→ v3×−→ v1, donde ×denota producto vectorial. Los dos ángulos que definen las coordenadas ecuatoriales absolutas de un astro son la ascensión recta α, que toma sus valores en [0,2π), y la declinación δ, que toma sus valores en [−π 2,π 2). Podemos ver la interpretación gráfica de estas coordenadas y de los conceptos anteriores en la Figura 2. Es habitual dar ambos valores en grados sexagesimales. La ascensión recta también se suele expresar en horas, minutos y segundos (360 grados se corresponden a 24 horas). De aquí en adelante, cada vez que hablemos de la posición de una estrella sin especificaciones mayores, supondremos que es la dada por las coordenadas ecuatoriales absolutas. El otro tipo de coordenadas que vamos a utilizar son las coordenadas galácticas. Éstas toman como referencia el plano de simetría de nuestra galaxia, denominado plano galáctico. La dirección perpendicular al plano galáctico se denomina polo galáctico. Determinar la ubicación de ambos no es una cuestión trivial, y se ha podido conseguir gracias a millones de observaciones astronómicas. El sistema de referencia galáctico, (−→ g1,−→ g2,−→ g3) es tal que que el plano galáctico contiene a a los vectores −→ g1y−→ g2, mientras −→ g3va en la dirección del polo norte galáctico; el vector −→ g1se elige en el sentido del centro galáctico. Las coordenadas galácticas son la latitud galáctica b, con valores en [−π 2,π 2)y la longitud galáctica l, con valores en [0,2π). Ambos sistemas de coordenadas, tanto las ecuatoriales absolutas como las galácticas, son sistemas de coordenadas dextrógiros, es decir, es el sentido de las agujas del reloj el que define a los giros positivos. Para una mejor comprensión de cualquier tipo de coordenadas astronómicas y conceptos relacionados, véase Abad et al. (2002) [1]. . xvi INTRODUCCIÓN Figura 2: Coordenadas ecuatoriales absolutas en la esfera celeste. 0.2.2. Movimientos propios de las estrellas Durante el curso de los siglos, las estrellas han aparentando mantener posiciones fijas unas con respecto a otras, formando siempre las mismas constelaciones. Sin embargo, actualmente se sabe que las constelaciones sí cambian su forma, pero tan lentamente que, en algunos casos, hacen falta miles de años para percibir esa diferencia. Esto es a consecuencia de que cada estrella tiene un movimiento intrínseco que llamamos movimiento propio. Este movimiento es interpretado como un movimiento relativo de las estrellas respecto al Sistema Solar. El Sol se mueve respecto al centro de la Vía Láctea a una velocidad aproximada de 220 km s−1, en una órbita cuasicircular, de un radio aproximado de 26000 años luz. El movimiento propio de una estrella (o de un astro) está caracterizado por dos cantidades, las variaciones angulares de la estrella en la ascensión recta αy en la declinación δ por unidad de tiempo. Si una estrella se mueve de la posición (α1, δ1)a la posición (α2, δ2) en un intervalo de tiempo ∆t, estas cantidades serán µα=α2−α1 ∆tyµδ=δ2−δ1 ∆t. Definimos el movimiento propio −→ µcomo un vector cuyo módulo es µ=µ2 δ+µ2 αcos2δ1. El factor cos2δ1tiene en cuenta el hecho de que la distancia lineal desde el eje del mundo a la esfera celeste varía con el factor cos δ1(será por ejemplo 0 en el polo norte celeste). Así, el movimiento propio será un vector velocidad (angular) en el plano, cuyas componentes serán 0.3. PARALAJE ESTELAR xvii Figura 3: Movimiento propio y sus componentes. Se han omitido las flechas en las magnitudes vectoriales. µαcos δ, usualmente denotada por µα∗, y µδ, llamadas movimiento propio en la ascensión recta ymovimiento propio en la declinación, respectivamente. Esta última componente paralela al ecuador tanto más pequeña cuanto más esté el objeto próximo a un polo celeste y constituye una medida de «cómo se aleja el astro de nosotros». La otra componente, en cambio, es una medida de «cómo cambia el objeto en nuestro plano de visión». Son estos dos valores los que suelen aparecer en las bases de datos estelares, usualmente medidos en mas año−1, donde mas denota milisegundos de arco. En el Apéndice Aencontramos una descripción completa de los acrónimos usados en el texto, basados casi todos ellos en el término inglés. La Figura 3 ilustra los conceptos previos. Una explicación más detallada acerca de los movimientos propios y de cómo se llega a la expresión de −→ µ se puede encontrar nuevamente en Abad et al. (2002) [1]. 0.3. Paralaje Estelar Otro concepto muy importante en el estudio de las estrellas y astros en general es el de paralaje anual (el término es indistintamente tanto masculino como femenino), que describimos a continuación. xviii INTRODUCCIÓN Figura 4: Ángulo de paralaje. Si consideramos la eclíptica, esta órbita es por supuesto elíptica, pero su excentricidad es tan próxima a cero que la podemos considerar circular sin cometer demasiado error. Dada una estrella Edel espacio, siempre podemos considerar el radiovector que une el Sol y la estrella E, cuya norma euclídea será la distancia del Sol a la estrella dada. Ahora, si consideramos una posición arbitraria T1de la Tierra sobre la eclíptica, y «miramos» a la estrella, la veríamos en cierta posición y respecto a un fondo de estrellas mas allá de la misma. Al pasar aproximadamente la mitad de un año, la Tierra estará en una posición T2, antipodal a T1, y si en este momento observamos la estrella, detectaríamos que el fondo de estrellas que vemos mas allá es distinto; para luego al transcurrir un año y llegar nuevamente a T1, volver a la misma situación del principio. Este movimiento aparente que describe la estrella debido al reflejo del movimiento orbital la Tierra define una elipse, llamada elipse de paralaje. Entre las posiciones T1yT2, este movimiento aparente dado por la elipse de paralaje define un cierto ángulo. Definimos el ángulo de paralaje, paralaje estelar o también paralaje anual como la mitad de dicho ángulo. De aquí en adelante, cada vez que hablemos de paralaje, nos referiremos al anterior, que es un tipo de paralaje trigonométrico. Existen otros tipo de paralajes en Astronomía, como los dinámicos y los espectroscópicos, en los que no entraremos. En la Figura 4 podemos ver que Pes el ángulo de paralaje. Este ángulo será menor cuanto más alejada esté la estrella del Sol. Al considerar la eclíptica como circular, la distancia de la Tierra al Sol no varía en el 0.4. ESTRELLAS DOBLES xix transcurso de la órbita, y esta distancia es por definición una unidad astronómica de distancia (au)(equivale aproximadamente a 1495978707 km). En la Figura 2 podemos ver claramente que tan P=au d, por lo que el paralaje está estrictamente relacionado con la distancia dde la estrella Ea la Tierra. Esto es muy importante, pues el conocimiento de los paralajes nos permitirá conocer distancias estelares, por lo que los paralajes serán un concepto que aparezca frecuentemente en los análisis estadísticos de datos astronómicos. Obviamente, para cualquier estrella, podemos reducir la situación a la de la Figura 4, por lo que a cada estrella le corresponde un y sólo un valor del paralaje. Un detalle importante es que la estrella más cercana al Sol (Próxima Centauri) tiene un valor del paralaje de P= 0.768700; este valor es, por tanto, el máximo valor de un posible paralaje estelar. Como estamos hablando entonces de ángulos próximos a cero, por truncamiento de series de Taylor podemos realizar la aproximación tan P=P, con lo que la expresión dada anteriormente resulta más sencilla, P=au d. Cuando el paralaje se da en segundos de arco (es lo habitual), se suele escribir P00. Se verifica que P00 =au dN00, donde N00 es el número de segundos de arco que tiene un radián (N00 =180 Π·60 ·60). Totalmente relacionado con ésto está el concepto de pársec (pc). Un pársec es la distancia a la que debería estar una estrella de la Tierra para que su paralaje fuese 1 segundo de arco. Como ya hemos visto que el paralaje de la más cercana es menor que tal paralaje, todas las estrellas estarán a más de 1 pc de nosotros. Se puede ver fácilmente que 1pc=N00 au. Otra equivalencia que puede resultar útil es 1 pc = 3.26 años luz. La relacíón importantísima que hay, pues, entre el paralaje de una estrella y su distancia en pc es P=1 d. Esta relación tiene importantes consecuencias estadísticas, que serán explicadas con detalle en el Capítulo 3. Debido a que es la notación utilizada en la mayoría de artículos y publicaciones astronómicas que utilizamos como referencia, de aquí en adelante denotaremos el paralaje anual de un astro por ω. 0.4. Estrellas Dobles En la actualidad se cree que la mayoría de las estrellas del Universo se encuentran asociadas en grupos de 2, 3, 4, etc , constituyendo lo que se llaman sistemas estelares múltiples. De éstos, los más abundantes son los de 2 componentes, que constituyen según la mayoría más del 80 % de los sistemas estelares, si bien esta cifra es siempre objeto de debate en Astronomía. xx INTRODUCCIÓN Una estrella doble (o estrella binaria) puede definirse como un par de estrellas físicamente ligadas por su mutua atracción gravitatoria, lo que, en virtud de El Problema de los Dos Cuerpos, tiene como consecuencia que cada una de ellas describa una órbita periódica respecto al centro de masas del sistema. Antes de introducirnos en aspectos puramente ténicos vamos a dar una breve reseña histórica de como estos sistemas estelares fueron descubiertos. Se le atribuye a J.B.Riccioli el haber descubierto la primera estrella doble. Fue en 1650 en el Observatorio de Palermo y se trata de la estrella Mizar, aunque no es descartable que ya Galileo se hubiese dado cuenta del caráter doble de este objeto. Para el conocimiento de los movimientos relativos de las estrellas dobles tenemos que irnos a 1718, cuando Edmund Halley detectó los movimientos propios de las estrellas y vio que no afectaban por igual a todas ellas, por lo que supuso que esto era a causa de que unas estrellas estaban a más distancia que otras de nosotros. Esto trajo como consecuencia intentar medir la distancia a las estrellas, descubriendo el movimiento paraláctico que acabamos de ver en la sección previa. Este problema del paralaje estelar preocupó a los astrónomos durante mucho tiempo, lo que hizo que el músico y astrónomo británico William Herschel (descubridor de Urano) también se involucrara en este campo. Con un telescopio reflector construído por el mismo, intentó observar la elipse de paralaje de ciertas estrellas; para su sorpresa, el período orbital de esta supuesta elipse de paralaje era mucho mayor que un año, que es lo que tardaría en completarse el movimiento paraláctico. Lo que en realidad estaba viendo Herschel eran los movimientos orbitales de las estrellas dobles. Las observaciones de Herschel eran la primera prueba de que la ley de gravitación de Newton era realmente universal. Para profundizar en la historia de las estrellas dobles, remitimos al lector a Coteau (2013) [7]. Hay multitud de clasificaciones para las estrellas dobles, pero quizá la más conocida sea la que divide a éstas en tres grandes grupos, a seguir: Estrellas dobles visuales: son pares de estrellas que pueden ser identificadas mediante métodos ópticos, ya sea usando telescopios, detectores electrónicos, etc. Habitualmente tomamos como origen de coordenadas la estrella más brillante. En consecuencia, la otra estrella describirá una órbita relativa (elipse) respecto a la primera. Sin embargo, esta órbita relativa no es la que nosotros observamos desde la Tierra, pues observamos la proyección de esta órbita sobre nuestro plano de visión, un plano perpendicular a nuestra visual (nuestra visual la marca un vector perpendicular a la esfera terrestre en el punto en el que nos encontramos), que es lo que denominamos órbita aparente. Haciendo un sencillo análisis geométrico se obtiene que la proyección de la elipse original sobre el plano 0.4. ESTRELLAS DOBLES xxi Figura 5: Ángulo de posición θy separación angular ρsobre la órbita aparente de una estrella doble. Nótese que la órbita no es en general una circunferencia. de la órbita aparente es otra elipse que conserva su centro, pero no en general sus focos. La forma que tenemos de medir esta órbita aparente es la siguiente. Denotamos por E la estrella que elegimos como principal y por E0la proyección de la secundaria, o estrella satélite. Damos la posición de la segunda respecto de la primera con un par de coordenadas polares sobre el plano de la órbita aparente. Tomamos como eje polar la dirección Norte y a partir de ahí medimos el ángulo de posición,θ, en sentido Norte-Este-Sur-Oeste; se mide en grados y decimal de grado. La segunda coordenada polar es ρ, la separación angular, y se trata del ángulo con vértice el observador según el cual se ven separadas las dos estrellas. Debido a que toma valores muy pequeños, es habitual dar esta coordenada en segundos de arco. Ilustramos ambas en la Figura 5. Estrellas binarias espectroscópicas: son estrellas binarias en las que sus dos componentes están tan próximas entre sí que no pueden ser resueltas por la vista, y en su mayoría ni siquiera usando poderosos telescopios. Existe un tipo de binarias espectroscópicas, las llamadas binarias espectro-interferométricas, para las que sí existen técnicas que permiten desdoblar su carácter doble. El carácter doble de las binarias espectroscópicas puede establecerse por el desplazamiento Doppler-Fizeau de sus líneas espectrales. Una breve explicación de en qué consiste es la siguiente: como ambas componentes del par estelar están describiendo una órbita respecto al centro de masas de sistema, cada una de ellas xxii INTRODUCCIÓN está alejándose y acercándose periódicamente al observador. En consecuencia, las líneas espectrales se desplazan en el espectro respecto a una posición que se corresponde a la del reposo relativo. La longitud de onda de una raya epectral, λ0, viene dada por λ0=λ(1+ Vr c); donde ces la velocidad de la luz, Vrla componente radial de la velocidad (en la dirección del observador) del vector velocidad del sistema con respecto al Sol y λla longitud de onda recibida o medida en el espectómetro. Por convenio, consideramos Vrpositiva cuando existe alejamiento (entonces λ0> λ), y negativa en caso de acercamiento (luego λ0< λ). En el primer caso las líneas espectrales se desplazan hacia el rojo, y en el segundo caso hacia el violeta. Dicho esto, si ambas estrellas del par no tienen una diferencia muy pronunciada en su luminosidad, podemos diferenciar en el espectro líneas correspondientes a las dos estrellas y, por tanto, cuando las de una componente se desplazan hacia el rojo, las de la otra lo hacen hacia el violeta. Estrellas binarias eclipsantes o fotométricas: son estrellas binarias cuyo plano orbital está orientado próximo al plano de visión del observador, de tal forma que desde nuestra perspectiva se producen eclipses entre ellas, que pueden parciales o totales. Su detección se basa en la búsqueda de patrones en las curvas de luz que recibimos de las mismas. En una binaria eclipsante, la curva de luz debe tener dos mínimos, el principal y el secundario. Aunque esta clasificación de las estrellas dobles es la dada de manera clásica, gracias a las modernas técnicas de visualización, podemos tener sistemas estelares que pertenezcan a los dos o incluso a los tres tipos. Hay un tipo de estrellas dobles, aquellas cuyas separaciones angulares son muy pequeñas, que se conocen como dobles cerradas (o con el término inglés close doubles). Sea cual sea el tipo al que pertenezcan, lo esencial es que cualquier sistema estelar doble encaja perfectamente en El Problema de los Dos Cuerpos. 0.4.1. Elementos orbitales un Sistema Binario Cuando intentamos estudiar un sistema estelar doble, se trata de determinar ciertos parámetros que denominamos elementos orbitales, cuyo cálculo nos permite obtener masas y distancias estelares. Vamos a suponer que estamos ante una binaria visual (pues si es de otro tipo varía la forma en la que tomamos los ejes del sistema de referencia, pero el procedimiento es análogo). Empecemos introduciendo un sistema de referencia espacial dextrógiro que llamaremos observable, (−→ s1,−→ s2,−→ s3), tal que el vector −→ s3está en la dirección de la visual, desde la estrella principal hacia el observador; el vector −→ s1es la dirección a 0.5. MAGNITUD DE UNA ESTRELLA xxiii partir de la cual se miden los ángulos de posición (dirección Norte) y −→ s2=−→ s3×−→ s1. En consecuencia, el plano −→ s1−→ s2es el que contiene la órbita aparente. El conjunto de elementos orbitales de un sistema binario es una 7-upla de valores (P, T, e, a, I, Ω, ω)medidos sobre la órbita relativa y definidos como sigue: P: período orbital (años). T: época de paso por el periastro (años). e: excentricidad de la elipse. a: semieje mayor de la elipse. I: inclinación. Es el ángulo diedro formado por los planos de la órbita relativa y aparente. Será I∈[0,90o)si el movimiento es directo y I∈(90o,180o]si el movimiento es retrógrado. Ω: ángulo del nodo. Es el ángulo formado por la dirección Norte con la línea de los nodos (intersección de los planos de las órbitas relativa y aparente). Definimos los nodos como la intersección de la línea de los nodos con la propias órbitas aparente y relativa. Por convenio se toma Ωen el intervalo [0,180o)mientras no se pueda precisar con medidas de velocidad radial. Se cuenta en sentido directo sobre el plano de la órbita aparente. ω: argumento del periastro. Es al ángulo medido sobre la órbita relativa y que va desde la posición del nodo hasta el periastro. Se cuenta en el sentido del movimiento. Es inmediato probar que al haber velocidad areolar constante en la órbita relativa, también la tenemos en la órbita aparente. 0.5. Magnitud de una estrella Para cuantificar el brillo de las estrellas se usa lo que denominamos magnitud. No entraremos en detalles sobre la definición de esta cantidad. Nos conformaremos con saber que su valor es tanto mayor cuanto menor es el brillo de la estrella. Además, la magnitud se puede medir en cierto rango de longitudes de onda, como puede ser el espectro visible, o en todo el rango espectral, definiendo lo que se llama magnitud bolométrica. Una amplia explicación sobre este parámetro, así como sobre todo lo relativo a estrellas dobles y al paralaje estelar, se puede encontrar en Abad et tal. 2002 [1]. 6CAPÍTULO 1. LA MISIÓN GAIA de dos catálogos estelares, el Tycho-Gaia Astrometric Solution (TGAS) y el Gaia Data Release 1 Catalogue (GDR1C). El primer catálogo contiene aproximadamente 2 millones de sources y tiene como objetivo identificar los objetos que ha detectado Gaia en sus primeras observaciones con dos catálogos previos, el Tycho-2 (basado en el satélite Hipparcos) y el propio Hipparcos. Para esto ofrece una columna con los identificadores de los sources en el catálogo Hipparcos y otra con los mismos en el catálogo Tycho-2. El segundo catálogo, GDR1C, contiene aproximadamente 1140 millones de sources. Los principales contenidos de GDR1 son: posiciones para los sources de los dos catálogos y cantidades astrométricas para los sources del TGAS. GDR2 se refiere a la segunda publicación de datos de Gaia. Corresponde a observaciones llevadas a cabo entre Julio de 2014 y Mayo de 2016. Publicado en Abril de 2018, consta de un catálogo, el Gaia Data Release 2 Catalogue (GDR2C). Éste contiene aproximadamente 1700 millones de sources, dando información acerca de sus magnitudes astrométricas, velocidades y otras muchas propiedas como la temperatura, el color, etc. En la Tabla 1.1 se detallan las cuestiones anteriores. Publicación Catálogo Número de sources Número de columnas GDR1 TGAS 2057050 59 GDR1C 1140622719 57 GDR2 GDR2C 1692919135 96 Tabla 1.1: Información relativa a los catálogos de Gaia. Ninguno de los tres catálogos anteriores es capaz por sí sólo de desdoblar estrellas dobles, en el sentido de que cada source debe ser tratado individualmente. Éste podría constituir, por tanto, una única estrella, un sistema estelar u otro objeto cualquiera. Si el satélite detecta un sistema estelar, lo traduce como un único source, y ésta es la única información que se obtendrá en estos tres primeros catálogos. Futuras publicaciones de datos resolverán esta cuestión. A su vez, el satélite necesita saber con cierta «seguridad» que un determinado source constituye un sistema estelar antes de trasnmitir la información. Ésto ha sido implementado en los algoritmos del DPAC. Por otra parte, los identificadores en GDR1C y GDR2C serán en principio diferentes, es decir, un determinado source que en GDR1C tiene un número de identificación, puede formar parte de GDR2C con otro número de identificación. Se puede obtener una descripción completa del proceso de identificación 1.3. CATÁLOGOS DE GAIA Y LENGUAJE DE CONSULTAS 7 de sources en Arenou el al. (2017) [2] para GDR1, y en Arenou et al. (2018) [3] para GDR2. Los catálogos estelares anteriores, y los que se publicarán en un futuro próximo, se pueden consultar en la página web de la ESA, concretamente en el Gaia Archive (GA) [20], sección dedicada exclusivamente a los datos recogidos por Gaia y a su procesamiento. 1.3.2. Lenguaje ADQL Debido a la enorme cantidad de datos almacenados en los catálogos estelares de las publicaciones GDR1 y GDR2, es impensable tener guardados todos los datos en un único archivo. La ESA utiliza una base de datos externa para el almacenamiento de los mismos. A través del GA, cualquier usuario puede realizar consultas y obtener los resultados relativos a los sources que cumplen una determinada condición y a las propiedades de éstos que más interesan. El lenguaje de consultas se denomina Astronomical Data Query Language (ADQL), y es una variante para Astronomía del famoso lenguaje de consultas Structured Query Language (SQL). Las tres sentencias principales de ambos son SELECT, donde se indican las variables que nos interesa seleccionar del catálogo estelar, FROM, donde se indica el catálogo estelar de donde queremos extraer los datos, y WHERE, donde se dan ciertas condiciones que queremos que cumplan nuestros datos de salida. El archivo de salida tendrá el mismo formato que el catálogo del cual ha sido extraído. Una explicación detallada del funcionamiento de este lenguaje de consultas viene dada en el mismo GA (apartado de ayuda). Ilustramos, a modo de ejemplo, una consulta concreta realizada en el GA, a seguir: SELECT parallax, ra, dec FROM gaiadr1.tgas_source WHERE parallax < 0.001 La consulta anterior pide los paralajes, ascensión recta y declinación de los sources del catálogo TGAS que verifican que su paralaje es menor que 0.001 mas. 1.3.3. Cross-match Aunque los catálogos de Gaia son en sí mismos una poderosa herramienta para la investigación astronómica, es su combinación con otros catálogos conocidos lo que verdaderamente permite exprimir todo su potencial. Aquí es donde surge el concepto de crossmatch. 8CAPÍTULO 1. LA MISIÓN GAIA Un cross-match (XM) entre dos catálogos estelares es un procedimiento que permite identificar sources de ambos catálogos. En líneas generales, los pasos a seguir para realizar un XM son los siguientes: 1) Se parte de un catálogo estelar, llamado catálogo estelar base, del cual se quieren identificar sus sources con otro catálogo estelar, llamado catálogo estelar de llegada. 2) Se define un algoritmo apropiado para asociar, en caso de que sea posible, a cada source del catálogo base un conjunto de posibles equivalentes en el catálogo de llegada. 3) Se define otro algoritmo para decidir, entre el conjunto de posibles candidatos, el que tiene más probabilidad de «ser» el source del que se parte. El procedimiento es muy complejo y requiere tomar ciertas hipótesis iniciales para desarollar el algoritmo. Existen distintos enfoques para realizarlo. El DPAC ha realizado diversos XM entre los catálogos estelares más famosos en el mundo de la Astronomía y los catálogos GDR1C y GDR2C; aunque todavía ninguno entre estos dos últimos. Los detalles de como han sido realizados los XM que involucran a GDR1C y a GDR2C pueden encontrarse, respectivamente, en Marrese et al. (2017) [16] y Marrese et al. (2019) [17]. El XM es, pues, un elemento clave cuando se trata de «cruzar» información entre varios catálogos estelares. Una situación habitual en la que se recurre a él es cuando tenemos un determinado conjunto de estrellas en un catálogo y necesitamos conocer acerca de ellas una propiedad sobre la que el catálogo dado no informa, pero sí lo hace otro. 1.4. Códigos de R y consultas ADQL En los análisis estadísticos sucesivos se usará el software libre R. Para permitir una lectura más cómoda, sólo vamos a incluir en los capítulos detalles relevantes del código de R que ha sido utilizado para obtener los resultados. Una puntualización que hacemos es que, cuando utilizemos conjuntamente datos relativos a la misma medida pero en distintos catálogos, redonderamos al número de cifras decimales del valor que menos cifras decimales posea. El código completo de R usado en cada sección, así como la correspondiente salida, se encuentra en el Apéndice B. De igual forma, son necesarias varias consultas ADQL y SQL con el objetivo de obtener los datos apropiados. Éstas también serán incluidas en el Apéndice B. Capítulo 2 Estudios de completitud y contrastes de distribuciones 2.1. Introducción Una manera sencilla de tener una idea de la completitud de un catálogo estelar recién publicado consiste en tomar otro catálogo estelar ya validado como auxiliar y estudiar que proporción de objetos de este catálogo aparecen en el primero. Al hablar de completitud nos referimos, pues, a la proporción de objectos que un catálogo ha sido capaz de detectar tomando cierto conjunto como referencia. Lo que pretendemos hacer ahora es, de alguna manera, medir la resolución angular que han tenido los telescopios del satélite Gaia en lo que respecta a las observaciones relativas a los dos primeros catálogos lanzados, el GDR1C y el GDR2C. Al hablar de resolución angular nos estamos restringiendo a las estrellas dobles y a su separación angular ρ, concepto que ya ha sido explicado previamente. No se trata, por tanto, de un análisis de completitud global, sino que queremos ver la capacidad de detección del satélite en función de la separación angular entre el par estelar. Aunque, como ya hemos dicho, no hay información acerca del carácter múltiple de los sources en estos dos primeros catálogos, consideraremos que una doble ha sido detectada en los mismos si existe en ellos un source que se identifica con una doble del catálogo auxiliar. La separación angular, por otra parte, es variable, pues como ya hemos visto estamos ante órbitas elípticas, pudiendo cambiar notablemente en función de la excentricidad. Sin embargo, es habitual en los catálogos de estrellas dobles dar únicamente un valor de ρ, la mayoría de los casos para una época determinada de observación que aparece especificada, sesgo que tendremos que asumir. ¿Detectarán las primeras observaciones de Gaia mejor las estrellas dobles que están 9 10 CAPÍTULO 2. COMPLETITUD Y CONTRASTES Figura 2.1: Tabla del catálogo DMSA. muy próximas?, ¿por el contrario, será el satélite capaz de distinguir mejor las estrellas que están muy separadas?, ¿o nos encontraremos con un nivel de detección similar en ambos casos?. A continuación intentamos dar respuesta a esta cuestión. 2.2. Toma de datos y filtrado del DMSA para los análisis de completitud Nuestro análisis va a ser llevado a cabo tomando como auxiliar el catálogo Double and Multiple System Annex (DMSA), producido por el satélite Hipparcos. El DMSA es un catálogo únicamente de sistemas estelares dobles y múltiples, como su nombre indica, y ha sido ampliamente estudiado. Este catálogo está disponible en forma de tabla en la base de datos astronómica VizieR, de libre acceso. Su formato es el siguiente: cada fila de la tabla se corresponde con una estrella perteneciente a cierto sistema estelar, y las estrellas del mismo sistema estelar se encuentran contiguas en la tabla. La tabla completa tiene 24588 filas. Cada columna de la tabla se corresponde, como en los catálogos de Gaia, con cierta propiedad de cada estrella. El número total de columnas de la tabla es 39. La base de datos permite al usuario descargar el catálogo seleccionando el número de filas deseado, así como las columnas que más interesen. Las estrellas que forman parte del mismo sistema estelar son aquellas que tienen el mismo identificador en la columna CCDM. Ilustramos una pequeña parte de esta tabla en la Figura 2.1. En ésta ya han sido elegidas las columnas que vamos a considerar al descargar el archivo completo, que son las que podemos necesitar para llevar a cabo nuestro análisis actual o posteriores. Las describimos a continuación: 2.2. TOMA Y FILTRADO DE DATOS 11 CCDM : Son las siglas del Catalog of Components of Double & Multiple stars publicado en Dommanget & Nys (2000) [8]. El número se corresponde con las coordenadas ecuatoriales absolutas aproximadas del sistema para el equinoccio 2000.0 en formato horas minutos y décimas de minutos +- grados y minutos (sexagesimales). QUAL: fiabilidad del sistema estelar. Es una variable que marca «en qué medida» podemos estar seguros de que el sistema estelar n-componente detectado constituye realmente tal sistema estelar. Toma 4 valores: A, B, C, D, correspondientes a las categorías «bueno», «aceptable», «pobre» e «incierto», respectivamente. Ncomp: número de componentes del sistema estelar. Nparm: número de parámetros libres. Se interpreta como el número de observaciones que el satélite ha realizado de un sistema estelar completo. Ncorr: número de registros de correlaciones. comp_id: identificación de componente en el sistema. Va tomando los valores A, B, C... con significado estrella principal, estrella secundaria, estrella terciaria... El catálogo toma como principal la estrella más brillante, como secundaria la segunda más brillante... y así sucesivamente. HIP: número de identificador en el catálogo Hipparcos. Hpmag: magnitud visual aparente, correspondiente a longitudes de onda entre 500 y 600 nm (mag). Es una medida que indica el grado de brillo en el rango visual que recibimos de una estrella en la Tierra. Se mide en unidades de magnitud visual aparente. RAICRS: ascensión recta para el equinoccio 1991.25 (grados sexagesimales). DEICRS: declinación para el equinoccio 1991.25 (grados sexagesimales). Plx: paralaje (mas). theta: ángulo de posición (grados sexagesimales). rho: separación angular (segundos de arco). _RA.icrs: ascensión recta para el equinoccio 2000.0 (grados sexagesimales). _DE.icrs: declinación para el equinoccio 2000.0 (grados sexagesimales). En la tabla de ejemplo ya se puede ver un primer problema con el que tendremos que enfrentarnos en el tratamiento de los datos, pues, como vemos, en el caso de sistemas estelares dobles, los valores del ángulo de posición y separación angular sólo aparecen en 12 CAPÍTULO 2. COMPLETITUD Y CONTRASTES Figura 2.2: Archivo de texto. la segunda fila de cada par. Lo primero que vamos a hacer es descargar desde VizieR el catálogo completo (todas las filas) con las columnas que hemos mencionado anteriormente. Lo descargaremos en forma de texto plano, de forma que luego lo podamos pasar fácilmente a un .txt y usar la función read.table de R, que permite cargar los datos desde un fichero .txt de una manera bastante amigable. Con unos comandos básicos del editor de texto, lo que hacemos es completar a cero los valores vacíos de las columnas theta yrho. Finalmente, tenemos un archivo de datos, llamado datos_dmsa.txt, fácilmente tratable con R y que presenta el aspecto de la Figura 2.2. Nos referiremos a este tipo de archivos, una vez que los hemos cargado con R, como data frames. Ahora, como hemos dicho, vamos a cargar los datos y a comprobar que la consola nos los ha cargado bien, comprobando que el data frame tenga 24588 filas y 15 columnas. A continuación, le ponemos nombre a las columnas con la función colnames. El comando head sirve para obtener una visión limpia de cuantas filas del data frame queramos, parámetro que le entra como argumento. > datos_dmsa<-read.table("datos_dmsa.txt", header=FALSE) > nrow(datos_dmsa) [1] 24588 > ncol(datos_dmsa) [1] 15 colnames(datos_dmsa)<-c("CCDM", "Qual", "Ncomp", "Nparm", "Ncorr", "comp_id", "HIP", 2.2. TOMA Y FILTRADO DE DATOS 13 "HPmag", "RA", "DE", "parallax", "theta", "rho", "RA2000", "DE2000") > head(datos_dmsa, 2) CCDM Qual Ncomp Nparm Ncorr comp_id HIP 1 00003-4417 A 2 11 1 A 25 2 00003-4417 A 2 11 1 B 25 HPmag RA DE parallax theta rho 1 6.894 0.07936537 -44.29030 13.74 0.0 0.000 2 7.551 0.07924029 -44.29021 13.74 315.8 0.463 RA2000 DE2000 1 0.07956353 -44.29056 2 0.07947489 -44.29047 Llegados a este punto, vamos a realizar la primera simplificación en nuestros datos. Como nuestro análisis va a estar centrado en estrellas dobles, necesitamos eliminar del data frame todas las filas correspondientes a estrellas pertenecientes a sistemas estelares de más de dos componentes. Para ello utilizamos la función subset, que como vemos tiene una sintaxis bastante intuitiva, pidiendóle a R que restringa nuestro data frame a aquellas filas donde Ncomp ≤2, obteniendo un nuevo data frame, datos_dmsa_dobles, cuyos datos ya estarán completamente restringidos a estrellas dobles. Vemos que apenas se han eliminado 578 filas respecto al data frame anterior, es decir, la proporción de estrellas del DMSA que pertencen a sistemas estelares de más de 2 componentes es de menos del 3 %. datos_dmsa_dobles<-subset(datos_dmsa, Ncomp<=2) > nrow(datos_dmsa_dobles) [1] 24010 Vamos a tener en cuenta ahora las recomendaciones que da la ESA y usadas por Arenou et al. (2017) [2] (sección 4.4.2) en cuanto a los análisis de completitud de catálogos estelares de estrellas dobles se refiere, con el fin de obtener un análisis lo más exacto posible. Tales son, descartar del catálogo auxiliar aquellos pares de estrellas dobles cuya magnitud de alguna de sus componentes sea mayor que una magnitud visual aparente de 20, así como aquellos pares cuya separación angular supere los 10 segundos de arco. Decidimos descartar tambíen aquellos sistemas estelares que no posean una etiqueta A o B en la columna Qual, 14 CAPÍTULO 2. COMPLETITUD Y CONTRASTES quedándonos así con los pares que son considerados como fiables de acuerdo con el catálogo auxiliar. Son tres, pues, los filtrados a llevar a cabo. Realizamos estos filtrados, en los que aplicamos entre otras cosas varios bucles y sentencias condicionales, y obtenemos el data frame datos_dmsa_dobles_filtrados2. Además, hemos adaptado este data frame para que cada fila se corresponda con una estrella doble, concretamente con la componente principal. Como hay solamente un valor de rho para cada sistema doble, esto no supone pérdida de información para nuestro análisis. Este último data frame presenta el aspecto siguiente. > head(datos_dmsa_dobles_filtrados2, 2) CCDM Qual Ncomp Nparm Ncorr comp_id HIP 2 00003-4417 A 2 11 1 B 25 4 00004-4711 A 2 9 1 B 37 HPmag RA DE parallax theta rho 2 7.551 0.07924029 -44.29021 13.74 315.8 0.463 4 11.745 0.10532213 -47.17955 3.74 332.0 0.230 RA2000 DE2000 2 0.07947489 -44.29047 4 0.10529738 -47.17953 Como el objetivo, como se ha dicho al principio del capítulo, es analizar la completitud en función de la separación angular, necesitamos agrupar de alguna manera los datos del data frame según pertenezcan a cierto intervalo de separaciones angulares ρ. Para esto solamente nos hacen falta dos columnas del data frame, el identificador de Hipparcos Hip yrho, por lo que generamos un nuevo data frame restringiendo el último a estas dos únicas columnas, llamado datos_dmsa_dobles_filtrados2_reducido, que es el que usaremos para llevar a cabo nuestro anáslisis. Habiendo descartado antes los valores de rho mayores que 1000, procedemos a hacer una partición del intervalo (0,10) (no hay ningún valor de separación angular que sea exactamente 1000, así que podemos considerar el intervalo abierto) formada por 50 intervalos con la misma longitud, con la que analizaremos la completitud. Guardamos los sistemas estelares de nuestro data frame en sus intervalos correspondientes en función de su valor de rho. La lista de 50 elementos part será tal que su componente j contendrá los valores de los Hip correspondientes a los pares estelares cuyo valor de rho cae en el intervalo jde la partición. Representamos finalmente un histograma del número de estrellas contenido en cada intervalo de la partición, que podemos ver en la Figura 2.3. El número mínimo de dobles contenidas en un intervalo es 24, lo que consideramos aceptable 2.3. ANÁLISIS DE COMPLETITUD DE GDR1 Y GDR2 15 para llevar a cabo un análisis de proporciones. Además, podemos apreciar claramente en el histograma que el número de estrellas a analizar en separaciones angulares pequeñas es mucho mayor que en separaciones angulares grandes. 2.3. Análisis de completitud de GDR1 y GDR2 2.3.1. Análisis de completitud de GDR1 Accedemos ahora al GA, concretamente a GDR1. Una vez en él, accedemos al TGAS, que es el catálogo que contiene los identificadores de Hipparcos. La única manera de cruzar información entre GDR1 e Hipparcos es a través del TGAS. Los sources de Hipparcos que han sido detectados en GDR1 aparecen en el TGAS. El identificador de Hipparcos de cada source del TGAS es el correspondiente a la columna hip. Llevamos a cabo 50 consultas ADQL, que consistirán en pedir que se nos devuelvan, para cada conjunto de identificadores de Hipparcos correspondientes a cada intervalo de la partición (cuyo número de elementos está almacenado en el vector totales, ordenado por intervalos), los que aparcen en el TGAS. El número de filas de cada archivo de salida en cada consulta será el número de dobles presentes en el TGAS (almacenado en el vector encontrados, también ordenado por intervalos), pues cada identificador de Hipparcos se ha asociado a una única fila del catálogo de GDR1 en la consulta ADQL. La completitud en cada intervalo simplemente la obtenemos realizando el cociente entre el número de identificadores HIP que ha entrado en la consulta ADQL y el número de filas del correspondiente archivo de salida. Una vez tenemos estos valores, estamos en condiciones de representar la completitud frente a la separación angular, que podemos ver en la Figura 2.4. Una presentación más detallada es la dada en la Tabla 3.3, donde se presentan los porcentajes de completitud en cada intervalo de separaciones angulares. Como podemos observar, el comportamiento general es que la completitud crece a medida que aumenta la separación angular, teniendo una caída considerable por debajo de separaciones de entre 1.5 a 2 segundos de arco. Esto concuerda con lo que dicen Makarov et al. (2017) [15], que fijan un límite de 1.5 segundos de arco para la detección de dobles en GDR1. Por encima de estos valores vemos que lo habitual es que GDR1 detecte como mínimo un 70% de dobles, llegando en muchos casos, para separaciones angulares por encima de 5 segundos de arco, a superar el 80 % . Hay un comportamiento algo anormal en el 22 CAPÍTULO 2. COMPLETITUD Y CONTRASTES Intervalo (00) Completitud ( %) Intervalo (00) Completitud ( %) [0, 0.2) 55% [5.0, 5.2) 62% [0.2, 0.4) 43% [5.2, 5.4) 72% [0.4, 0.6) 42% [5.4, 5.6) 63% [0.6, 0.8) 50% [5.6, 5.8) 68% [0.8, 1.0) 55% [5.8, 6.0) 66% [1.0, 1.2) 64% [6.0, 6.2) 75% [1.2, 1.4] 65% [6.2, 6.4) 78% [1.4, 1.6) 67% [6.4, 6.6) 70% [1.6, 1.8) 67% [6.6, 6.8) 74% [1.8, 2.0) 68% [6.8, 7.0) 51% [2.0, 2.2) 67% [7.0, 7.2) 59% [2.2, 2.4) 74% [7.2, 7.4) 68% [2.4, 2.6) 63% [7.4, 7.6) 69% [2.6, 2.8) 73% [7.6, 7.8) 70% [2.8, 3.0) 70% [7.8, 8.0) 74% [3.0, 3.2) 67% [8.0, 8.2) 63% [3.2, 3.4) 66% [8.2, 8.4) 63% [3.4, 3.6) 71% [8.4, 8.6) 70% [3.6, 3.8) 71% [8.6, 8.8) 58% [3.8, 4.0) 70% [8.8, 9.0) 51% [4.0, 4.2) 56% [9.0, 9.2) 72% [4.2, 4.4) 60% [9.2, 9.4) 63% [4.4, 4.6) 62% [9.4, 9.6) 58% [4.6, 4.8) 76% [9.6, 9.8) 63% [4.8, 5.0) 71% [9.8, 10) 52% Tabla 2.2: Completitud por intervalos de separación angular en GDR2. 2.4. CONTRASTES DE DISTRIBUCIONES EN GDR2 23 Figura 2.6: Completitud frente a separación angular en ambas publicaciones. Los triángulos azules corresponden a GDR2, mientras que los puntos rojos a GDR1. 24 CAPÍTULO 2. COMPLETITUD Y CONTRASTES paralaje es una variable aleatoria unidimensional. Comencemos por los paralajes, pues los conceptos serán más fáciles de entender en el caso univariante. Fijada una estrella doble, susceptible de que su paralaje sea medido tanto por Hipparcos como por Gaia (medida relativa a GDR2C), veamos lo que pasa con su diferencia de paralajes medidos. Las medidas del paralaje de esta estrella en ambos satélites son sendas variables aleatorias unidimensionales (de las que se tiene un valor muestral que es el valor que aparece en cada catálogo), que además suponemos independientes, pues la medición concreta del paralaje de un sistema estelar por Hipparcos no interfiere en la obtenida por Gaia. Una explicación más detalla de esta independencia puede verse en Marrese et al. (2019) [17]. Denotemos pues, como ωGa la variable aleatoria «medida del paralaje de la estrella relativa a GDR2C» y como ωHa su análoga para Hipparcos. Su diferencia será la variable aleatoria ∆ω=ωG−ωH. Suponemos que ambas variables están normalmente distribuidas con medias el paralaje verdadero, que denotamos por ωT, y respectivas varianzas σ2 ωGyσ2 ωH. Esta suposición es, como hemos visto en la sección 1.2.2, intrínseca al proceso de medida implementado en el satélite de Gaia, y lo mismo ocurre para Hipparcos; será de donde parta nuestro modelo. En Makarov et al. (2017) [15], Bailer-Jones (2015) [5], Astraatmadja & Bailer-Jones (2016b) [4] o Luri et al. (2018) [14] podemos ver como se hace uso de la misma. Denotamos lo anterior como ωG∼N(ωT, σ2 ωG)yωH∼N(ωT, σ2 ωH). Los valores de las desviaciones típicas σωGyσωHvienen dados en los catálogos estelares para cada source. Necesitamos recordar ahora ciertos resultados básicos del análisis multivariante. Los vectores serán en todo momento vectores columna. Teorema 2.1. Si Xes un vector aleatorio de dimensión m con distribución normal multivariante con vector de medias µy matriz de covarianzas Σ(denotamos X∼Nm(µ, Σ) ) yCes una matriz p×mde rango p, con p≤m, entonces: CX ∼Np(Cµ, CΣC0) Teorema 2.2. Si X1,X2. . . Xmson variables aleatorias normales univariantes y son mutuamente independientes, entonces el vector aleatorio (X1,X2,. . . ,Xm)0sigue una distribución normal multivariante. Si volvemos ahora a nuestras variables aleatorias, notemos que, por ser ωGyωHinde- 2.4. CONTRASTES DE DISTRIBUCIONES EN GDR2 25 pendientes, el vector (ωG, ωH)0será un vector aleatorio con distribución normal bivariante por el Teorema 2.2. Si tomamos ahora en el Teorema 2.1 la matriz C= (1 −1), obtenemos tras cálculos triviales que ∆ω=ωG−ωH∼N(0, σ2 ωG+σ2 ωH). Otro resultado importante que relaciona distribuciones es el siguiente. Teorema 2.3. Si X∼Nm(µ, Σ), entonces (X−µ)0Σ−1(X−µ)∼χm2, siendo χm2la distribución chi cuadrado con m grados de libertad. Dado la distribución de ∆ω, y el teorema anterior, tenemos que: (∆ω−0)(σ2 ωG+σ2 ωH)−1(∆ω−0) = (ωG−ωH)2 σ2 ωG+σ2 ωH∼χ2 1 El hecho de que esa variable aleatoria, que de aquí en adelante denotaremos por V, siga una distribución teórica tan sencilla es lo que nos va a permitir estudiar el comportamiento de los valores muestrales, y es aquí donde va a aparecer el concepto de Probability- Probability plot (P-P plot). Hemos visto que la variable aleatoria Vtiene esa distribución teórica y además se supone que tendremos «muchos» valores muestrales de esa variable. Es esperado, entonces, que los valores muestrales se comporten según esa distribución de probabilidad y lo que hace un P-P plot es analizar si realmente es así. De manera general, un P-P plot es un gráfico que representa dos funciones de distribución cualesquiera, una frente a la otra (recordemos que la función de distribución de una variable aleatoria unidimensional Xes la función que, para cada valor x, devuelve la probabilidad de que la variable aleatoria tome un valor menor o igual que x). Entonces, dadas dos funciones de distribución FyG, de ciertas variables aleatorias, el P-P plot representaría los puntos del plano de la forma (F(z), G(z)), con ztomando cualquier valor real. Por tanto, se trata de un gráfico parámetrico de dominio (−∞,∞)y de rango el cuadrado unidad [0,1] ×[0,1], producto de rangos de dos funciones de distribucion arbitrarias. El segmento que tenemos que tomar como referencia para la comparación es el que une los puntos (0,0) y(1,1), la diagonal del cuadrado unidad, pues es claro que las dos funciones de distribución son iguales si y sólo si el gráfico completo cae sobre la diagonal. Cualquier desviación indica diferencia entre las distribuciones. Siguiendo este razonamiento teórico, un P-P plot también puede ser utilizado para realizar comparaciones entre dos muestras (ver si ambas muestras proceden de una población con idéntica distribución), así como comparar una muestra frente a una distribución 26 CAPÍTULO 2. COMPLETITUD Y CONTRASTES teórica. Este tercer caso es el que nos ocupa. Hemos visto que la variable V=(ωG−ωH)2 σ2 ωG+σ2 ωH posee la distribución dada χ2 1. Si evaluamos los valores de la misma que cada estrella nos proporciona, tendríamos una muestra de esta variable V, susceptible de ser enfrentada en un P-P plot frente a la distribución teórica χ2 1, pudiendo extraer valiosas conclusiones según sea el comportamiento del gráfico. La idea es la siguiente: tenemos m valores muestrales de V:V1, . . . , Vm; ordenamos estos valores, obteniendo la muestra ordenada V(1), . . . , V(m), con V(1) ≤. . . ≤V(m). Definimos ahora la función de distribución muestral (también conocida como función de distribución empírica) Fm, que nos proporciona, para cada elemento de la muestra, la frecuencia relativa de los datos muestrales que son menores o iguales que éste. Es decir, será tal que Fm(V(i)) = i m. Si χes la función de distribución de χ2 1, los puntos a representar vienen dados por el conjunto {χ(V(i)), Fm(V(i))}, para i= 1, . . . , m. Si el conjunto de estos puntos cae «cerca» de la diagonal, aceptamos que la distribución de la muestra coincide con la teórica, mientras que si ocurre lo contrario, tendremos una discrepancia entre la distribución teórica y la muestral, estando obligados a analizar las causas de este desajuste, que pueden ser varias. Evidentemente, existen varios tests que miden «a partir de cuándo» podemos considerar que hay discrepancia, siempre fijando un nivel de significación, claro está. Aplicaremos algún test concreto más adelante. Las otras dos cantidades astrométricas que tenemos son los movimientos propios y las posiciones, que son vectores aleatorios bidimensionales. Escogemos los movimientos propios para hacer el análisis, siendo análogo para las posiciones. Razonando como en la situación previa, supongamos que tenemos sendos vectores independientes, uno representando la medida del movimiento propio de una estrella dada llevada a cabo por el satélite de Gaia y otra por Hipparcos, es decir, pg= (µα∗, µδ)0yph= (µα∗, µδ)0y consideramos el vector aleatorio, también bidimensional, definido por las diferencias entre ellos ∆p=pg−ph= (µα∗−µα∗, µδ−µδ)0. Supongamos ahora que pgyphson normales bidimensionales con el mismo vector de medias, el valor del movimiento propio real, digamos pT, y con matrices 2×2de covarianzas respectivas CgyCh. En nuestra notación, pg∼N2(pT, Cg)yph∼ N2(pT, Ch). Encontrar una distribución sobre la que apoyarse no es aquí tan fácil como en el caso unidimensional de los paralajes, y necesitamos recurrir a otro concepto básico en Estadística, que introducimos a continuación. Definición 2.4. Se llama función generadora de momentos de una variable aleatoria 2.4. CONTRASTES DE DISTRIBUCIONES EN GDR2 27 Xcon función de distribución Fa la función real MX(t) = E(etX ), (donde Edenota el operador esperanza), siempre que tal esperanza sea finita. Análogamente, si X es un vector aleatorio, se define su función generadora de momentos usando el producto escalar, como MX(t) = E(et0X). Ahora el domino será vectorial. Recordamos que la función generadora de momentos, en caso de existir, es única, y además caracteriza la distribución de probabilidad del vector aleatorio. Enunciamos otros dos resultados importantes relativos a la función generadora de momentos. Proposición 2.5. Si XeYson vectores aleatorios independientes, entonces MX+Y(t) = MX(t)MY(t) Proposición 2.6. Si Xes un vector aleatorio con distribución normal m−dimensional con vector de medias µy matriz de covarianzas Σ, entonces su función generadora de momentos es MX(t) = exp(t0µ+1 2t0Σt) Ahora el objetivo es aplicar las consideraciones anteriores al vector ∆p. Por definición, este vector es la diferencia de dos variables aleatorias normales bivariantes con los parámetros dados, que además son independientes. Teniendo en cuenta esto, la Proposición 2.5 y que −ph∼N2(−pT, Ch), la función generadora de momentos de ∆pserá M∆p(t) = Mpg(t)M−ph(t). Ahora, aplicando la Proposición 2.6, tenemos que M∆p= exp(t0pT+1 2t0Cgt) exp(−t0pT+1 2t0Cht) = exp(1 2t0(Cg+Ch)t), que es, precisamente, la función generadora de momentos de una distribución normal bidimensional, con vector de medias nulo y matriz de covarianzas Cg+Ch, que denotamos por C. Es decir, ∆p∼N2(−→ o , C). Finalmente, el Teorema 2.3 nos permite concluír que ∆0 p(Cg+Ch)−1∆p∼χ2 2, siendo χ2 2la distribución chi cuadrado con 2 grados de libertad. Denotamos por Wa esta variable aleatoria. Nuevamente, nos encontramos con una distribución sencilla con la que realizar un contraste con un P-P plot, de forma análoga a como lo hemos explicado antes. 2.4.2. Datos para los contrastes de distribuciones Los datos que vamos a usar son los correspondientes a las componentes de sistemas dobles extraídas del DMSA utilizadas en las secciones anteriores, después de haberles aplicado el tercer filtrado que se describe en la sección 2.2. El data frame es da- 28 CAPÍTULO 2. COMPLETITUD Y CONTRASTES tos_dmsa_dobles_filtrado_reducido. Para nuestro análisis vamos a necesitar ciertas variables que marcan la covarianza entre parámetros, y éstas no aparecen como variables en el catálogo DMSA, por lo que tendremos que descargar estos datos del archivo principal del catálogo Hipparcos y hacer la identificación utilizando los números HIP (recordemos que el DMSA es un suplemento del catálogo Hipparcos). Lo primero que hacemos es utilizar el XM de Hipparcos con el catálogo GDR2C e identificar cuáles de nuestras dobles del DMSA se encuentran en GDR2C. La consulta ADQL será la misma que las que hemos hecho antes diferenciando los identificadores de Hipparcos según a qué intervalo de separaciones angulares perteneciesen, sólo que ahora los introducimos todos en una única consulta. Pedimos como variables de salida los mismos identificadores de Hipparcos y los identificadores source_id de GDR2C. El objetivo es tener un data frame de dos columnas, que nos sirva para identificar los mismos pares en ambos catálogos, y que usaremos como «puente». Éste es gdr2_hipp. Necesitaremos en este momento los valores de la columna de identificadores en GDR2C, los source_id, para poder realizar una nueva consulta ADQL para conocer los valores de nuestras tres cantidades astrométricas, que se encuentran en el catálogo principal. Además, el propio GDR2C también nos permite conocer los valores de las desviaciones típicas de cada parámetro, y las correlaciones entre ellos. Recordemos que los parámetros astrómetricos que estamos estudiando son cinco: αyδ, que definen la posición; µα∗yµδ, que definen el movimiento proio, y el paralaje ω, siempre referidas a un source y a un satélite concre- to. Por tanto, podemos resumir toda la variabilidad astrométrica de un source en 15 cantidades: las varianzas de cada uno de estos parámetros y las 10 covarianzas existentes entre cada posible par. Muchas veces, como podemos ver en Makarov et al. (2017) [15], se da esta información en una matriz 5×5simétrica, que contiene todos estos valores. Para obtener la matriz de covarianzas relativa a los parámetros que nos interesen no hay más que tomar la submatriz adecuada. Tener toda la información astrométrica de un source relativa a un catálogo implica conocer los 5 parámetros que definen las magnitudes astrómetricas y los 15 que definen la variabilidad entre ellas, en ese catálogo; es decir, 20 cantidades. Habiendo realizado las consultas ADQL necesarias y tras varios procesos de filtrado, en los que eliminamos las estrellas de las que no se dispone del valor de ciertos parámetros, guardamos toda la información astrométrica de nuestras estrellas relativa a GDR2C y relativa a Hipparcos en un único data frame, que hemos tenido que obtener mediante un procedimiento muy cuidadoso. Este último contiene información relativa a 2045 pares 2.4. CONTRASTES DE DISTRIBUCIONES EN GDR2 29 estelares. 2.4.3. Contraste de distribución para los paralajes Para llevar a cabo este proceso vamos a usar el data frame parallax_data, extraído del total, y que contiene la información astrométrica relativa a los paralajes. Este data frame presenta la siguiente estructura: > head(parallax_data, 3) HIP Plx e_Plx source_id parallax parallax_error 1 40 -3.40 4.25 5.285634e+17 1.0788363 0.05140893 2 229 2.16 2.57 3.872488e+17 3.8659872 0.05280850 3 250 5.19 1.99 3.957315e+17 0.8343263 0.31346937 Para cada doble, Plx se corresponde con ωH,e_Plx se corresponde con σωH,parallax se corresponde con ωGypallax_error con σωG. Teniendo en cuenta eso y que lo que queremos contrastar es que los valores muestrales de Vproceden de una distribución χ2 1, el siguiente código nos dará un P-P plot relativo a los paralajes basándonos en los 2045 valores muestrales. Hemos respetado en la medida de la posible la notación usada en la sección 2.4.1. > omega_h <- parallax_data$Plx > omega_g <- round(parallax_data$parallax, 2) > sigma_h <- parallax_data$e_Plx > sigma_g <- round(parallax_data$parallax_error, 2) > var_h <- sigma_h^2 > var_g <- sigma_g^2 > vector_muestral <- (omega_g-omega_h)^2/(var_g+var_h) > vector_muestral1 <- sort(vector_muestral) > discreto <- seq(1, 2045) > acumulado <- discreto/2045 > puntos <- pchisq(vector_muestral1, df=1) > x <- seq(0, 1, length=10) > y <- x > plot(acumulado ~ puntos, type=’l’, main="P-Pplot para los paralajes" 30 CAPÍTULO 2. COMPLETITUD Y CONTRASTES Figura 2.7: P-P plot correspondiente a la variable V. + , xlab="Funcion de distribución teorica", + ylab="Funcion de distribucion muestral") > points(x,y, type=’l’, col=’blue’) La idea es la siguiente: 1) Calculamos los m= 2045 valores muestrales de Vy los ordenamos, generando el conjunto V(1), ..., V(m). 2.4. CONTRASTES DE DISTRIBUCIONES EN GDR2 31 2) Calculamos la función de distribución muestral en estos valores. Estos valores de la función de distribución muestral serán 1 m,2 m··· ,1. 3) Calculamos la función de distribución de χ2 1en cada uno de los valores muestrales ordenados, usando el comando pchisq, al que se le indican los grados de libertad. 4) Representamos los puntos del conjunto {χ(V(i)), Fm(V(i))}, para i= 1,...,2045. El P-P plot obtenido es el que se muestra en la Figura 2.7 A simple vista, parece que hay una desviación entre la distribución teórica y la distribución muestral, sobre todo en lo que respecta a la parte superior del gráfico. De todas formas, para «medir» esa desviación, vamos a utilizar un test de Kolmogorov-Smirnov (KS). Este test toma como hipótesis nula que la verdadera distribuciónde de la que han sido extraídos los datos es una χ2 1y como hipótesis alternativa que éstos provienen de cualquier otra distribución. La salida del test es lo que llamamos p-valor, que se puede interpretar como el mayor nivel de significación que nos permite aceptar la hipótesis nula. Es decir, un p-valor muy pequeño nos indica que debemos rechazar la hipótesis nula. En R, la sentencia para este test es ks.test, y tiene como argumentos de entrada el vector muestral a contrastar y la distribución respecto a la cual queremos hacer el contraste. Una justificación teórica sobre este test se puede encontrar en Vélez Ibarrola & García Pérez (1997) [19] (pág 454). > ks.test(vector_muestral1, "pchisq", 1) One-sample Kolmogorov-Smirnov test data: vector_muestral1 D = 0.042124, p-value = 0.00141 alternative hypothesis: two-sided La salida de R nos devuelve un p-valor de 0.00141. Esto nos dice que hay pruebas significativas para rechazar la hipótesis nula. Luego, la conclusión es que los datos no provienen de la distribución teórica dada. Tendremos que o bien buscar una causa para explicar esta desviación o bien cuestionar el modelo de partida. Vamos a intentar lo primero. Observamos en el gráfico que, para aproximadamente el 70 % de los valores ordenados 38 CAPÍTULO 2. COMPLETITUD Y CONTRASTES > head(proper_motion_data, 2) HIP source_id pmra pmdec pmra_error 1 40 5.285634e+17 -1.721576 -2.388676 0.07204935 2 229 3.872488e+17 18.683311 -4.237698 0.07839015 pmdec_error pmra_pmdec_corr pmRA pmDE e_pmRA 1 0.08050603 -0.2809894 -2.99 -3.18 4.14 2 0.05618416 -0.3010555 18.80 -5.10 1.59 e_pmDE pmDE.pmRA 1 3.75 -0.10 2 1.59 -0.08 > pmra_g_vector <- round(proper_motion_data$pmra, 2) > pmde_g_vector <- round(proper_motion_data$pmdec, 2) > pmra_sd_g_vector <- round(proper_motion_data$pmra_error, 2) > pmde_sd_g_vector <- round(proper_motion_data$pmdec_error, 2) > pmra_pmde_g_corr <- round(proper_motion_data$pmra_pmdec_corr, 2) > pmra_var_g_vector <- round(pmra_sd_g_vector^2, 2) > pmde_var_g_vector <- round(pmde_sd_g_vector^2, 2) > pmra_pmde_g_cov <- round(pmra_pmde_g_corr *pmra_sd_g_vector*pmde_sd_g_vector, 2) > pmra_h_vector <- proper_motion_data$pmRA > pmde_h_vector <- proper_motion_data$pmDE > pmra_sd_h_vector <- proper_motion_data$e_pmRA > pmde_sd_h_vector <- proper_motion_data$e_pmDE > pmra_pmde_h_corr <- proper_motion_data$pmDE.pmRA > pmra_var_h_vector <- pmra_sd_h_vector^2 > pmde_var_h_vector <- pmde_sd_h_vector^2 > pmra_pmde_h_cov <- pmra_pmde_h_corr*pmra_sd_h_vector*pmde_sd_h_vector > La obtención de los valores muestrales de la variable Wla realizamos mediante un bucle cuya implementación hace uso de los dos siguientes resultados: Lema 2.7. Si v= (a, b)0y M es la matriz simétrica de orden 2 y de entradas c11,c12 =c21 yc22, entonces la expresión v0Mv viene dada por a2c11 + 2abc12 +b2c22 Lema 2.8. Si Mes la matriz simétrica de orden 2 con entradas c11,c12 =c21 yc22 y 2.4. CONTRASTES DE DISTRIBUCIONES EN GDR2 39 det M6= 0, entonces su inversa M−1puede obtenerse como: M−1="c22 −c12 −c12 c11 #1 c11c22 −c2 12 El código de R que nos proporciona el P-P plot para Wes el siguiente: > muestra_p_m <- numeric(2045) > for (i in seq(1, 2045)) { + inv_det <- 1/(C[i,1]*C[i,2]-C[i, 3]^2) muestra_p_m[i]=delta_p[i,1]^2*inv_det*C[i,2]+delta_p[i, 1]*delta_p[i,2]*2*inv_det*(-C[i,3])+delta_p[i,2]^2*inv_det*C[i,1] + } > muestra_p_m_ordenada <- sort(muestra_p_m) > discreto1 <- seq(1, 2045) > acumulado1 <- discreto/2045 > puntos1 <- pchisq(muestra_p_m_ordenada, df=2) > x <- seq(0, 1, length=10) > y <- x > plot(acumulado1 ~ puntos1, type=’l’, + main="P-Pplot para los movimientos propios" + , xlab="Funcion de distribucion teorica", + ylab="Funcion de distribucion muestral") > points(x,y, type=’l’, col=’blue’) La Figura 2.11 nos muestra el P-P plot de los movimientos propios. Parece todavía más claro que en el caso anterior que hay una desviación respecto a la distribución teórica. Aplicamos, en todo caso, el test de KS. > ks.test(muestra_p_m_ordenada, "pchisq", 2) One-sample Kolmogorov-Smirnov test 40 CAPÍTULO 2. COMPLETITUD Y CONTRASTES Figura 2.11: P-P plot para la variable W. 2.4. CONTRASTES DE DISTRIBUCIONES EN GDR2 41 data: muestra_p_m_ordenada D = 0.10971, p-value < 2.2e-16 alternative hypothesis: two-sided El comportamiento de este gráfico parece ser parcialmente diferente al anterior. En primer lugar, no vemos un primer conjunto de puntos claramente por encima de la diagonal, como teníamos antes. Aunque en un principio la hipótesis de las varianzas y covarianzas sobrestimadas debería dar lugar a valores más pequeños de lo esperado según la forma cuadrática W, esto no se manifiesta prácticamente en este P-P plot. La única posible explicación que encontramos es que el efecto del fotocentro, de existir, sea en este caso (por cuestiones físicas que desconocemos) más influyente en la medida de los movimientos propios que en la de los paralajes y oculte la sobreestimación de las varianzas. Vamos a buscar indicios nuevamente de esta hipótesis. Necesitamos encontrar una relación similar a la que nos ofrecía el modelo lineal de antes, pero en este caso entre el cuadrado de la norma euclídea vector ∆p,||∆p||2(parámetro que haría el papel del numerador de antes, ahora en la forma cuadrática W) y el valor absoluto de la diferencia de magnitudes. De nuevo necesitamos aplicar un logaritmo para obtener resultados más claros. El diagrama de dispersión tras la transformación, así como la recta ajustada por un modelo de regresión lineal, se pueden ver en la Figura 2.12. De nuevo, vemos que el modelo lineal ajusta una recta con pendiente negativa, posible indicio del efecto de fotocentro. Además, el estimador de la pendiente en este caso es también significativo, y vale aproximadamente -0.19, menos que en el caso anterior, lo que podría interpretarse como que la influencia del efecto del fotocentro es mayor en el caso de los movimiento propios que en el de los paralajes. Esto explicaría el comportamiento de un mayor número de puntos por debajo de la diagonal, y de manera más pronunciada, en el segundo P-P plot con respecto al primero. Volvemos a remarcar que desconocemos las cuestiones por las que esto puede ser debido. En definitiva, parece que no hay motivos suficientes para cuestionar el modelo de partida ni en lo que respecta a los paralajes ni en lo que respecta a los movimiento propios. La hipótesis del fotocentro (parcialmente comprobada) podría explicar el comportamiento de ambos P-P plots, aunque de una forma preliminar y sin demasiado rigor. Teniendo en cuenta esto, podemos concluir entonces que el rendimiento astrométrico en GDR2 ha sido aceptable en lo que respecta a estrellas dobles. 42 CAPÍTULO 2. COMPLETITUD Y CONTRASTES Figura 2.12: Diagrama de dispersión del valor absoluto de la diferencia de magnitudes en el rango relativo a Hipparcos y ||∆p||2. Se presenta también la recta ajustada por un modelo de regresión lineal. Capítulo 3 Inferencia de distancias estelares a partir de los paralajes 3.1. Generalidades sobre el tratamiento de los paralajes Nuestro tema central en este capítulo van a ser principalmente los paralajes ω. De los mismos sabemos que están estrictamente relacionados con las distancias a las que se encuentran las estrellas de nosotros. Como bien sabemos, la relación teórica entre el paralaje y la distancia es que ambos son inversos el uno del otro. Esta cuestión, aunque parece banal, va a dar mucho juego en lo que respecta a los asuntos de inferencia. Supongamos que conocemos el valor del paralaje de una estrella dado en un catálogo estelar y nuestro objetivo es conocer la distancia a la estrella. Parece claro que no tenemos más que invertir el valor del paralaje, obteniendo el valor en las unidades correspondientes de la distancia buscada. Desafortunadamente, este enfoque es equivocado si lo que queremos es trabajar con paralajes medidos por un satétilite, debido a multitud de razones, como pueden ser las peculiaridades del propio aparato de medida o la incertidumbre de la medición. Cuando Gaia realiza mediciones astrométricas de los paralajes de los astros, lo hace basándose en la dirección que sigue el cuerpo en el firmamento, intentando modelar ésta como una función del tiempo; esto se realiza teniendo en cuenta tanto el propio movimiento del objeto en el espacio como el movimiento del propio satélite, en un proceso muy complejo. De modo resumido, sin entrar en detalles, presentamos aquí el modelo que utiliza Gaia para describir la dependencia temporal del movimiento de un objeto fuera del Sistema Solar, tomando como referencia la dirección hacia el observador. Ésta viene dada por el vector 43 44 CAPÍTULO 3. INFERENCIA DE DISTANCIAS ESTELARES unitario (donde norm denota normalización): −→ u(t) = norm(−→ r+ (tB−tep)(−→ p µα∗+−→ q µδ+−→ r µr)−ω−→ b(t) au )(3.1) En la expresión anterior, tes el tiempo de observación, tep es un tiempo de referencia, ambos medidos en unidades de Tiempo Baricéntrico Coordinado (TCB). Ésta es la escala de tiempo astronómico en el Sistema de Referencia Celeste Baricéntrico (BCRS), definido en el contexto de la relatividad general en la XXI Asamblea General de la Unión Astronómica Internacional. −→ p , −→ qy−→ rson los vectores unitarios en la dirección creciente de la ascensión recta, la dirección creciente de la declinación y hacia la posición del astro, respectivamente; tBes el tiempo de observación al que se le ha aplicado una corrección específica; −→ b(t)es la posición baricéntrica del observador (el satélite) en el tiempo de la observación (concepto relacionado con la posición respecto al baricentro del Sistema Solar, un punto que se puede interpretar como el centro de masas del sistema); au es la unidad astronómica de distancia. Las componentes del movimiento propio asociadas a −→ p y a −→ qson respectivamente µα∗yµδ,ωes el paralaje y µres el movimiento propio radial, que tiene en cuenta que la distancia al objeto cambia como consecuencia de su movimiento radial, que a su vez afecta a su movimiento propio y su paralaje. Este último término es normalmente despreciado sin afectar al comportamiento del vector resultante. El satélite obtiene los paralajes mediante un ajuste de este modelo a las observaciones. Las hipótesis de normalidad que han sido usadas para el paralaje en el capítulo previo están estrictamente relacionadas con este modelo. El modelo anterior predice un movimiento de patrón ondulatorio para el movimiento aparente de un astro dado. Podemos obtener una descripción completa del modelo y de su tratamiento en Luri et al. (2018) [14]. El ajuste del mismo a observaciones con mucha incertidumbre puede llevar a la obtención de paralajes sin sentido físico. El paralaje aparece en la ecuación 3.1 en el factor −ω au acompañando a la posición baricéntrica del observador, lo que significa que, para cada objeto, su movimiento paraláctico tendrá un sentido, que reflejará el sentido del movimiento del observador en torno al Sol. Si nuestras observaciones tienen alta cantidad de ruido, es decir, una alta variabilidad (lo que puede ocurrir muy fácilmente tratándose de observaciones astronómicas en un campo tan amplio como es la Vía Láctea), es completamente posible que el valor estimado del paralaje para un determinado astro sea nulo o negativo. Esto sería interpretado como el hecho de que la medida está siendo consistente con el movimiento del cuerpo «en la dirección incorrecta» sobre el firmamento. 3.2. EL PROBLEMA DE LA ESTIMACIÓN DE LA DISTANCIA 45 Dos importantes hechos pueden ser extraídos de lo que acabamos de exponer. Por un lado, aunque parezca contradictorio, los paralajes medidos pueden tomar valores negativos, y por otro, el paralaje en este contexto no es una medida directa de la distancia a un determinado objeto. En consecuencia, la distancia, y cualquier cantidad obtenida a través de la misma, debe ser estimada dado el paralaje observado (es decir, el obtenido por el satélite), teniendo en cuenta la incertidumbre o variabilidad en el proceso de medición. En conclusión, el tratamiento de los paralajes, lejos de lo que pudiese parecer, debe ser un proceso muy cuidadoso. 3.2. El problema de la estimación de la distancia Creemos conveniente introducir antes de continuar algunos conceptos básicos de la inferencia estadística, a seguir. Definición 3.1. Un estadístico (muestral) es una medida cuantitativa, derivada de un conjunto de datos de una muestra, con el objetivo de estimar o inferir características de una población. En otras palabras, un estadístico es una función de la muestra, que tiene como objetivo extraer información de una población. Definición 3.2. Un estimador ˆ θes cualquier estadístico que intente estimar un parámetro desconocido θde la población. Nótese que un estimador es una variable aleatoria. Definición 3.3. Para un estimador ˆ θde un parámetro poblacional θ, definimos el sesgo de ˆ θcomo: sesgo(ˆ θ)=E(ˆ θ)-θ. Definición 3.4. Decimos que un estimador ˆ θes insesgado si su sesgo es nulo. Definición 3.5. Dada una variable aleatoria X, se llama coeficiente de variación de X a la cantidad √V ar(X) E(X). El hecho de que un estimador sea insesgado quiere decir que su esperanza coincide con el parámetro que intenta estimar, lo que desde luego es una propiedad deseable para un estimador. Otra propiedad deseable para un estimador es que su varianza sea lo más baja posible, pues, a menor varianza, mejor estimará al parámetro que se quiere determinar. 46 CAPÍTULO 3. INFERENCIA DE DISTANCIAS ESTELARES Volviendo a nuestro contexto, denotaremos, a partir de ahora, para cada estrella, como ωTal valor del verdadero paralaje y como ωal valor del paralaje observado o medido por el satélite, supuesto estimador del parámetro ωT. Análogamente, denotaremos por r=1 ωT a la verdadera distancia a la que la estrella está situada, parámetro que queremos estimar, y como ρal valor 1 ω. Denotaremos por fTal coeficiente de variación de la variable que mide los paralajes, fT=σω ωT=σωr, y por fa su análoga muestral, f=σω ω. Interpretamos fTcomo una medida de la incertidumbre relativa (a la media). En la práctica, tanto fTcomo ωTson desconocidos. Los valores de los paralajes siempre se supondrán dados, mientras no se diga lo contrario, en segundos de arco y los valores de las distancias en pc. Las suposiciones de normalidad en la medida del paralaje observado nos permiten escribir su función de densidad gen función del paralaje verdadero ωTy la incertidumbre o desviación típica de la medición σω(σω≥0), como: g(ω) = 1 √2πσω exp(−(ω−ωT)2 2σω2) = 1 √2πσω exp(−(ω−1 r)2 2σω2)(3.2) Los procedimientos clásicos de inferencia en poblaciones normales llevan a que los siguientes intervalos de confianza para el parámetro 1 r, [ω−2σω, ω]y[ω, ω + 2σω](3.3) tengan, cada uno de ellos, una probabilidad de aproximadamente 0.477 de contener al parámetro 1 r. Este valor lo podemos obtener fácilmente con la función pnorm de R. La transformación de 1 rares monótona, y por lo tanto conserva probabilidades. Entonces, los intervalos, [1 ω,1 ω−2σω ]y[1 ω+ 2σω ,1 ω](3.4) tendrán también, cada uno, una probabilidad de aproximadamente 0.477 de contener a la verdadera distancia r. Sin embargo, mientras que los intervalos para 1 rson del mismo tamaño (la normal es simétrica), no lo son para r. Por ejemplo, para valores de ω=0.100 y σω=0.0200, estos intervalos son [10, 16.7] y [7.14, 10]. La razón está en que la transformación de 1 rares no lineal, y en consecuencia se pierde la simetría. Pero, ¿qué pasa si los errores de medición son «grandes»?, por ejemplo con f=1 2. Si nos fijamos en el intervalo de la izquierda en 3.4, sería [1 ω,∞). Este intervalo permite 3.2. EL PROBLEMA DE LA ESTIMACIÓN DE LA DISTANCIA 47 distancias tan grandes como queramos. Además, seguiríamos teniendo una probabilidad de 0.477 de que el verdadero valor de la distancia estuviese ahí dentro. Si el error fuese aún mayor, con f > 1 2, el intervalo no estaría definido, pues aparecería en el extremo superior un valor negativo. En este caso además, parece que «perdemos» alguna probabilidad, pues se deberían respetar los valores de 0.477 de probabilidad para cada intervalo. Vemos, por tanto, que la estimación de la distancia rpresenta una serie de problemas, sobre todo al aumentar el valor de f. Imaginemos por un momento que estamos ante una situación de ausencia de incertidumbre en la medida. Está claro que, conocido el paralaje ωT, (ωT=ωen este caso), obtendríamos trivialmente la distancia como r=1 ω. Parece que el enfoque más simple consistiría, pues, en estimar la distancia mediante la inversión del paralaje observado, es decir, mediante ρ. Tomar ρcomo estimador de la verdadera distancia r. Evidentemente, el uso de este estimador nos llevaría a distancias carentes de sentido físico en el caso de que el paralaje medido fuese un valor negativo (casos en los que el estimador de máxima verosimilitud sería el valor cero). Sin embargo, podríamos seguir considerando el uso de ρcomo estimador para valores positivos del paralaje medido, por ejemplo, en el caso de una muestra en la que la mayoría de los valores observados fuesen positivos, o incluso una muestra formada por un único valor positivo. Por tanto, nos interesa conocer las propiedades estadísticas del estimador ρ. Si tuviese «buenas propiedades», su uso limitado a valores positivos podría estar justificado. Para estudiarlas, necesitamos obtener la densidad de ρa partir de la densidad de ω.Teniendo cuenta que ρ=1 ωy aplicando un cambio de variable, tenemos que la densidad hde ρsigue la expresión: h(ρ) = g(ω)|dω dρ |=1 ρ2√2πσω exp(−(1 ρ−ωT)2 2σω2) = 1 ρ2√2πσω exp(−(1 ρ−1 r)2 2σω2),(3.5) que no es la densidad de una distribución normal. Veamos cómo se comporta esta función de densidad. Supongamos que tenemos dos estrellas situadas a distancias reales r1=50 pc yr2= 1000 pc, lo que corresponderían a valores de paralajes (reales) respectivos aproximados de ωT1=0.0200yωT2=0.00100. A su vez, supongamos que en ambos casos σω= 0.000300. Notemos que el coeficiente de variación fTes muy distinto en cada una de estas situaciones, siendo aproximadamente fT1=0.015 y fT2=0.30. Hemos elegido los valores de modo que las densidades otorguen una probabilidad muy pequeña a los valores negativos. Representemos ahora las densidades de ρen ambos 54 CAPÍTULO 3. INFERENCIA DE DISTANCIAS ESTELARES GDR2C e Hipparcos y obtenemos los valores ω. A continuación, imaginemos que estamos interesados en estudiar los sistemas estelares más lejanos de la muestra (recordemos que distancias grandes se corresponden teóricamente con paralajes pequeños), es decir, estamos interesados en los sistemas cuyo paralaje está por debajo de un determinado límite, fijado por nosotros. Fijamos este límite en 1.8 mas y nos quedamos con los datos correspondientes. A continuación, realizamos un truncamiento en los paralajes, quedándonos sólo con aquellos que son positivos, que aproximadamente corresponden al 88 % de la muestra de partida. En la Figura 3.4 podemos ver sendos histogramas, uno correspondiente a los datos iniciales y otro correspondiente a los datos truncados. La aparición de valores negativos con cierta frecuencia produce una notable variación en la distribución de los paralajes. La media de los paralajes iniciales vale aproximadamente 0.0006200, y la de los paralajes después de aplicar el truncamiento vale aproximadamente 0.000100, una diferencia claramente significativa. Por otro lado, desprendernos de los paralajes negativos implicaría rechazar parte de la información dada en el catálogo y renunciar a la estimación de la distancia a los respectivos sources. Para averiguar aproximadamente qué porcentaje de sources en GDR2C tienen un valor negativo del paralaje vamos a realizar cuatro consultas ADQL en distintas regiones de la galaxia (lo más alejadas posible) según la longitud galáctica l, y ver cuál es la proporción de paralajes negativos de cada una de las consultas. Esto es porqué realizar una única consulta ADQL pidiendo los sources con paralaje negativo en todo el catálogo no es factible debido al enorme volumen de datos; ni siquiera una consulta variando los valores de len un intervalo de longitud 1 grado sería posible. Los resultados de nuestras consultas se dan en la Tabla 3.1. A la vista de los datos, podemos estimar una proporción general de sources con paralajes negativos en GDR2C por encima del 15 %. No vemos como una opción factible desprendernos de casi una quinta parte del catálogo. El truncamiento de los datos queda entonces totalmente descartado. A modo de resumen, acabamos de ver que el estimador ρno es en general un estimador insesgado, tiene mucha varianza, y además no nos permite tratar con los valores negativos de los paralajes medidos, que representan un porcentaje significativo en los datos de GDR2C que no nos podemos permitir el lujo de descartar. La conclusión es clara, debemos buscar otro enfoque para la estimación de las distancias estelares. 3.3. EL PROBLEMA DE LOS PARALAJES NEGATIVOS EN GDR2 55 Intervalo de l(grados) Node sources Node sources con ωG<0Proporción (%) (0, 0.02) 313838 69737 22% (90, 90.02) 57588 10307 18% (180, 180.02) 25715 3937 15% (270, 270.02) 47398 9040 19% Tabla 3.1: Proporción de sources con valores negativos de los paralajes presentes en GDR2C, según valores de la longitud galáctica l. Figura 3.4: Izquierda: histograma de los paralajes antes del truncamiento. Derecha: histograma de los paralajes después de haber eliminado los valores negativos. 56 CAPÍTULO 3. INFERENCIA DE DISTANCIAS ESTELARES 3.4. Inferencia bayesiana de distancias estelares 3.4.1. Planteamiento del problema Antes hemos visto que el estimador ρno es adecuado para inferir distancias estelares, pero nuestro problema empezaba mucho antes. Hemos estado tratando de inferir el valor de rusando simplemente la ecuación 3.2, aún cuando ésta define la distribución de ω, no la de r. Una posible solución pasa por plantear el problema desde otro punto de vista, usando lo que se denomina Estadística Bayesiana. Vamos a empezar ilustrando este enfoque probabilístico con un ejemplo muy fácil de entender. Imaginemos que estamos pensando en la probabilidad de padecer determinada enfermedad a lo largo de nuestra vida, y la única información de la que disponemos es que la prevalencia global de esa enfermedad es del 0.01. En esta situación, podríamos decir, pues, que tenemos un probabilidad del 0.01 de padecer esa enfermedad a lo largo de nuestra vida, y eso es lo que hace el enfoque frecuentista. Sin embargo, consultando información acerca del trastorno, encontramos que la práctica de deporte previene esa enfermedad, y el consumo de alcohol y tabaco aumenta las posibilidades de padecerla. Está claro que, si somos deportistas y no fumadores, la probabilidad de que contrayamos esta enfermedad, aunque no la podamos dar explícitamente, es menor que el 0.01 de prevalencia general. En base a la observación, hemos conseguido un valor más «exacto» de la probabilidad buscada. La idea clave de la situación anterior es que, en vez de considerar la probabilidad de un determinado suceso como algo fijo e inamovible (como hace el enfoque frecuensita), podemos considerar la probabilidad como algo flexible, que va variando según aumenta la información de las que disponemos. Es decir, la probabilidad aquí representaría un grado de creencia en un determinado suceso, que puede cambiar. La herramienta fundamental en la que se basa este enfoque probabilístico es el Teorema de Bayes, aquí enunciado. Teorema 3.6. Sean AyBdos sucesos tales que P(B)6= 0, entonces se verifica P(A|B) = P(B|A)P(A) P(B), donde P(A|B)se interpreta como «la probabilidad de Asabiendo que Bes cierto», y análogamente P(B|A). En Estadística Bayesiana, normalmente Arepresenta una proposición (por ejemplo sobre un parámetro desconocido) y Blas «observaciones» relativas a A, como puede ser la información extraída de una muestra. El término genérico que se usa para referirse a B es el de datos. El término P(A)se conoce como probabilidad a priori, y se interpreta justamente de esa manera, es el conocimiento que tenemos de Asin tener en cuenta la 3.4. INFERENCIA BAYESIANA DE DISTANCIAS ESTELARES 57 evidencia ni ninguna nueva información. P(B)representa la evidencia, o nueva información que es tenida en cuenta. El término P(B|A)se denomina verosimilitud, y puede ser interpretada como la probabilidad de los datos sabiendo que Aes verdad; es una manera de cuantificar el grado en que la evidencia apoya a la proposición A. Finalmente, P(A|B)es la probabilidad a posteriori, es decir, la probabilidad de la proposición Adespués de tener en cuenta los datos. Dado que si consideramos variables continuas obtener el factor P(B) implica muchas veces difíciles cálculos integrales, muchas veces se obvia, ya que no varía, pues no depende de A, y se consideran simplemente los términos de la probabilidad a priori y la verosimilitud. Por último, notemos que, en el análisis que acabamos de hacer, si lo que queremos es inferir el valor de un cierto parámetro, dado el valor de otros, las probabilidades en el teorema de Bayes pueden ser interpretadas como funciones de distribución o de densidad. Así, nos referiremos a ellas como distribución a priori (PD) y distribución a posteriori (POD). Como hemos visto en el ejemplo, podemos interpretar el Teorema de Bayes como una forma de ir actualizando las probabilidades. Es decir, el término de la izquierda, la probabilidad buscada, es actualizado según vamos obteniendo nueva información, a través la verosimilitud y los datos. Volviendo a la cuestión de la estimación de la distancia, vamos a tratar de estimar nuevamente el valor de ra partir del valor del paralaje observado ω. Formalmente, lo que queremos conocer es la distribución de probabilidad de la distancia (término a posteriori) dado el valor observado de ωy la incertidumbre σω. Usando el teorema 3.6, el problema puede ser planteado de la siguiente manera. P(r|ω, σω) = 1 ZP(ω|r, σω)P(r)(3.7) El término verosimilitud P(ω|r, σω)es la distribución del paralaje observado dado el parámetro ry viene dado por la densidad de la ecuación 3.2. La PD P(r)contiene nuestros supuestos, y Z es lo que se conoce como una constante de normalización, que puede ser interpretado como P(ω), la distribución de la evidencia en el Teorema 3.6. En todo caso, no depende de ry es el término que permite convertir la POD en una densidad propia, de integral unidad. En este caso Z=Rr=∞ r=0 P(ω|r, σω)P(r)dr. Dos importantes elecciones tenemos que hacer al llevar a cabo la inferencia sobre r: la elección de la PD y la elección del estimador en la POD. Hay múltiples PD en las que podemos pensar. Un buena PD será aquella que concuerde tanto como sea posible con 58 CAPÍTULO 3. INFERENCIA DE DISTANCIAS ESTELARES los datos. En la práctica, debe contener toda la información relevante que tengamos y ser independiente de las medidas individuales. Ésta podría basarse en una combinación de la distribución esperada de las estrellas en la galaxia y de cómo éstas son seleccionadas por el satélite (resolución del satélite respecto a varias magnitudes, magnitud visual límite etc). Las posibilidades son infinitas, y probablemente lo mejor sea escoger una PD cada vez, según el problema concreto en el que queramos llevar a cabo la estimación, la zona de La galaxia donde queramos estimar las distancias y muchas más variables. A continuación se analizan tres casos de PD propuestos en Luri et al. (2018) [14], señalando las POD correspondientes, sus propiedades y sus limitaciones. Finalmente, se aplican los conceptos a la estimación de distancias de sources concretos de GDR2C. 3.4.2. PD Uniforme Impropia y Uniforme Propia Una primera PD en la que podríamos pensar sería una distribución uniforme no acotada para la distancia r, llamada distribución uniforme impropia (IUP). Esto tendría cierto sentido, pues estaríamos considerando que ren principio podría tomar cualquier valor sin favorecer a ninguno de ellos. Debido a esta razón, este tipo de PD se conoce en Estadística Bayesiana como distribución a priori no informativa. La distribución vendría dada este caso por: P∗ iu(r) =        1si r > 0 0si r ≤0 (3.8) Hemos introducido el símbolo * para indicar que la distribución define una densidad que no es normalizable (normalizar simplemente significa dividir entre el valor de la integral para conseguir un valor unitario al integrar), es decir, obtenemos un valor infinito al integrarla sobre r. Tales distribuciones son una extensión de las distribuciones de probabilidad finitas, y son conocidas como impropias (razón del nombre de esta distribución). De la ecuación 3.7 podemos deducir que la POD vendrá dada en este caso por la verosimilitud, pero considerándola ahora como función de ren vez de como función de ω, por lo que hemos de añadir de nuevo la restricción de ra valores positivos: P∗ iu(r|ω, σω) =        P(ω|r, σω)si r > 0 0si r ≤0 (3.9) 3.4. INFERENCIA BAYESIANA DE DISTANCIAS ESTELARES 59 Figura 3.5: POD correspondiente a IUP para diferentes valores de f. Podemos visualizar esta POD en la Figura 3.5 para un valor de ω= 0.0100 y varios valores de f. Vemos claramente que a medida que aumenta fse va formando una asimetría. Examinando la ecuación 3.2 se ve que, l´ım r→∞ P∗ iu(r|ω, σω) = cte En consecuencia, la POD no converge y su densidad encierra un área infinita, por lo que tampoco es normalizable (en realidad la definición de función de densidad implica que su integral valga uno, pero el lector notará usaremos el término también en estos casos). Por lo tanto, no tiene media, ni mediana, ni ningún otro cuantil. El único estimador razonable en 60 CAPÍTULO 3. INFERENCIA DE DISTANCIAS ESTELARES esta POD es la moda, que, como vemos en la Figura 3.5, está definida para cualquier valor de f, y coincide con el estimador 1 ωpara ω > 0. Para los valores ω≤0, si nos fijamos de nuevo en la ecuación 3.2, la POD crece desde r= 0, acercándose asintóticamente a cierto valor, con lo que la moda estaría en r=∞, careciendo de sentido físico. Por tanto, sólo podríamos tratar con paralajes positivos. Además, al estimar mediante la moda estaríamos usando de nuevo el estimador ρ. Usar como PD una IUP nos lleva, como vemos, al mismo resultado que ha sido analizado en la sección 3.2. Para lidiar con las limitaciones de la IUP sin dejar de utilizar una PD uniforme, lo ideal sería introducir un valor límite en las distancias, rlim (esto parece apropiado si tratamos de estimar distancias de estrellas en una «zona concreta» de la galaxia y disponemos de la información adecuada). La PD correspondiente sería una distribución uniforme propia (PUP): Pu(r) =        1 rlim si 0< r ≤rlim 0en otro caso (3.10) La POD correspondiente seguiría el mismo comportamiento que las de la Figura 3.5, pero anulándose para r > rlim, obteniendo: Pu(r|ω, σω) =        1 rlim P(ω|r, σω)si 0< r ≤rlim 0en otro caso (3.11) Esta POD no tiene primitivas elementales, sin embargo, al haber fijado un rlim encierra un área finita y podemos aproximar la integral numéricamente, pudiendo normalizarla. Por otro lado, tomar como PD una PUP nos permite tratar con valores ω≤0, pues al limitar el rango de valores de rno tenemos el problema de la asíntota que se presentaba en el caso de una IUP, alcanzándose el valor máximo en r=rlim. Varias POD normalizadas son presentadas en la Figura 3.6 para distintos valores de fy de ω. Se puede ver claramente como la POD es una combinación de la verosimilitud y la PD. Cuando «los datos son buenos», es decir, para valores pequeños de f, la verosimilitud domina la POD. En cambio, en el caso contrario, para valores grandes de f, el factor 1 rlim influye más en la POD, siendo ésta muy próxima a la anterior constante en todo su dominio. Como estimador podemos nuevamente tomar la moda de la POD, que estará definida de la siguiente manera: 3.4. INFERENCIA BAYESIANA DE DISTANCIAS ESTELARES 61 Figura 3.6: POD normalizada correspondiente a PUP para diferentes valores de f. El color negro se corresponde con un valor de ω=−0.0100. El resto, con un valor de ω= 0.0100. Se ha tomado un valor de rlim = 1000 pc. 62 CAPÍTULO 3. INFERENCIA DE DISTANCIAS ESTELARES Moda(P OD) =                  1 ωsi 0<1 ω≤rlim rlim si 1 ω> rlim rlim si ω ≤0 (3.12) El hecho de que las distancias asociadas a paralajes negativos sean estimadas mediante rlim es consistente con el proceso de medida de los paralajes, como podemos ver en Bailer- Jones (2015) [5]. Aún así, aunque estimar la distancia mediante la moda en esta situación presenta mayores ventajas que usando una IPUP, los resultados no son todavía satisfactorios, pues seguimos estimando mediante ρmuchas de las distancias. Por otra parte, ahora que tenemos una POD normalizada podemos pensar en estimar rtambién mediante la media o la mediana. Sin embargo, esta idea no parece demasiado factible, ya que, como vemos en la Figura 3.6, éstas estarán fuertemente influenciadas por la elección de rlim para valores grandes de f, que es justo donde necesitamos un estimador más robusto que la moda. Necesitamos buscar una PD que nos proporcione más y mejores opciones en las estimaciones. 3.4.3. PD de decaimiento exponencial de la densidad de volumen estelar Si queremos obtener una POD que nos proporcione estimaciones precisas para cualquier valor de fy de ω, podemos pensar en sustituír una PD caracterizada por una distancia límite como la del caso anterior, por una PD que decaiga ainstóticamente hacia 0 cuando r→ ∞. Vamos a investigar ahora una PD que supone una decaimiento exponencial en la densidad de volumen estelar (EDVDP) y que aparece definida en Luri et al. (2018) [14]. Se define de la siguiente manera, Pv(r) =        1 2L3r2exp(−r L)si r > 0 0si r ≤0 (3.13) donde Les una distancia. Éste parámetro está relacionado con el tamaño de la región galáctica en la que estemos considerando el decaimiento de la densidad volumétrica estelar. Una buena elección del mismo será crucial para la calidad de las estimaciones. Este tipo de distribución está encuadrada en las llamadas distribuciones gamma, muy importantes 3.4. INFERENCIA BAYESIANA DE DISTANCIAS ESTELARES 63 en Estadística. La definición de este tipo de distribuciones se puede consultar en Vélez- Ibarrola & García Pérez (1997) [19]. La POD correspondiente tendrá (salvo constantes), la siguiente forma: Pv(r|ω, σω) =        r2exp(−r L) σωexp( −1 2σ2 ω (ω−1 r)2)si r > 0 0si r ≤0 (3.14) Ofrecemos ejemplos de esta POD en la Figura 3.7. Dependiendo del valor de f, vemos que podemos tener 1 o 2 modas. En este ejemplo, para paralajes positivos, tenemos una única moda para 0 <f<0.30, dos modas para 0.30 ≤f < 0.373 y una moda para f≥0.373. A mayor valor de f, «menos información nos dan los datos» y por tanto, la POD se convierte en la PD. Podemos ver en el gráfico que para f= 0.5, la POD se acerca ya bastante a la PD. Para valores de fcercanos a 1, la POD será casi indistingible de la PD. La curva de color negro en la figura se corresponde a la POD para un paralaje negativo ω=−0.0100 y con un valor de |f|= 0.25. Si |f|fuese mayor, la curva se desplazaría hacia la derecha. Esto tiene sentido, ya que un valor de |f|más pequeño significa que «estamos más seguros» de que el verdadero valor del paralaje es cercano a cero. A medida que |f|aumenta, la curva se desplaza hacia la izquierda, pareciéndose cada vez más a la PD. Parece que obtenemos un comportamiento muy razonable en lo que respecta a la estimación de paralajes negativos. Para obtener analíticamente las modas, derivamos la ecuación 3.14 e igualamos a cero, obteniendo, r3 L−2r2+ω σ2 ω r−1 σ2 ω = 0 (3.15) que es una ecuación cúbica que tendrá en general 3 raíces complejas. Estas raíces serán función tanto de fcomo de L. Dado que ambos valores son positivos, una inspección de la ecuación nos permite concluír que existen dos casos posibles: que haya tres raíces reales, correspondientes a dos modas y un mínimo, o que haya una única raíz real, correspondiendo a una única moda. Esto concuerda perfectamente con lo que pasa en la Figura 3.7. Un análisis de las raíces nos permite llegar a la siguiente estrategia para estimar la distancia mediante la moda: 1) Si hay solamente una raíz real, es el máximo, así que ése será nuestro estimador, denotado rest. 70 CAPÍTULO 3. INFERENCIA DE DISTANCIAS ESTELARES > for (i in seq(1, 100000)) {areas[i]=integrate(densidad3, + lower=0, upper=prueba_error[i]) + + + } > indices <- which(areas > 0.049999) > indices[1] [1] 90966 > prueba_error[indices[1]] [1] 72.77273 > areas <- numeric(100000) > for (i in seq(1, 100000)) {areas[i]=integrate(densidad3, + lower=0, upper=prueba_error[i]) + + + } > indices <- which(areas > 0.949999) > prueba_error[indices[1]] [1] 73.65434 >1/omega4 [1] 119.19 Bibliografía [1] Abad, A., Docobo, J.A., Elipe, A. 2002, Curso de Astronomía, Prensas de la Universidad de Zaragoza [2] Arenou, F., Luri, X., Babusiaux, C., et al. 2017, Gaia Data Release 1: Catalogue validation, A&A, 599, A50 [3] Arenou, F., Luri, X., Babusiaux, C., et al. 2018, Gaia Data Release 2: Catalogue validation, A&A, 616, A17 [4] Astraatmadja, T.L., & Bailer-Jones, C.A.L. 2016b, Estimating distances from parallaxes II, ArXiv e-prints, [arXiv: 1609.03424] [5] Bailer-Jones, C.A.L. 2015, Estimating distances from parallaxes, Arxiv e-prints [ar- Xiv:1507.02105] [6] Bailer-Jones, C.A.L., Rybizki, et al. 2018, Estimating distances from parallaxes IV, ApJ, 156, 58 [7] Coteau, Paul. 2013, Estos astrónomos locos por el cielo o la historia de la observación de las estrellas dobles, editorial USC [8] Dommanget, J., Nys, O. 2000, The visual double stars observed by the Hipparcos satellite, A&A, v.363, p.991-994 [9] Gaia Collaboration., Prusti, T., de Bruijne, J. H. J., et al. 2016, The Gaia mission, A&A, 595, A1 [10] Galadí-Enríquez, D., & Ribas, I. 1999, Manual práctico de astrometría con CCD, editorial Omega [11] Holl, B., & Lindegren, L. 2012, Error characterization of the Gaia astrometric solution I, A&A, 543, A14 71 72 BIBLIOGRAFÍA [12] Holl, B., Lindegren, L., & Hobbs, D. 2012, Error characterization of the Gaia astrometric solution II, A&A, 543, A15 [13] Lindegren, L., Hernández, J., et al. 2018, Gaia Data Release 2: The astrometric solution, A&A, 616, A2 [14] Luri, X., Brown, A.G.A., et al. 2018, Gaia Data Release 2: Using Gaia parallaxes, A&A, 616, A9 [15] Makarov, V., Fabricius, C., & Frouard, J. 2017, Double Stars and Astrometric Uncertainties in Gaia Data Release 1, ApJ, 840, L1 [16] Marrese, P.M., Marinoni, S., et al. 2017, Gaia Data Release 1: Cross-match with external catalogues. Algorithm and results, A&A, 607, A105 [17] Marrese, P.M., Marinoni, S., et al. 2019, Gaia Data Release 2: Cross-match with external catalogues. Algorithm and results, A&A, 621, A144 [18] Smith, H. Jr., & Eichhorn, H. 1996, Statistical effect from Hipparcos Astrometry , Highlights of Astronomy, Vol.12 [19] Vélez Ibarrola, R., & García Pérez, A. 1997, Principios de inferencia estadística, editorial UNED [20] https://gea.esac.esa.int/archive/ (Consultado en Enero de 2019) Apéndice A Acrónimo Descripción A&A Astronomy & Astrophysics ADQL Astronomical Data Query Language ApJ Astrophysical Journal au Astronomical Unit DPAC Data Processing & Analysis Consortium EDVDP Exponentially Decreasing Volume Density Prior ESA European Space Agency GA Gaia Archive GDR1 Gaia Data Release 1 GDR2 Gaia Data Release 2 GDR1C Gaia Data Release 1 Catalogue GDR2C Gaia Data Release 2 Catalogue IUP Improper Uniform Prior KS Kolmogorov-Smirnov mas Milliarsecond PD Prior Distribution POD Posterior Distribution P-P plot Probability-Probability Plot PUP Proper Uniform Prior SQL Structured Query Language TCB Temps-coordonnée barycentrique BCRS Barycentric Celestial Reference System XM Cross-Match Tabla 3.3: Acrónimos 73 74 APÉNDICE A Apéndice B Código para 2.3 > datos_dmsa<-read.table("datos_dmsa.txt", header=FALSE) > nrow(datos_dmsa) [1] 24588 > ncol(datos_dmsa) [1] 15 > head(datos_dmsa, 4) V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 1 00003-4417 A 2 11 1 A 25 6.894 0.07936537 -44.29030 2 00003-4417 A 2 11 1 B 25 7.551 0.07924029 -44.29021 3 00004-4711 A 2 9 1 A 37 10.966 0.10536643 -47.17960 4 00004-4711 A 2 9 1 B 37 11.745 0.10532213 -47.17955 V11 V12 V13 V14 V15 1 13.74 0.0 0.000 0.07956353 -44.29056 2 13.74 315.8 0.463 0.07947489 -44.29047 3 3.74 0.0 0.000 0.10534168 -47.17959 4 3.74 332.0 0.230 0.10529738 -47.17953 > colnames(datos_dmsa)<-c("CCDM", "Qual", "Ncomp", "Nparm", + "Ncorr", "comp_id", "HIP", + "HPmag", "RA", "DE", "parallax", + "theta", "rho", "RA2000", "DE2000") > head(datos_dmsa, 2) CCDM Qual Ncomp Nparm Ncorr comp_id HIP HPmag 1 00003-4417 A 2 11 1 A 25 6.894 2 00003-4417 A 2 11 1 B 25 7.551 RA DE parallax theta rho RA2000 75 76 APÉNDICE B 1 0.07936537 -44.29030 13.74 0.0 0.000 0.07956353 2 0.07924029 -44.29021 13.74 315.8 0.463 0.07947489 DE2000 1 -44.29056 2 -44.29047 > length(parallax) [1] 24588 > datos_dmsa_dobles<-subset(datos_dmsa, Ncomp<=2) > head(datos_dmsa_dobles, 4) CCDM Qual Ncomp Nparm Ncorr comp_id HIP HPmag 1 00003-4417 A 2 11 1 A 25 6.894 2 00003-4417 A 2 11 1 B 25 7.551 3 00004-4711 A 2 9 1 A 37 10.966 4 00004-4711 A 2 9 1 B 37 11.745 RA DE parallax theta rho RA2000 1 0.07936537 -44.29030 13.74 0.0 0.000 0.07956353 2 0.07924029 -44.29021 13.74 315.8 0.463 0.07947489 3 0.10536643 -47.17960 3.74 0.0 0.000 0.10534168 4 0.10532213 -47.17955 3.74 332.0 0.230 0.10529738 DE2000 1 -44.29056 2 -44.29047 3 -47.17959 4 -47.17953 > nrow(datos_dmsa_dobles) [1] 24010 > datos_dmsa_dobles_fiables<-subset(datos_dmsa_dobles, Qual==’A’ | Qual==’B’) > nrow(datos_dmsa_dobles_fiables) [1] 20696 > HPmagdobles<-datos_dmsa_dobles_fiables$HPmag > length(HPmagdobles) [1] 20696 > head(datos_dmsa_dobles_fiables, 2) CCDM Qual Ncomp Nparm Ncorr comp_id HIP HPmag 77 1 00003-4417 A 2 11 1 A 25 6.894 2 00003-4417 A 2 11 1 B 25 7.551 RA DE parallax theta rho RA2000 1 0.07936537 -44.29030 13.74 0.0 0.000 0.07956353 2 0.07924029 -44.29021 13.74 315.8 0.463 0.07947489 DE2000 1 -44.29056 2 -44.29047 > auxiliar<-numeric(0) > vectorbucle<-seq(2, length(HPmagdobles), by=2) > for (i in vectorbucle) + { HPmagdoblesi<-HPmagdobles[i] + HPmagdoblesi_1<-HPmagdobles[i-1] + if (HPmagdoblesi<=20 & HPmagdoblesi_1<=20) + auxiliar<-c(auxiliar, i)} > datos_dmsa_dobles_filtrados1<-subset(datos_dmsa_dobles_fiables [auxiliar,]) > nrow(datos_dmsa_dobles_fiables) [1] 20696 > nrow(datos_dmsa_dobles_filtrados1) [1] 10348 > head(datos_dmsa_dobles_filtrados1, 2) CCDM Qual Ncomp Nparm Ncorr comp_id HIP HPmag 2 00003-4417 A 2 11 1 B 25 7.551 4 00004-4711 A 2 9 1 B 37 11.745 RA DE parallax theta rho RA2000 2 0.07924029 -44.29021 13.74 315.8 0.463 0.07947489 4 0.10532213 -47.17955 3.74 332.0 0.230 0.10529738 DE2000 2 -44.29047 4 -47.17953 > datos_dmsa_dobles_filtrados2<-subset(datos_dmsa_dobles_filtrados1, rho<10) > nrow(datos_dmsa_dobles_filtrados2) [1] 9345 > head(datos_dmsa_dobles_filtrados2, 4) 78 APÉNDICE B CCDM Qual Ncomp Nparm Ncorr comp_id HIP HPmag 2 00003-4417 A 2 11 1 B 25 7.551 4 00004-4711 A 2 9 1 B 37 11.745 6 00005+6713 A 2 9 1 B 40 11.176 8 00005-7212 A 2 9 1 B 45 11.954 RA DE parallax theta rho RA2000 2 0.07924029 -44.29021 13.74 315.8 0.463 0.07947489 4 0.10532213 -47.17955 3.74 332.0 0.230 0.10529738 6 0.11781651 67.21518 -3.40 224.9 8.200 0.11779774 8 0.13192459 -72.20307 15.10 242.5 2.830 0.13162877 DE2000 2 -44.29047 4 -47.17953 6 67.21517 8 -72.20308 > max(datos_dmsa_dobles_filtrados2$rho) [1] 9.99 > min(datos_dmsa_dobles_filtrados2$rho) [1] 0.09 > vec<-seq(0, 10, 0.2) > length(vec) [1] 51 > vec1<-vec[1:length(vec)-1] > vec2<-vec[2:length(vec)] > part<-numeric(length(vec)-1) > datos_dmsa_dobles_filtrados2_reducido<- subset(datos_dmsa_dobles_filtrados2, select=c(’HIP’, ’rho’)) > head(datos_dmsa_dobles_filtrados2_reducido) HIP rho 2 25 0.463 4 37 0.230 6 40 8.200 8 45 2.830 10 50 1.700 12 55 3.810 79 > all.equal(nrow(datos_dmsa_dobles_filtrados2_reducido), nrow(datos_dmsa_dobles_filtrados2)) [1] TRUE > for ( i in 1:(length(vec)-1)) + {part[i]<-subset(datos_dmsa_dobles_filtrados2_reducido , rho<=vec2[i] & rho>vec1[i])} > sum=0 > for (i in 1: (length(vec)-1)) + {sum=sum+length(part[[i]])} > print(sum) [1] 9345 > nrow(datos_dmsa_dobles_filtrados2_reducido) [1] 9345 > for (i in 1: (length(vec)-1)) + {write.table(part[[i]], paste0(i, ".txt"))} 50 búsquedas ADQL para GDR1C SELECT * FROM gaiadr1.tgas_source WHERE HIP in (vector de HIP) > totales<-numeric(50) > for (i in 1:(length(vec)-1)) + {totales[i]=length(part[[i]])} 86 APÉNDICE B FROM gaiadr2.hipparcos2_best_neighbour WHERE original_ext_source_id in (valores de identificadores de hipparcos) ORDER BY original_ext_source_id > gdr2_hipp <- read.csv("../Archivos TFG/2.4/gdr2_hipp.csv" + , header=T) > gdr2_vector <- write.table(gdr2_hipp$source_id, + file="../Archivos TFG/2.4/gdr2_vector.txt") Segunda consulta ADQL SELECT HIP, "20 cantidades astrometricas" FROM gdr2.gaia_source WHERE source_id in ("gdr2_vector") > gdr2 <- read.csv("../Archivos TFG/2.4/gdr2.csv", header=TRUE) > gdr2_parallax <- subset(gdr2, select=c("source_id","parallax", + "parallax_error" )) > head(gdr2_parallax, 5) source_id parallax parallax_error 1 3.245105e+15 10.517934 0.05829324 2 7.001557e+15 16.624334 0.28851213 3 7.264203e+15 6.048140 0.11639701 4 7.633158e+15 4.236715 0.07368882 5 7.870579e+15 4.842846 0.04417935 > nrow(gdr2_parallax) 87 [1] 2359 > gdr2_parallax <- subset(gdr2_parallax, gdr2_parallax$parallax!="NA" + & gdr2_parallax$parallax_error!="NA") > nrow(gdr2_parallax) [1] 2045 > matching <- gdr2_hipp$source_id %in% gdr2_parallax$source_id > head(matching) [1] FALSE TRUE FALSE FALSE FALSE FALSE > matching1 <- which(matching==TRUE) > length(matching1) # para excluir los valores NA [1] 2045 > previo_hipparcos <- subset(gdr2_hipp[matching1,]) > nrow(previo_hipparcos) [1] 2045 > hipparcos_parallaxes <- previo_hipparcos$original_ext_source_id > hipparcos_parallaxes_ordenados <- sort(hipparcos_parallaxes) > length(hipparcos_parallaxes_ordenados) [1] 2045 > write.table(hipparcos_parallaxes_ordenados, + file="../Archivos TFG/2.4/hipparcos_parallaxes_ordenados.txt") > hipparcos_astrometric_parallax <- +read.csv(file="../Archivos TFG/2.4/hipparcos_astrometric_parallax.csv", + header=TRUE) > nrow(hipparcos_astrometric_parallax) [1] 2045 > hipparcos_astrometric_parallaxes <- na.exclude(hipparcos_astrometric_parallax) > nrow(hipparcos_astrometric_parallax) [1] 2045 > previo_hipparcos_1 <- previo_hipparcos[with(previo_hipparcos, + order(previo_hipparcos$source_id)),] > gdr2_parallax1 <- gdr2_parallax[with(gdr2_parallax, + order(gdr2_parallax$source_id)),] > head(previo_hipparcos, 3) 88 APÉNDICE B source_id original_ext_source_id 2 5.285634e+17 40 8 3.872488e+17 229 9 3.957315e+17 250 > data_auxiliar <- merge(previo_hipparcos_1, gdr2_parallax1 ) > head(data_auxiliar, 2) source_id original_ext_source_id parallax parallax_error 1 3.245105e+15 15253 10.51793 0.05829324 2 7.001557e+15 14075 16.62433 0.28851213 > colnames(data_auxiliar) <- c("source_id", "HIP", + "parallax", "parallax_error") > data_auxiliar1 <- data_auxiliar[with(data_auxiliar, + order(data_auxiliar$HIP)),] > hipparcos_astrometric_parallax1 <- + hipparcos_astrometric_parallax[with(hipparcos_astrometric_parallax, + order(hipparcos_astrometric_parallax$HIP)),] > parallax_data <- merge(hipparcos_astrometric_parallax1, data_auxiliar1) > head(parallax_data, 4) HIP Plx e_Plx source_id parallax parallax_error 1 40 -3.40 4.25 5.285634e+17 1.0788363 0.05140893 2 229 2.16 2.57 3.872488e+17 3.8659872 0.05280850 3 250 5.19 1.99 3.957315e+17 0.8343263 0.31346937 4 261 8.21 2.41 2.339776e+18 5.1722473 0.04473758 > omega_h <- parallax_data$Plx > omega_g <- round(parallax_data$parallax, 2) > sigma_h <- parallax_data$e_Plx > sigma_g <- round(parallax_data$parallax_error, 2) > var_h <- sigma_h^2 > var_g <- sigma_g^2 > vector_muestral <- (omega_g-omega_h)^2/(var_g+var_h) > vector_muestral1 <- sort(vector_muestral) > discreto <- seq(1, 2045) > acumulado <- discreto/2045 > puntos <- pchisq(vector_muestral1, df=1) 89 > x <- seq(0, 1, length=10) > y <- x > plot(acumulado ~ puntos, type=’l’, main="P-Pplot para los paralajes" + , xlab="Funcion de distribucion teorica", + ylab="Funcion de distribucion muestral") > points(x,y, type=’l’, col=’blue’) > ks.test(vector_muestral1, "pchisq", 1) One-sample Kolmogorov-Smirnov test data: vector_muestral1 D = 0.042124, p-value = 0.00141 alternative hypothesis: two-sided > gdr2_proper_motion <- subset(gdr2, select=c("source_id", "pmra", "pmdec", + "pmra_error", "pmdec_error", "pmra_pmdec_corr")) > head(gdr2_proper_motion, 2) source_id pmra pmdec pmra_error pmdec_error pmra_pmdec_corr 1 3.245105e+15 24.50467 -269.44675 0.09872993 0.08102237 0.07712530 2 7.001557e+15 36.78335 -85.73746 0.51253975 0.42438689 -0.07672557 > gdr2_proper_motion <- subset(gdr2_proper_motion, + gdr2_proper_motion$pmra!="NA" & gdr2_proper_motion$pmdec!="NA" + & gdr2_proper_motion$pmra_error!="NA" & + gdr2_proper_motion$pmdec_error!="NA" & gdr2_proper_motion$pmra_pmdec_corr!="NA") > head(gdr2_proper_motion, 2) source_id pmra pmdec pmra_error pmdec_error pmra_pmdec_corr 1 3.245105e+15 24.50467 -269.44675 0.09872993 0.08102237 0.07712530 2 7.001557e+15 36.78335 -85.73746 0.51253975 0.42438689 -0.07672557 > nrow(gdr2_proper_motion) [1] 2045 > matching_proper_motion <- gdr2_hipp$source_id %in% + gdr2_proper_motion$source_id > comprobacion <- matching == matching_proper_motion 90 APÉNDICE B > which(comprobacion =="FALSE") integer(0) > previo_proper_motion <- previo_hipparcos > head(previo_proper_motion, 2) source_id original_ext_source_id 2 5.285634e+17 40 8 3.872488e+17 229 > nrow(previo_proper_motion) [1] 2045 Búsqueda SQL en Hipparcos SELECT ("5 magnitudes relativas a los movimientos propios") FROM I/239/hip_main WHERE HIP in (vector de identificadores de Hipparcos) > hipparcos_astrometic_p_m <- read.csv("2.4/hipparcos_astrom_p_m.csv" + , header=T) > head(gdr2_proper_motion, 2) source_id pmra pmdec pmra_error pmdec_error pmra_pmdec_corr 1 3.245105e+15 24.50467 -269.44675 0.09872993 0.08102237 0.07712530 2 7.001557e+15 36.78335 -85.73746 0.51253975 0.42438689 -0.07672557 > head(previo_proper_motion, 2) source_id original_ext_source_id 2 5.285634e+17 40 8 3.872488e+17 229 > head(hipparcos_astrometic_p_m, 2) HIP pmRA pmDE e_pmRA e_pmDE pmDE.pmRA 1 40 -2.99 -3.18 4.14 3.75 -0.10 2 229 18.80 -5.10 1.59 1.59 -0.08 91 > previo_proper_motion1 <- previo_proper_motion[with(previo_proper_motion, + order(previo_proper_motion$source_id)),] > gdr2_proper_motion1 <- gdr2_proper_motion[with(gdr2_proper_motion, + order(gdr2_proper_motion$source_id)),] > auxiliar_p_m <- merge(gdr2_proper_motion1, previo_proper_motion1) > head(auxiliar_p_m, 2) source_id pmra pmdec pmra_error pmdec_error pmra_pmdec_corr 1 3.245105e+15 24.50467 -269.44675 0.09872993 0.08102237 0.07712530 2 7.001557e+15 36.78335 -85.73746 0.51253975 0.42438689 -0.07672557 original_ext_source_id 1 15253 2 14075 > > colnames(auxiliar_p_m) <- c("source_id", "pmra", "pmdec", + "pmra_error","pmdec_error", "pmra_pmdec_corr", "HIP") > auxiliar_p_m1 <- auxiliar_p_m[with(auxiliar_p_m, order(auxiliar_p_m$HIP)),] > hipparcos_astrometic_p_m1<- hipparcos_astrometic_p_m[with(hipparcos_astrometic_p_m, + order(hipparcos_astrometic_p_m$HIP)),] > proper_motion_data <- merge(auxiliar_p_m1, hipparcos_astrometic_p_m1) > head(proper_motion_data, 2) HIP source_id pmra pmdec pmra_error pmdec_error 1 40 5.285634e+17 -1.721576 -2.388676 0.07204935 0.08050603 2 229 3.872488e+17 18.683311 -4.237698 0.07839015 0.05618416 pmra_pmdec_corr pmRA pmDE e_pmRA e_pmDE pmDE.pmRA 1 -0.2809894 -2.99 -3.18 4.14 3.75 -0.10 2 -0.3010555 18.80 -5.10 1.59 1.59 -0.08 > pmra_g_vector <- round(proper_motion_data$pmra, 2) > pmde_g_vector <- round(proper_motion_data$pmdec, 2) > pmra_sd_g_vector <- round(proper_motion_data$pmra_error, 2) > pmde_sd_g_vector <-round(proper_motion_data$pmdec_error, 2) > pmra_pmde_g_corr <- round(proper_motion_data$pmra_pmdec_corr, 2) > pmra_var_g_vector <- round(pmra_sd_g_vector^2, 2) > pmde_var_g_vector <- round(pmde_sd_g_vector^2, 2) 92 APÉNDICE B > pmra_pmde_g_cov <- (pmra_pmde_g_corr*pmra_sd_g_vector*pmde_sd_g_vector, 2) > pmra_h_vector <- proper_motion_data$pmRA > pmde_h_vector <- proper_motion_data$pmDE > pmra_sd_h_vector <- proper_motion_data$e_pmRA > pmde_sd_h_vector <- proper_motion_data$e_pmDE > pmra_pmde_h_corr <- proper_motion_data$pmDE.pmRA > pmra_var_h_vector <- pmra_sd_h_vector^2 > pmde_var_h_vector <- pmde_sd_h_vector^2 > pmra_pmde_h_cov <- pmra_pmde_h_corr*pmra_sd_h_vector*pmde_sd_h_vector > comp1_vector <- pmra_g_vector-pmra_h_vector > comp2_vector <- pmde_g_vector-pmde_h_vector > delta_p <- matrix(0, nrow=2045, ncol=2) > #Creamos el vector delta_p_m para cada estrella > for (i in seq(1, length(comp1_vector))) { + delta_p[i, 1]=comp1_vector[i] + delta_p[i, 2]=comp2_vector[i] + } > delta_p <- cbind(comp1_vector, comp2_vector) > head(delta_p, 2) comp1_vector comp2_vector [1,] 1.2684239 0.7913241 [2,] -0.1166886 0.8623024 > C_g <- cbind(pmra_var_g_vector, pmde_var_g_vector, pmra_pmde_g_cov ) > C_h <- cbind(pmra_var_h_vector, pmde_var_h_vector, pmra_pmde_h_cov ) > C <- C_g + C_h > head(C, 2) pmra_var_g_vector pmde_var_g_vector pmra_pmde_g_cov [1,] 17.144791 14.068981 -1.5541299 [2,] 2.534245 2.531257 -0.2035739 > colnames(C) <- c("pmra_v", "pmde_v", "pmra_pmde_cov") > muestra_p_m <- numeric(2045) > for (i in seq(1, 2045)) { + inv_det <- 1/(C[i,1]*C[i,2]-C[i, 3]^2) + muestra_p_m[i]=delta_p[i,1]^2*inv_det*C[i,2]+ delta_p[i,1]*delta_p[i,2]*2*inv_det*(-C[i,3])+ 93 delta_p[i,2]^2*inv_det*C[i,1] + + + + } > head(muestra_p_m) [1] 0.1528151 0.2946435 2.9473157 3.1405530 14.5896575 0.7103292 > length(muestra_p_m) [1] 2045 > muestra_p_m_ordenada <- sort(muestra_p_m) > discreto1 <- seq(1, 2045) > acumulado1 <- discreto/2045 > puntos1 <- pchisq(muestra_p_m_ordenada, df=2) > x <- seq(0, 1, length=10) > y <- x > plot(acumulado1 ~ puntos1, type=’l’, main="P-Pplot para los movimientos propios" + , xlab="Funcion de distribucion teorica", + ylab="Funcion de distribucion muestral") > points(x,y, type=’l’, col=’blue’) > ks.test(muestra_p_m_ordenada, "pchisq", 2) One-sample Kolmogorov-Smirnov test data: muestra_p_m_ordenada D = 0.10971, p-value < 2.2e-16 alternative hypothesis: two-sided Búsqueda SQL en el catálogo de Van Leeuwen SELECT HIP, Plx, e_Plx FROM "I/311/hip2" 94 APÉNDICE B WHERE HIP in (vector de identificadores de Hipparcos) > vanlewen <- read.csv("../Archivos TFG/2.4/vanlewen1.csv", header=T) > head(vanlewen, 2) HIP Plx e_Plx 1 40 -2.26 3.22 2 229 2.85 1.52 > colnames(vanlewen) <-c(’HIP’, ’paralaje’, ’error’) > head(vanlewen, 2) HIP paralaje error 1 40 -2.26 3.22 2 229 2.85 1.52 > vanlewen_analisis <- merge(vanlewen, parallax_data) > head(vanlewen_analisis, 6) HIP paralaje error Plx e_Plx source_id parallax parallax_error 1 40 -2.26 3.22 -3.40 4.25 5.285634e+17 1.0788363 0.05140893 2 229 2.85 1.52 2.16 2.57 3.872488e+17 3.8659872 0.05280850 3 250 4.01 1.05 5.19 1.99 3.957315e+17 0.8343263 0.31346937 4 261 4.98 1.61 8.21 2.41 2.339776e+18 5.1722473 0.04473758 5 274 1.31 0.38 0.93 0.57 4.316121e+17 0.2339128 0.48813736 6 316 5.67 0.65 5.48 1.25 2.766121e+18 5.4225870 0.06165799 > omega_h_v <- vanlewen_analisis$paralaje > omega_g_v <- round(vanlewen_analisis$parallax, 2) > sigma_h_v<- vanlewen_analisis$error > sigma_g_v <- round(vanlewen_analisis$parallax_error, 2) > var_h_v <- sigma_h_v^2 > var_g_v <- round(sigma_g_v^2, 2) > vector_muestral_v <- (omega_g_v-omega_h_v)^2/(var_g_v+var_h_v) > vector_muestral1_v <- sort(vector_muestral_v) > discreto_v <- seq(1, 2045) > acumulado_v <- discreto/2045 > puntos_v <- pchisq(vector_muestral1_v, df=1) > x <- seq(0, 1, length=10) > y <- x 95 > plot(acumulado_v ~ puntos_v, type=’l’, main="P-Pplot para los paralajes" + , xlab="Funcion de distribucion teorica", + ylab="Funcion de distribucion muestral") > points(x,y, type=’l’, col=’blue’) > ks.test(acumulado_v, ’pchisq’, 1) One-sample Kolmogorov-Smirnov test data: acumulado_v D = 0.31731, p-value < 2.2e-16 alternative hypothesis: two-sided > identificadores <- parallax_data$HIP > length(identificadores) [1] 2045 > write.table(identificadores, ’../Archivos TFG/2.4/identificadores.txt’) > # bucle > length(datos_dmsa1_dobles_fiables$HPmag) [1] 20696 > auxiliar<-numeric(0) > vectorbucle<-seq(2, 20696, by=2) > for(i in vectorbucle) + {auxiliar<-c(auxiliar, i)} > datos_comprobacion <- subset(datos_dmsa1_dobles_fiables[auxiliar,]) > head(datos_comprobacion, 6) CCDM Qual Ncomp Nparm Ncorr comp_id HIP HPmag RA 2 00003-4417 A 2 11 1 B 25 7.551 0.07924029 4 00004-4711 A 2 9 1 B 37 11.745 0.10532213 6 00005+6713 A 2 9 1 B 40 11.176 0.11781651 8 00005-7212 A 2 9 1 B 45 11.954 0.13192459 10 00006-5306 A 2 9 1 B 50 9.962 0.14241738 12 00006-6641 A 2 9 1 B 55 9.499 0.15516515 DE parallax pmRA pmDE theta rho RA2000 DE2000 102 APÉNDICE B + ylab=’numero de estrellas’, main=’’, xlim=c(-1, 1) ) > hist(difdis1 , breaks=80, col=’blue’, + xlab=’errores en la estimacion de la distancia(kpc)’ , + ylab=’numero de estrellas’, main=’’, xlim=c(-1,0.5)) 103 Código para 3.3 > datos<-read.csv(’../Archivos TFG/3.3/Makarov.csv’, header=TRUE) > head(datos) source_id 1 4.995999e+18 2 3.872488e+17 3 3.957315e+17 4 2.444844e+18 5 2.333016e+18 6 3.839349e+17 > source_id<-datos$source_id > write.table(source_id, file=’Makarov1.txt’) > head(source_id) [1] 4.995999e+18 3.872488e+17 3.957315e+17 2.444844e+18 2.333016e+18 [6] 3.839349e+17 Búsqueda ADQL SELECT parallax FROM gaiadr2.gaia_source WHERE parallax < 0 and source_id in ("vector de source_id") > datos_parallaxes<-read.csv(’../Archivos TFG/3.3/parallaxes.csv’) > head(datos_parallaxes) parallax 1 3.397189 2 6.048140 3 5.334009 4 7.551024 5 12.311582 6 10.523360 > parallaxes<-datos_parallaxes$parallax 104 APÉNDICE B > small_parallaxes<-parallaxes[parallaxes<1.8] > small_parallaxes_truncated<-small_parallaxes[small_parallaxes>0] > length(small_parallaxes_truncated)/length(small_parallaxes) [1] 0.8809524 > mean(small_parallaxes) [1] 0.6171852 > mean(small_parallaxes_truncated) [1] 1.011445 > par(mfrow=c(1,2)) > hist(small_parallaxes, breaks = 60,col=’green’, main=’Histograma de los + paralajes de estrellas lejanas’,xlab=’paralajes (mas)’, ylab=’numero de estrellas’) > hist(small_parallaxes_truncated, breaks=30, xlim=c(-2,2), col=’orange’, main=’Histograma + con los paralajes truncados’, xlab=’paralajes truncados (mas)’, ylab=’numero de estrellas’) Búsquedas ADQL según longitud galáctica SELECT * FROM gaiadr2.gaia_source WHERE l > 0 and l < 0.02 SELECT * FROM gaiadr2.gaia_source WHERE l > 0 and l < 0.02 $ parallax < 0 SELECT * 105 FROM gaiadr2.gaia_source WHERE l > 90 and l <90.02 SELECT * FROM gaiadr2.gaia_source WHERE l > 90 and l < 90.02 $ parallax < 0 SELECT * FROM gaiadr2.gaia_source WHERE l > 180 and l <180.02 SELECT * FROM gaiadr2.gaia_source WHERE l > 180 and l < 180.02 $ parallax < 0 SELECT * FROM gaiadr2.gaia_source WHERE l > 270 and l <270.02 SELECT * 106 APÉNDICE B FROM gaiadr2.gaia_source WHERE l > 2700 and l < 270.02 $ parallax < 0 107 Código para 3.4 > posterior1 <- function(r) { 1/(sqrt(2*pi)*0.001)*exp((-1/(2*0.001^2))*(0.01-1/r)^2)} > posterior2 <- function(r) {1/(sqrt(2*pi)*0.002)*exp((-1/(2*0.002^2))*(0.01-1/r)^2)} > posterior3 <- function(r) {1/(sqrt(2*pi)*0.005)*exp((-1/(2*0.005^2))*(0.01-1/r)^2)} > posterior4 <- function(r) {1/(sqrt(2*pi)*0.01)*exp((-1/(2*0.01^2))*(0.01-1/r)^2)} > r <- seq(45, 150, length=50) > y1 <- posterior1(r) > y2 <- posterior2(r) > y3 <- posterior3(r) > y4 <- posterior4(r) > plot(r, y1, type=’l’, col=’blue’, lwd=3, main=’POD para diferentes valores de f’, xlab=’r(pc)’, ylab=’POD’) > lines(r, y2, col=’green’, lwd=3) > lines(r, y3, col=’red’, lwd=3) > lines(r, y4, col=’orange’, lwd=3) > legend("topleft",col=c("blue","green", ’red’, ’orange’), legend =c("f=0.1","f=0.2", ’f=0.5’,’ f=1’), lwd=3, bty = "n") > rlim <- 1000 > posterior1 <- function(r) {(1/rlim)*1/(sqrt(2*pi)*0.001)*exp((-1/(2*0.001^2))*(0.01-1/r)^2)} > posterior2 <- function(r) {(1/rlim)*1/(sqrt(2*pi)*0.002)*exp((-1/(2*0.002^2))*(0.01-1/r)^2)} 108 APÉNDICE B > posterior3 <- function(r) {(1/rlim)*1/(sqrt(2*pi)*0.005)*exp((-1/(2*0.005^2))*(0.01-1/r)^2)} > posterior4 <- function(r) {(1/rlim)*1/(sqrt(2*pi)*0.01)*exp((-1/(2*0.01^2))*(0.01-1/r)^2)} > posterior5 <- function(r) {(1/rlim)*1/(sqrt(2*pi)*0.0025)*exp((-1/(2*0.0025^2))*(-0.01-1/r)^2)} > integrate(posterior1, lower=0, upper=1000) 10.31616 with absolute error < 4.2e-05 > integral1 <- 10.316 > integrate(posterior2, lower=0, upper=1000) 11.55355 with absolute error < 4.3e-06 > integral2 <- 11.554 > integrate(posterior3, lower=0, upper=1000) 27.61169 with absolute error < 0.00016 > integral3 <- 27.612 > integrate(posterior4, lower=0, upper=1000) 28.9847 with absolute error < 0.0013 > integral4 <- 28.985 > integrate(posterior5, lower=0, upper=1000) 0.002932285 with absolute error < 1.4e-07 > integral5 <- 0.00293 > posterior1_norm <- function(r) {(1/integral1)*(1/rlim)*1/(sqrt(2*pi)*0.001)* exp((-1/(2*0.001^2))*(0.01-1/r)^2)} > posterior2_norm <- function(r) 109 {(1/integral2)*(1/rlim)*1/(sqrt(2*pi)*0.002)* exp((-1/(2*0.002^2))*(0.01-1/r)^2)} > posterior3_norm <- function(r) {(1/integral3)*(1/rlim)*1/(sqrt(2*pi)*0.005)* exp((-1/(2*0.005^2))*(0.01-1/r)^2)} > posterior4_norm <- function(r) {(1/integral4)*(1/rlim)*1/(sqrt(2*pi)*0.01)* exp((-1/(2*0.01^2))*(0.01-1/r)^2)} > posterior5_norm <- function(r) {(1/integral5)*(1/rlim)*1/(sqrt(2*pi)*0.0025)* exp((-1/(2*0.0025^2))*(-0.01-1/r)^2)} > r <- seq(0, 1000, length=10000) > y1 <- posterior1_norm(r) > y2 <- posterior2_norm(r) > y3 <- posterior3_norm(r) > y4 <- posterior4_norm(r) > y5 <- posterior5_norm(r) > plot(r, y1, type=’l’, col=’blue’, main=’POD para diferentes valores de f’, xlab=’r(pc)’, ylab=’POD’) > lines(r, y2, col=’green’) > lines(r, y3, col=’red’) > lines(r, y4, col=’orange’) > lines(r, y5, col=’black’) > legend("topright",col=c("blue","green", ’red’, ’orange’, ’black’),legend =c("f=0.1","f=0.2", ’f=0.5’,’ f=1’,’|f|=0.25’), lwd=3, bty = "n") 110 APÉNDICE B > posterior1 <- function(r) {((r^2*exp(-r/L))/0.001)* exp((-1/(2*0.001^2))*(0.01-1/r)^2)} > posterior2 <- function(r) {((r^2*exp(-r/L))/0.002)* exp((-1/(2*0.002^2))*(0.01-1/r)^2)} > posterior3 <- function(r) {((r^2*exp(-r/L))/0.0029)* exp((-1/(2*0.0029^2))*(0.01-1/r)^2)} > posterior4 <- function(r) {((r^2*exp(-r/L))/0.0031)* exp((-1/(2*0.0031^2))*(0.01-1/r)^2)} > posterior5 <- function(r) {((r^2*exp(-r/L))/0.0033)* exp((-1/(2*0.0033^2))*(0.01-1/r)^2)} > posterior6 <- function(r) {1/2.5*((r^2*exp(-r/L))/0.005)* exp((-1/(2*0.005^2))*(0.01-1/r)^2)} > posterior7 <- function(r) {1/4*((r^2*exp(-r/L))/0.01)* exp((-1/(2*0.01^2))*(0.01-1/r)^2)} > posterior8 <- function(r) 111 {200*((r^2*exp(-r/L))/0.0025)* exp((-1/(2*0.0025^2))*(-0.01-1/r)^2)} > prior <- function(r) {13.5*r^2*exp(-r/L)} > L=10^3 > r <- seq(0,3000, length=200000) > y1 <- posterior1(r) > #y2 <- posterior2(r) > y3 <- posterior3(r) > #y4 <- posterior4(r) > y5 <- posterior5(r) > y6 <- posterior6(r) > #y7 <- posterior7(r) > y8 <- posterior8(r) > p <- prior(r) > par(ann=F) > plot(r, y1, type=’l’, col=’blue’ , lwd=1.5, main=’POD para diferentes valores de f’, xlab=’r(pc)’, ylab=’POD’) > #lines(r, y2, lwd=1.5, col=’green’) > lines(r, y3 , lwd=1.5, col=’red’) > #lines(r, y4, lwd=1.5) > lines(r, y5, lwd=1.5, col=’orange’) > lines(r, y6, lwd=1.5, col=’green’) > #lines(r, y7, lwd=1.5) > lines(r, y8, col=’black’, lwd=1.5) > lines(r, p, col=’yellow’, lwd=1.5) > legend("topright",col=c("blue","red", ’orange’ , ’green’, ’black’, ’yellow’),legend =c("f=0.1" ,"f=0.29", ’f=0.33’,’ f=0.5’,’|f|=0.25’, ’PD’), lwd=3, bty = "n") 118 APÉNDICE B [1] 73.65434 > sigma4 <- 0.00102 > omega4<- 0.00839 > L4 <- 720 > > integrand4 <- function(x) {(1/(sigma4*sqrt(2*pi)*2*L4^3))*(x^2)* + exp((-((omega4)-(1/x))^2)/(2*sigma4^2)-x/L4)} > > > > integrate(integrand4, lower=0, upper=120) 0.08844813 with absolute error < 2.1e-07 > integrate(integrand4, lower=20, upper=130) 0.1514656 with absolute error < 4.6e-09 > integrate(integrand4, lower=0, upper=250) 0.2654866 with absolute error < 7.2e-08 > integrate(integrand4, lower=0, upper=260) 0.2655031 with absolute error < 5.1e-07 > integrate(integrand4, lower=0, upper=270) 0.2655123 with absolute error < 2e-07 > integrate(integrand4, lower=0, upper=280) 0.2655176 with absolute error < 1e-06 > > norm4 <- 0.265518 > densidad4 <- function(x) {(1/norm4)* (1/(sigma4*sqrt(2*pi)*2*L4^3))*(x^2)* + exp((-((omega4)-(1/x))^2)/(2*sigma4^2)-x/L4)} > integrate(densidad4, lower=0, upper=300) 1.000017 with absolute error < 6.2e-06 > > integrand_mean4 <- function(x) {x*(1/norm4)* (1/(sigma4*sqrt(2*pi)*2*L4^3))*(x^2)* 119 + exp((-((omega4)-(1/x))^2)/(2*sigma4^2)-x/L4)} > integrate(integrand_mean4, lower=0, upper=300) 129.3829 with absolute error < 0.00045 > > > > prueba_error <- seq(0, 300, length=100000) > modas <- densidad4(prueba_error) > indice <- which.max(modas) > indice [1] 40867 > prueba_error[indice] [1] 122.5992 > > > > areas <- numeric(100000) > for (i in seq(1, 100000)) {areas[i]=integrate(densidad4, + lower=0, upper=prueba_error[i]) + + + } > > indices <- which(areas > 0.4999999) > indices[1] [1] 42305 > prueba_error[indices[1]] [1] 126.9133 > > > areas <- numeric(100000) > for (i in seq(1, 100000)) {areas[i]=integrate(densidad4, + lower=0, upper=prueba_error[i]) + + 120 APÉNDICE B + } > > indices <- which(areas > 0.049999) > indices[1] [1] 34704 > prueba_error[indices[1]] [1] 104.11 > > > areas <- numeric(100000) > for (i in seq(1, 100000)) {areas[i]=integrate(densidad4, + lower=0, upper=prueba_error[i]) + + + } > > indices <- which(areas > 0.949999) > prueba_error[indices[1]] [1] 144.7918 > sigma5 <- 0.00066 > omega5<- -0.00201 > L5 <- 496.31 > > integrand5 <- function(x) {(1/(sigma5*sqrt(2*pi)*2*L5^3))*(x^2)* + exp((-((omega5)-(1/x))^2)/(2*sigma5^2)-x/L5)} > > > > integrate(integrand5, lower=0, upper=100000) 121 0.2843173 with absolute error < 1e-04 > integrate(integrand5, lower=20, upper=110000) 0.2843174 with absolute error < 5.2e-07 > integrate(integrand5, lower=0, upper=120000) 0.2843174 with absolute error < 1.7e-06 > norm5 <- 0.284317 > densidad5 <- function(x) {(1/norm5) *(1/(sigma5*sqrt(2*pi)*2*L5^3))*(x^2)* + exp((-((omega5)-(1/x))^2)/(2*sigma5^2)-x/L5)} > integrate(densidad5, lower=0, upper=100000) 1.000001 with absolute error < 4.1e-07 > > integrand_mean5 <- function(x) {x*(1/norm5) *(1/(sigma5*sqrt(2*pi)*2*L5^3))*(x^2)* + exp((-((omega5)-(1/x))^2)/(2*sigma5^2)-x/L5)} > integrate(integrand_mean5, lower=0, upper=7200) 2676.42 with absolute error < 0.00035 > > > > prueba_error <- seq(0, 120000, length=100000) > modas <- densidad5(prueba_error) > indice <- which.max(modas) > indice [1] 1869 > prueba_error[indice] [1] 2241.622 > > > > areas <- numeric(100000) > for (i in seq(1, 100000)) {areas[i]=integrate(densidad5, + lower=0, upper=prueba_error[i]) + + + } 122 APÉNDICE B > > indices <- which(areas > 0.4999999) > indices[1] [1] 2114 > prueba_error[indices[1]] [1] 2535.625 > > > areas <- numeric(100000) > for (i in seq(1, 100000)) {areas[i]=integrate(densidad5, + lower=0, upper=prueba_error[i]) + + + } > > indices <- which(areas > 0.049999) > indices[1] [1] 1178 > prueba_error[indices[1]] [1] 1412.414 > > > areas <- numeric(100000) > for (i in seq(1, 100000)) {areas[i]=integrate(densidad5, + lower=0, upper=prueba_error[i]) + + + } > indices <- which(areas > 0.949999) > prueba_error[indices[1]] [1] 3280.771 123 Búsqueda SQL en Vizier para obtener los restantes valores de la Tabla 3.2 SELECT RAICRS, DEICRS, pmRA, pmDE, HIP FROM "I/239/hip_main" WHERE HIP in (2814, 274, 404, 760, 5844)