Siguiendo con la idea de obtener polinomios que se ajusten a los datos para calcular la pendiente de dicho polinomio y hallar el valor cero de la pendiente (lo cual indica que la curva ha pasado por un máximo o un mínimo), rebusqué por Internet y hallé el programa Tablecurve, el cual es todo lo que necesitamos para obtener la curva que se ajusta mejor a nuestros datos fotométricos. En el caso de estrellas eclipsantes con simetría perfecta o cuasiperfecta, basta usar el método de Kwee van Voerden para obtener un valor muy exacto, pero para curvas de estrellas pulsantes donde la simetría no es perfecta entre la subida y bajada no da muy buenos resultados.
Así que mi idea de derivada no era descabellada y aunque intenté escribir un programa en qbasic comprobé que ya existía todo lo que necesitaba.
El modelado escogido de la curva es por medio del algoritmo de Savitzky Golay, que consiste en ajustar una recta de regresión móvil (polinomio de grado 1) o un polinomio de grado 2, 4, etc... por mínimos cuadrados a los datos.
Los parámetros a combinar son los siguientes:
El número de puntos que usaremos para hallar la regresión. A más puntos implicados más suave es la curva que sale, pues los puntos dispares pierden fuerza frente a una mayoría homogénea. Ésto se ve en los siguientes gráficos, donde se han escogido 5, 10 y 15 puntos para hacer la regresión.
Otra opción es cambiar el orden del polinomio de regresión (usar polinomios de segundo, cuarto grado, etc...). Así podemos ver lo que pasa al aumentar el grado de polinomio de ajuste. Aquí vemos la curva con 10 puntos y grado dos de polinomio.
Y con 10 puntos y grado cuatro para el polinomio.
Al aumentar el grado del polinomio aparecen más curvas que se ajustan mejor a los datos.
Otra opción es repetir la regresión varias veces, de modo que hacemos una segunda regresión sobre los datos obtenidos de la primera regresión, suavizándose cada vez más la curva resultante.
Aquí vemos un ajuste con una , dos y tres pasadas, manteniendo el número de puntos (15 puntos) en la regresión móvil y el mismo grado de polinomio de ajuste (cuatro).
Partimos de que nuestro datos tienen una dispersión notable, por lo que más que ajustarse a todos los saltos que se deriven de la dispersión de los datos buscamos una curva "suave" que se ajuste a la tendencia "general" de la curva observada. Una vez que decidimos que un "grupo" de curvas se ajustan aceptablemente a los datos, podemos buscar la derivada de la curva para obtener el instante del extremo.
La curva de la pendiente (derivada) de la curva polinómica ajustada a los datos tiene un aspecto como la siguiente.
Cuando la curva pasa por el valor cero tenemos el instante exacto del extremo de la curva. Dicho valor se puede obtener numéricamente.
El problema es que dicho valor va cambiando según los ajustes que modifiquemos, pero podemos hacer finalmente un ejercicio de estadística para calcular el valor medio de todas las curvas (o aquellas que no den muchos saltos por ajustarse demasiado a los errores de las medidas).
En mi caso obtuve los siguientes valores:
Spline Savitzky Golay 19 points order 4 passes 3 7.4865827
Spline Savitzky Golay 19 points order 2 passes 3 7.4865359
Spline Savitzky Golay 19 points order 2 passes 3 7.4865359
Spline Savitzky Golay 15 points order 2 passes 2 7.4861165
Spline Savitzky Golay 15 points order 2 passes 3 7.4859521
Spline Savitzky Golay 12 points order 2 passes 3 7.4866057
Spline Savitzky Golay 12 points order 2 passes 2 7.4868928
Spline Savitzky Golay 12 points order 2 passes 1 7.4867535
Spline Savitzky Golay 10 points order 2 passes 1 7.4890655
Spline Savitzky Golay 10 points order 2 passes 2 7.4883725
Spline Savitzky Golay 12 points order 2 passes 1 7.4867535
Spline Savitzky Golay 10 points order 2 passes 1 7.4890655
Spline Savitzky Golay 10 points order 2 passes 2 7.4883725
Spline Savitzky Golay 10 points order 2 passes 3 7.4880479
Spline Savitzky Golay 08 points order 2 passes 2 7.4901273
Spline Savitzky Golay 08 points order 2 passes 3 7.4897913
El máximo oscila entre 7.486 y 7.490 , de estos valores podría calcular el valor medio y la dispersión ya que no existe una curva que podamos decir que es la verdadera. El coeficiente de ajuste R que da el programa para cada ajuste va variando (0.94, 0.95, ...) y se aproxima a uno para polinomios de grado 4 , pero esos polinomios van siguiendo las dispersiones que podamos obtener al tomar las medidas por lo que no es fiable al 100% ese parámetro para elegir la "mejor" curva.
En general parece que usar polinomios de grado dos es más efectivo que el grado cuatro y que realmente debemos usar del orden de 8 a 15 puntos para obtener un valor promedio (si disparamos una foto cada minuto, el valor medio sale de 8 a 15 fotos, o lo que es lo mismo obtenermos un valor medio a intervalos de 8 a 15 minutos). Y sobre esos valores medios se va construyendo la curva suavizada.
En fin, ésto de los extremos está entre lo puramente científico y el arte del término medio. Al final nuestro ojo debe determinar que curvas son más "reales".
El programa tiene muchísimas opciones y yo poco tiempo de estudiarlo pero creo que es bastante efectivo para lo que nos interesa (encontrar extremos de conjunto de valores).
Como otros observadores en otras áreas astronómicas, quizás sea más relevante seguir los pasos para obtener medidas con dispersiones muy bajas (usar estrellas desenfocadas, una buena altura sobre el horizonte de la variable , etc...). Esto permitiría obtener curvas que no se distinguan apenas entre sí, dando un valor de extremo muy parecido todas ellas.
No hay comentarios:
Publicar un comentario