Dibujando la Aleatoriedad

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")