top of page

COSAS VEREDES

Continuando la parte del anterior post véase:



Se había conseguido con los programas de R proporcionados generar 2 vectores con distribuciones pseudo-aleatorias a lo largo de sus componentes de los números 1, 2 y 3.


En este sentido se había observado que mientras la distribución dada por R permitía dibujar el fractal conocido como triángulo de Sierpinski, sin embargo, se observó que las distribuciones de los vectores que finalmente generan las figuras, el basado en el algoritmo Mersenne – Twister que usa R por defecto y el basado en la sucesión de Fibonacci visualmente mostraban vectores con los números 1, 2 y 3 aparentemente bien distribuidos.


Por tanto, la pregunta es ¿Pueden los test estadísticos decidir con un número la presencia de aleatoriedad o no? Si se observa, por ejemplo, la distribución de los vectores generados, se obtienen los siguientes diagramas de barras:



Salvo un “pelín” de acumulación en el 2 cuando se trata de los números pseudo-aleatorios tipo Fibonacci básicamenonte la distribución es análoga.


No obstante, para contrastar realmente la aleatoriedad debe aplicarse algún test de rachas, que en R se podría hacer (de modo no paramétrico) mediante el uso de la función runs.test() de la librería tseries. Para aplicar este test de rachas se escogen los datos de los vectores que generan los anteriores gráficos y se aplica el modulo 2 a sendos vectores ya que esta función analiza rachas en variables binarias, por lo que si se hace tal transformación se alcanzan los siguientes resultados:



Donde básicamente se obtiene aleatoriedad en ambas secuencias de números al ser el p-valor muy elevado y no poderse rechazar la hipótesis de no aleatoriedad.


Otra comprobación adicional que se podría hacer es una basada en correlaciones mediante el establecimiento de un modelo de regresión lineal, en este sentido cabe generar una variable que se puede denominar t y que va desde 1 hasta 100000 (dado este caso) y que serviría para ordenar las observaciones. Así pues, si se ponen los datos obtenidos (números del 1 al 3) que genera el triángulo de Sierpinski y el otro en función de la variable t debería pasar que el la pendiente del modelo debería ser en todo caso nula y con mínima o nula significatividad, mientras que la constante debería tomar un valor próximo a 2 como de hecho ocurre:


Lo anterior nos muestra, que hay que tener mucho cuidado con los test estadísticos porque claramente no nos están dando una información suficientemente precisa y hay que investigar siempre más.


También cabe indicar que a la luz de los anteriores test son lo mismo las 2 sucesiones (incluso en esta última prueba la distribución tipo Fibonacci podría decirse que es incluso más aleatoria), pero gráficamente hemos visto que no tienen nada que ver por el modo de escoger los 3 vértices para generar los puntos del triángulo de sierpinski, por tanto, y para concluir, ojito que no todo lo que hay que saber de estadística (o sobre data scientist) está en un único libro (ni mucho menos en Chat GPT) y cosas aparentemente inconexas pueden estar más conectadas de lo que creemos.


PD Intenté usar Chat GPT para esto, pero si en algo me ayudó, se confunde un montón con el tema de los test y los resultados dando respuestas falsas que continuamente corrige y al final te hace ir a las fuentes y a libros fiables.

  • fjroar
  • 25 ene
  • 4 Min. de lectura

Investigando y repasando el mundo de los fractales, me viene a la mente algunas ideas donde se mezcla por un lado la frialdad del mundo de la aleatoriedad y por otro la belleza de la geometría de Mandelbrot.


En este sentido, acabando el libro de Fractals for de Classroom, me encuentro en su capítulo 6 observaciones sobre el cuidado que hay que poner a la hora de elegir un generador de números aleatorios, elemento éste muy utilizado en aplicaciones reales cuando se quieren realizar, por ejemplo, simulaciones tipo Montecarlo en riesgo de mercado.


En este sentido, no sé hasta donde se tiene idea de cómo se generan los números aleatorios en un ordenador, pero desde luego son todo menos aleatorios y se suelen denominar más propiamente como “pseudo-aleatorios”.


