pi[1] 3.141593format(pi, digits = 8) # Con más dígitos[1] "3.1415927"Una visualización para aproximar el valor de pi
Y para aprender a programar funciones
El problema
En A First Course in Machine Learning (Rogers and Girolami 2016, 164) hay un ejercicio bonito:
El ejercicio consiste en ingeniar un método que aproxime el valor de \(\pi\) a partir de una muestra de datos (puntos, tal como veremos pronto); estos datos deben seguir una distribución uniforme.
¿Ninguna experiencia previa programando en R?
Aquí publiqué una pre-introducción a R para personas que nunca antes han programado.
El valor de \(\pi\) es \(3.1415927\). pi existe en R:
Pero no podremos usar pi directamente en nuestro método pues, de hecho, el objetivo es aproximarlo.
Los paquetes que necesitamos para llevar a cabo este ejercicio y reproducir todo el código son los siguientes:
library(dplyr)
library(ggplot2)
library(glue)
library(patchwork)Recordemos que este ejercicio requiere producir datos que estén uniformemente distribuidos.
Es relevante saber que, por ejemplo, las funciones rnorm() y runif() permiten generar datos que sigan una distribución normal y una distribución uniforme, respectivamente:
## Producir datos que sigan un distribución normal
normal <- rnorm(n = 1000000)
## Visualizar (histograma)
ggplot() +
geom_histogram(
aes(x = normal),
fill = "#fff2f9",
color = "#c90076"
)
## Producir datos que sigan una distribución uniforme
uniforme <- runif(n = 1000000, min = -10, max = 10)
## Visualizar (histograma)
ggplot() +
geom_histogram(
aes(x = uniforme),
fill = "#fff2f9",
color = "#c90076"
)
Pensemos en los datos como coordenadas. Para ubicar puntos en un plano necesitamos sus coordenadas en los ejes X y Y:
## Establecer el número de puntos
n <- 3
## Producir coordenadas en eje X
(x <- runif(n = n, min = -10, max = 10))[1] -5.184785 -7.565158 1.161610## Producir coordenadas en eje Y
(y <- runif(n = n, min = -10, max = 10))[1] 3.578837 3.916641 9.540349## Visualizar puntos
ggplot() +
geom_point(
aes(x = x, y = y),
shape = 23,
size = 4,
fill = "#fff2f9",
color = "#c90076"
) +
xlim(-10, 10) +
ylim(-10, 10)
Los tres puntos que miramos arriba los situamos en el plano en función de distribuciones uniformes.
Esto es todo lo necesario para iniciar este ejercicio.
El método
El método consiste en crear una visualización que ilumine una manera de aproximar el valor de \(\pi\) a partir de datos distribuidos uniformemente a lo largo y ancho de un plano con coordenadas polares.
Sabemos, para empezar, que \(\pi\) interviene en el cálculo del área del círculo: \(\pi r^2\). O sea, intuimos que de alguna forma tiene que haber un círculo involucrado en la solución de este ejercicio.
Vamos a dibujar un círculo. Imaginemos tener un compás en nuestras manos mientras miramos la siguiente visualización:
plot <- ggplot() +
## Segmento vertical
geom_segment(
aes(
x = 0, y = -10,
xend = 0, yend = 10
),
color = "black",
linetype = "dashed"
) +
## Segmento horizontal
geom_segment(
aes(
x = -10, y = 0,
xend = 10, yend = 0
),
color = "black",
linetype = "dashed"
) +
## Punto central
geom_point(
aes(x = 0, y = 0),
size = 3,
color = "#00c953",
shape = 19
) +
## Punto izquierdo
geom_point(
aes(x = -10, y = 0),
size = 2,
color = "gold"
) +
## Punto derecho
geom_point(
aes(x = 10, y = 0),
size = 2,
color = "gold"
) +
## Punto bajo
geom_point(
aes(x = 0, y = -10),
size = 2,
color = "gold"
) +
## Punto alto
geom_point(
aes(x = 0, y = 10),
size = 2,
color = "gold"
) +
xlab("") + # Sin título en eje X
ylab("") # Sin título en eje Y
plot
Si esto fuese una manualidad, clavaríamos el compás en el punto verde, ubicado justo en el centro del plano, y dibujaríamos el círculo conectando los puntos amarillos con el lápiz del compás.
Para dibujar el círculo en R vamos a implementar una función publicada en Stack Overflow. Una vez trazado, el círculo se vería exactamente así:
## Crear la función plot_circle()
plot_circle <- function(
center = c(0,0), # El centro
diameter = 20, # El diámetro
npoints = 100 # El número de puntos
){
## Radio
r <- diameter / 2
## Coordenadas
tt <- seq(0,2 * pi, length.out = npoints)
xx <- center[1] + r * cos(tt)
yy <- center[2] + r * sin(tt)
## Data frame con las coordenadas
df <- data.frame(xx = xx, yy = yy)
## Output de la función: df
return(df)
}
## Crear data frame con las coordenadas del círculo
df_circle <- plot_circle()
## Imprimir el círculo
plot1 <- plot +
geom_path(
aes(
df_circle$xx,
df_circle$yy
),
color = "black"
)
plot1
Del ejemplo anterior debemos sacar tres importantes conclusiones:
Los segmentos punteados van de \(-10\) a \(10\), en ambos ejes;
Por tanto, el diámetro del círculo es \(20\);
- Por tanto, el radio del círculo es \(10\).
Ahora vamos a generar puntos a partir de distribuciones uniformes y con esos puntos vamos a poblar la visualización.
Programar funciones en R
Estamos por programar una función, plot_points(), que cumplirá varios propósitos:
Imprimir puntos.
Imprimir segmentos.
Imprimir el círculo.
Calcular la hipotenusa de cada punto.
Distinguir, con base en la hipotenusa, cuáles puntos están ubicados dentro del círculo.
Aproximar el valor de \(\pi\).
A efectos de este ejercicio vamos a programar la función plot_points(). Esta función reúne todo lo necesario para crear distintos tipos de visualizaciones:
## Crear la función plot_points()
plot_points <- function(
n = 5, # El número de puntos
hyp = FALSE,
seg = FALSE,
cir = FALSE
){
## Límites del diámetro
min <- -10
max <- 10
## La semilla, para replicar los resultados
set.seed(10)
## Coordenadas en eje X
x <- runif(n = n, min = min, max = max)
## Coordenadas en eje Y
y <- runif(n = n, min = min, max = max)
## Data frame con las coordenadas de los puntos
## También calcula la hipotenusa
df <- data.frame(
x = x, # Coordenadas en eje X
y = y, # Coordenadas en eje Y
hyp = sqrt((x ^ 2) + (y ^ 2)) # Hipotenusa
)
## Conteo de puntos ubicados dentro del círculo
## Usando la hipotenusa
inside <- df |>
dplyr::filter(hyp <= 10) |>
nrow()
## Proporción de puntos = (dentro / total) * 4
## Esto es, la aproximación de pi
prop <- format(
((inside / nrow(df)) * 4),
digits = 5
)
## Argumento: hyp
if(hyp == TRUE){
## Distinguir puntos dentro y fuera del círculo
## Usando la hipotenusa
plot <- ggplot(data = df, aes(x = x, y = y)) +
geom_point(
color = ifelse(df$hyp > 10, "#c90076", "#0076c9"),
shape = ifelse(df$hyp > 10, 4, 1)
) +
xlim(min, max) +
ylim(min, max) +
ggtitle( # Pi aproximado se imprime en el título
glue("N = {nrow(df)}, pi aproximado = {prop}")
)
} else{
## Imprimir baseline plot
plot <- ggplot(data = df, aes(x = x, y = y)) +
geom_point(shape = 15, color = "#c90076") +
xlim(min, max) +
ylim(min, max)
}
## Argumento: seg
if(seg == TRUE){
## Imprimir segmentos
plot <- plot +
## Segmento vertical
geom_segment(
aes(
x = 0, y = -10,
xend = 0, yend = 10
),
color = "black",
linetype = "dashed"
) +
## Segmento horizontal
geom_segment(
aes(
x = -10, y = 0,
xend = 10, yend = 0
),
color = "black",
linetype = "dashed"
) +
## Punto central
geom_point(
aes(x = 0, y = 0),
size = 3,
color = "#00c953",
shape = 19
) +
## Punto izquierdo
geom_point(
aes(x = -10, y = 0),
size = 2,
color = "gold"
) +
## Punto derecho
geom_point(
aes(x = 10, y = 0),
size = 2,
color = "gold"
) +
## Punto bajo
geom_point(
aes(x = 0, y = -10),
size = 2,
color = "gold"
) +
## Punto alto
geom_point(
aes(x = 0, y = 10),
size = 2,
color = "gold"
)
} else{
## Imprimir baseline plot
plot
}
## Argumento: cir
if(cir == TRUE){
## Imprimir círculo
plot <- plot +
geom_path(
data = plot_circle(),
aes(
xx,
yy
),
color = "black"
)
} else{
## Imprimir baseline plot
plot
}
## Output de la función: plot
return(plot)
}Creada la función plot_points(), empecemos por visualizar nada más que puntos:
## Visualizar 5 puntos
plot_points(n = 5)
## Visualizar 50 puntos
plot_points(n = 50)
## Visualizar 500 puntos
plot_points(n = 500)
## Visualizar 5000 puntos
plot_points(n = 5000)
## Visualizar 50000 puntos
plot_points(n = 50000)
Es bastante claro que los puntos van formando un cuadrado. Cuantos más puntos imprimimos, más evidente es el cuadrado.
Además, si volvemos a imprimir el círculo, nos resultará claro que algunos puntos han quedado fuera de él (los puntos ubicados en las cuatro esquinas, por decirlo así):
plot_points(n = 50000, seg = TRUE, cir = TRUE)
Lo que tenemos ahora es un círculo inscrito en un cuadrado: el área del cuadrado es el todo; el área del círculo inscrito en el cuadrado es una parte del todo.
El método para aproximar el valor de \(\pi\) que implementaremos a continuación explota precisamente la relación entre ese cuadrado y el círculo inscrito en él.
Estamos hablando de una proporción: el área del cuadrado es el todo; el área del círculo inscrito en el todo es -you guessed it- una parte del todo.
Por tanto, si logramos identificar y contar cuántos puntos están dentro del círculo, podremos entonces calcular una proporción: los puntos dentro del círculo en relación con el número total de puntos.
Estamos por aprender que la hipotenusa es la clave para determinar si un punto está situado dentro o fuera del círculo.
La hipotenusa
Según Wikipedia:
“La hipotenusa es el lado opuesto al ángulo recto en un triángulo rectángulo, y resulta ser su lado de mayor longitud.
Según el teorema de Pitágoras, el cuadrado de la longitud de la hipotenusa es igual a la suma de los cuadrados de las respectivas longitudes de los otros dos lados del triángulo rectángulo, denominados catetos.”
Pero vamos a iniciar con un ejercicio algo más simple, imprimiendo nada más dos puntos:
plot2 <- plot1 +
geom_point( # El punto rojo
aes(x = -5, y = 5),
color = "red",
shape = 15,
size = 3
) +
geom_point( # El punto púrpura
aes(x = 7.5, y = -7.5),
color = "purple",
shape = 15,
size = 3
)
plot2
En la visualización anterior, notemos que el punto rojo está dentro del círculo, mientras el punto púrpura está fuera del círculo.
Sabemos bien que las coordenadas en los ejes X y Y nos indican la posición de los puntos en el plano (en nuestro ejercicio original, estas coordenadas las generamos a partir de distribuciones uniformes).
Igualmente podríamos pensar en dichas coordenadas como si fuesen traslados a lo largo y ancho del plano, partiendo del centro:
plot3 <- plot2 +
## Cateto en Y, punto rojo
geom_segment(
aes(
x = -5, y = 5,
xend = -5, yend = 0
),
color = "red"
) +
## Cateto en X, punto rojo
geom_segment(
aes(
x = -5, y = 0,
xend = 0, yend = 0
),
color = "red"
) +
## Cateto en Y, punto púrpura
geom_segment(
aes(
x = 0, y = -7.5,
xend = 0, yend = 0
),
color = "purple"
) +
## Cateto en X, punto púrpura
geom_segment(
aes(
x = 0, y = -7.5,
xend = 7.5, yend = -7.5
),
color = "purple"
)
plot3
Arriba, al dibujar los traslados necesarios para alcanzar cada punto partiendo del centro, básicamente dibujamos los dos lados que formarían el ángulo recto de un triángulo rectángulo (los catetos).
El único lado que nos falta para, en efecto, conseguir triángulos rectángulos es la hiponetenusa.
Agreguemos entonces ambas hipotenusas. Las vamos a hacer ligeramente más gruesas, para distinguirlas de los catetos:
plot3 +
## Hipotenusa, punto rojo
geom_segment(
aes(
x = -5, y = 5,
xend = 0, yend = 0
),
color = "red",
linewidth = 1.5
) +
## Hipotenusa, punto púrpura
geom_segment(
aes(
x = 7.5, y = -7.5,
xend = 0, yend = 0
),
color = "purple",
linewidth = 1.5
) 
La visualización anterior arroja un hallazgo crucial:
- Si la hipotenusa es mayor que el radio del círculo (\(10\), en nuestro caso), ello quiere decir que el punto en cuestión está fuera del círculo.
En suma, si calculamos la hipotenusa de todos los puntos y las filtramos en función de \(10\), podremos asignar colores y formas diferentes a los puntos localizados dentro y fuera del círculo, según corresponda, y podremos también contarlos y, en última instancia, aproximar el valor de \(\pi\).
El método anterior está programado en la función plot_points(), que creamos anteriormente.
Si establecemos que el argumento hyp es igual a TRUE, en el título de cada visualización se imprimirá la cantidad de puntos y la aproximación de \(\pi\) que resulta de implementar el método:
plot_points(n = 10, hyp = TRUE) 
Arriba miramos \(10\) puntos, de los cuales sólo uno está localizado fuera del círculo. Con un total de apenas \(10\) puntos la aproximación de \(\pi\) la podemos calcular mentalmente, siguiendo esta fórmula:
\[ 4 * \frac{número\ de\ puntos\ dentro\ del\ circulo}{número\ total\ de\ puntos} \]
\[ = 4 * \frac{9}{10} \]
\[ = 3.6 \]
\(3.6\) es justamente el valor aproximado de \(\pi\) impreso en el título de la visualización anterior.
¿Por qué multiplicamos por \(4\)? Recordemos que este método para aproximar \(\pi\) se basa en una proporción entre el área del cuadrado (el todo) y el área del círculo inscrito en él (una parte del todo).
Más exactamente, la proporción entre las áreas podemos expresarla de este modo:
\[ \frac{área\ del\ círculo}{área\ del\ cuadrado} = \frac{\pi r^2}{(2r)^2} \]
En este ejemplo, \(r = 10\), de modo que:
\[ \frac{área\ del\ círculo}{área\ del\ cuadrado} = \frac{\pi (10)^2}{(2*10)^2} \]
\[ = \frac{100\pi}{400} \]
\[ = \frac{\pi}{4} \]
Es decir:
La proporción entre el número de puntos dentro del círculo y el número total de puntos se aproxima a la proporción entre el área del círculo inscrito en el cuadrado y el área del cuadrado.
A su vez, la proporción entre el área del círculo inscrito en el cuadrado y el área del cuadrado representa una cuarta parte de \(\pi\), por lo que debemos multiplicar ese resultado por \(4\) para aproximar la totalidad de \(\pi\).
Quarto
Este documento lo preparé en Quarto. Entre otras cosas, Quarto permite combinar código y texto, imprimir visualizaciones, incluir referencias bibliográficas y escribir fórmulas y ecuaciones matemáticas. Es un gran asset aprender a preparar informes en Quarto.
Con más y más puntos notaremos que el círculo que dibujamos al principio se va formando solito; también el valor de \(\pi\) se va aproximando más a su valor real, a saber, \(3.1415927\):
p1 <- plot_points(n = 100, hyp = TRUE, cir = TRUE)
p2 <- plot_points(n = 1000, hyp = TRUE, cir = TRUE)
p3 <- plot_points(n = 10000, hyp = TRUE)
p4 <- plot_points(n = 100000, hyp = TRUE)
(p1 + p2) / (p3 + p4)
Hemos ingeniado un método que aproxima el valor de \(\pi\) a partir de datos que siguen una distribución uniforme, y lo hemos programado en la función plot_points(). 😎
Es todo, por ahora.
Contact me
Edwin Alvarado-Mena. PhD Student & Data Science Ambassador at the University of Arizona, USA. Programming literacy advocate and Certified Carpentries Instructor, he is also pursuing the Master’s degree in Data Science and the Computational Social Science Graduate Certificate.
References
Rogers, Simon, and Mark Girolami. 2016. A First Course in Machine Learning. Chapman; Hall/CRC.