Así pues, una prueba geométrica para ver si un sistema cumple con esa condición de aleatoriedad, es ver si es capaz de generar o no un triángulo de Sierpinski tal y como se muestra a continuación:


La idea geométrica para realizar se basa en dar los siguientes pasos:


-        Se fijan 3 vértices de un triángulo que en este caso son el (0; 0), el (4; 0) y el (2; 3’46)

-        Posteriormente se elije un punto inicial como pueda ser el (2; 1)

-        Y se generan el resto de los puntos siguiendo la regla:

 


Donde V es uno de los anteriores 3 vértice que se escoge aleatoriamente en cada

iteración con n suficientemente grande, aquí se supone que n va desde 1 a 100.000

 

Pues bien, si os entretenéis en ejecutar el siguiente programa se obtiene la anterior figura:

 

set.seed(123)

n_points <- 100000

vertices <- matrix(c(0, 0,

                     4, 0, 

                     2, 3.46),

                   ncol = 2, byrow = TRUE)

current_point <- c(2, 1)

sierpinski_points <- matrix(NA, nrow = n_points, ncol = 2)

guardaVerticesAleatorio <- numeric(n_points)

 

for (i in 1:n_points) { 

  guardaVerticesAleatorio[i] <- sample(1:3, 1)

  chosen_vertex <- vertices[guardaVerticesAleatorio[i], ] 

  current_point <- (current_point + chosen_vertex) / 2

  sierpinski_points[i, ] <- current_point

}

plot(sierpinski_points, pch = ".", col = "blue", xlab = "", ylab = "", asp = 1,

     main = "Triángulo de Sierpinski")

 

La ejecución del anterior código sólo requiere Rbase y ante esto uno puede pensar que lo anterior es sencillo y puede que hasta haya escuchado hablar de generadores de números pseudoaleatorios basados en sucesiones con operaciones modulares del tipo:

                                                                                                                                                       

 

O incluso algunas de tipo más sofisticado como las de tipo Fibonacci:



Pensando en un m suficientemente grande como 2^18 para evitar ciclos sistemáticos.

 

Entonces ¿Qué sucedería si tomamos uno de tipo Fibonacci tal y como se generaría a continuación?:

 

# Generador de Fibonnacci

m <- 2^18         

r <- numeric(n_points)   

r[1] <- 5

r[2] <- 7

for(k in c(3:100000)){

  r[k] <- (r[k-2] + r[k-1]) %% m

}

# Transformar los valores al rango deseado {1, 2, 3}

r_mod <- (r %% 3) + 1

 

Pues bien, lo que sucedería es que a pesar de que esta secuencia tendría este aspecto:

Y la generada por R tendría este otro aspecto:


No haciendo pensar aparentemente en ninguna diferencia en cuanto a la “independencia” de estas generaciones.


Sin embargo, si se ejecuta el programa completo basado en el generador de Fibonacci (ver anexo), el resultado es el siguiente (nótese la diferencia en la altura que en este gráfico llega a 3):



Es decir, esta serie de números no es capaz de rellenar el triángulo de Sierpinski tal y como se indica el mencionado libro.


Por tanto resulta curioso cómo cuando se consigue estar cerca de la aleatoriedad con un algoritmo suficientemente fiable se es capaz de tener una figura de cierta estética, mientras que cuando nos apartamos del “caos” el resultado aunque ofrece cierta regularidad, no es capaz de aportar completitud, por lo que resulta importante conocer si los algoritmos que sustentan la generación de números pseudo-aleatorios son suficientemente fiables, en el caso de R el módulo del algoritmo (por defecto se usa el Mersenne – Twister) es igual a 2^19937 – 1 dando una idea de la complejidad del problema frente al anterior basado en Fibonacci.

 

Anexo: Segunda parte del código R al completo:


vertices <- matrix(c(0, 0,  # Vértice 1

                     4, 0,  # Vértice 2

                     2, 3.46), # Vértice 3 (altura de un triángulo equilátero)

                   ncol = 2, byrow = TRUE)

current_point <- c(2, 1)

sierpinski_points <- matrix(NA, nrow = n_points, ncol = 2)

r = 0

m <- 2^18         

r <- numeric(n_points)

r[1] <- 5

r[2] <- 7

for(k in c(3:100000)){

  r[k] <- (r[k-2] + r[k-1]) %% m 

}

r_mod <- (r %% 3) + 1

 

for (i in 1:n_points) {

  chosen_vertex <- vertices[r_mod[i], ]

  current_point <- (current_point + chosen_vertex) / 2

 

  sierpinski_points[i, ] <- current_point

}

plot(sierpinski_points, pch = ".", col = "blue", xlab = "", ylab = "", asp = 1,

     main = "Triángulo de Sierpinski")

Aprovechando que estoy investigando algoritmos que tienen que ver con la Luna, me encuentro con uno en particular, muy sencillo y que es la base de la mayoría de los algoritmos astronómicos como es el cálculo del día juliano.


Aunque no me voy a poner describir este algoritmo en este momento, lo que me llama la atención es el gran uso que hace de una función que se enseña en lo que ahora es segundo o tercero de la ESO (o al menos eso creo y espero) y antes era séptimo u octavo y que me parece incluso hasta más interesante que la fórmula de I = Crt/100 que para los preocupados en que se enseñe como pedir créditos a los chavales (véase algún comentario de linkedin anterior que hice), también se enseña por ese tramo de edad.

Pues bien, la gráfica de la función resulta ser la siguiente:


Imagen generada en R con plotly usando la función floor()
Imagen generada en R con plotly usando la función floor()

Así pues, como la mayoría de la gente sabe, por ejemplo, la parte entera de -1.65 es -2 y si se quiere mayor detalle se puede consultar en https://eliatron.blogspot.com/2012/05/cuanto-vale-la-parte-entera-de-165.html

Sin embargo, el problema que se tiene en Astronomía es que el grado de precisión importa y por ejemplo, si hago la pregunta de ¿Cuál es la parte entera de -2.001? Estoy seguro que casi todo el mundo acierta y me dice -3, de hecho, si en R hago eso, observo lo mismo con tan sólo teclear lo siguiente:

                                                           floor(-2.001)


O al que le guste Python (que también los hay y tendrán que llamar a una librería):


                                                           import math

math.floor(-2.01)

 

Sin embargo ¿Cuál sería el resultado de lo anterior con:


                                                           floor(-2.0000000000000001)

O con:

                                           import math

math.floor(-2.0000000000000001)

 

La respuesta sorprendente, al menos en mi máquina de de 64bit, hp de 32gb de ram es que el resultado es -2 y sólo hay 16 decimales tras la coma, lo cual en algoritmos de alta precisión os podéis imaginar qué sucede.


Lo anterior, me obligó en la función que tengo par el cálculo de días julianos en mi librería RMoon y que podéis encontrar en https://github.com/FJROAR/RMoon/blob/main/R/JulianDayMeeus.R  a definir una función especial (por recomendación del libro de Jean Meeus y del responsable del Grupo de Observación Lunar de la Agrupación Astronómica de Madrid https://www.aam.org.es/, el gran Alberto Martos a programar algo específico para evitar errores que diese problemas en el cálculo del día juliano y por tanto para ser preciso en el algoritmo encontráis la función goodInteger():

   

goodInteger <- function(number){

for (i in c(1: length(number))){

  number[i] <- ifelse(number[i] < 0, as.integer(number[i]) - 1, as.integer(number[i]))

}


return(number)


}

 

Que soluciona este problema y ofrece el valor correcto pero ¿Por qué es tan importante calcular bien la función parte entera de un número? Pues bien, en el algoritmo del cálculo del día juliano, cuyo concepto parte de Justus Scaliger (1540 – 1609) en 1583, existen algunas piezas claves, (y que incorporo en RMoon) siendo una de las importantes la siguiente:

 

JD <- goodInteger(365.25 (Y + 4716)) + goodInteger(30.6001 (M + 1)) + D + B - 1524.5

 

No voy a hablar del algoritmo en sí y de todo lo que significa la anterior línea, pero sí me voy a detener en el primer sumando y en la última resta, es decir, en el último número que se descompone como 1524.5 = 1461 + 62 + 1.5. Esta línea se puede re-escribir del siguiente modo para separar los efectos en 2 partes distintas:

 

JD1 <- goodInteger(365.25 * (Y + 4716)) – (1461 + 1.5)

JD2 <- goodInteger(30.6001 *(M + 1)) + D + B – 62

 

Fijémonos en JD1, ya que JD2 sirve para otros temas que tienen que ver con la reforma de Gregorio XIII de 1582.


Partiendo de la base de que Y es un año cualquiera. ¿Qué sucede cuando por ejemplo es el año 2022? Pues en ese caso el valor que toma JD1 sería:


[365.25 * (2022+4716)] – 1462.5 = [2461054.5] – 1462.5 = 2461054 – 1462.5 = 2459591.5


Es decir, al no ser 2022 bisiesto el decimal 0.5 no tiene efecto y se anula y no se cuenta y se acaban contando 365 días que se resta a una cantidad fija, pero es que si se considera por ejemplo el año 2023, el resultado sería:


[365.25 * (2023+4716)] – 1462.5 = [2461419.75] – 1462.5 = 2460689 – 1462.5 = 2459956.5


Y ahora, que el redondeo natural sería una cifra más por ser 0.75 lo que hay dentro de la función parte entera, también desaparece y se sigue contando esos 365 día, pero fijémonos que cuando se alcanza el año 2024, sucede lo siguiente:


[365.25 * (2024+4716)] – 1462.5 = [2461785] – 1462.5 = 2461785 – 1462.5 = 2460322.5


¡Ahora no se eliminan decimales! Porque cada 4 años, la parte 0.25 de 365.25 aumenta un día y ese día se mantiene inalterado mientras que, en el resto de los casos, ese 0.25 no tiene efecto y así con una simple línea donde la función parte entera es esencial se van contando los años (¡Ojo! Que como he mencionado hay más cosas en el algoritmo pero quería mostrar la ingeniosidad de la esta función).


Además, se puede comprobar que como no podía ser de otro modo que la diferencia de días entre los años se mantendría del siguiente modo:


                                             2460322.5 – 2459956.5 = 366

                                             2459956.5 – 2459336.5 = 365


Por tanto la función parte entera es muy útil para la configuración de modo eficiente e ingeniosa de algoritmos como el aquí presentado, pero que nos advierte de que hay que tener cuidado con las funciones que vienen de los lenguajes a nivel base o en librerías estándar y que todo el mundo asume que hace bien los cálculos pero que no alcanza el grado de exactitud que se requiere sobre todo en los puntos límite y vuelvo a recordar que en la época actual con el hype de la Inteligencia Artificial donde las funciones ReLu que componen muchos de nuestros algoritmos y que usan funciones similares a la parte entera.


Miedo me da, cuando hay miles de estas funciones entremezcladas donde no se ha hecho un estudio adecuado de si en casos límite, la cosa va o no, después vienen las “alucinaciones”, los errores de predicción inexplicables, … pero si fallamos en la base, después es difícil encontrar por qué todo fue mal.


Finalmente quisiera acabar con que hay que investigar y aprender y esto no va sólo apostar y arriesgar patrimonios y si no lo haces eres un looser, admiro a quien lo hace, pero también admiro al que estudia e indaga y no se queda con la primera respuesta, le da vueltas a las cosas y también se equivoca, este último es muy posible que no llegue a ser millonario, pero tiene altas probabilidades de disfrutar con lo que hace, e investigar el día juliano te lleva a analizar cómo funciona el tiempo, los calendarios, cómo se establece la hora civil, … y eso te llena el alma, aunque no la cuenta bancaria.

© 2021 by Francisco J. Rodríguez Aragón. Proudly created with Wix.com

bottom of page