R:/

Se generó una función que toma como parámetro un número n de puntos y calcula una estimación para Pi.

# Carga el paquete knitr y kable para la creación de tablas
library(knitr)
library(kableExtra)

# Inicializa el vector de resultados
results <- list()

# Función para estimar Pi
estimar_pi <- function(n) {
  
  # Generación de n puntos x e y uniformemente distribuidos entre 0 y 1
  x <- runif(n, 0, 1)
  y <- runif(n, 0, 1)
  
  # Cálculo de las distancias al centro y determinación de puntos dentro del círculo ed radio 0.5
  in_circle <- sum((x - 0.5)^2 + (y - 0.5)^2 < 0.25)
  
  # Estimación de π, dado que π/4 es aproximadamente in_circle/n
  pi_estimate <- (in_circle / n) * 4
  
  return(pi_estimate)
}

Con la función definida, la función es llamada con diferentes valores de n para observar las diferentes estimaciones de Pi producidas. Se ejecutó cada estimación 1000 veces.

# Realiza n_reps repeticiones de cada estimación
n_reps=1000
for(n in c(100, 1000, 10000, 100000, 1000000)) {
  pi_estimates <- numeric(n_reps) # Almacenar las estimaciones de π
  for(j in 1:n_reps) {
    pi_estimates[j] <- estimar_pi(n)
  }
  results[[as.character(n)]] <- pi_estimates
}

Se obtuvo el promedio para obtener una mejor idea del error medio de las estimaciones.

# Calcular la media y desviación estándar para las estimaciones
statistics <- data.frame(
  Sample_Size = character(),
  Mean_Estimate = numeric(),
  SD_Estimate = numeric()
)

for(n in names(results)) {
  statistics <- rbind(statistics, data.frame(
    Sample_Size = n,
    Mean_Estimate = mean(results[[n]]),
    SD_Estimate = sd(results[[n]])
  ))
}

# Ampliar el data.frame de estadísticas para incluir el Error del Promedio
statistics <- transform(statistics, Error_Promedio = numeric(nrow(statistics)))

for(i in 1:nrow(statistics)) {
  statistics$Error_Promedio[i] <- abs(statistics$Mean_Estimate[i] - pi)
}

# Mostrar la tabla
kable(statistics, "html", escape = FALSE, caption = "Media y Desviación Estándar de Estimaciones de π y Errores") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "bordered"), full_width = F)

Media y Desviación Estándar de Estimaciones de π y Errores
Sample_Size	Mean_Estimate	SD_Estimate	Error_Promedio
100	3.149760	0.1637228	0.0081673
1000	3.143744	0.0536972	0.0021513
10000	3.141126	0.0167177	0.0004663
1e+05	3.141520	0.0050859	0.0000727
1e+06	3.141615	0.0016618	0.0000228

Como se puede observar, la estimación se acerca cada vez mas al valor de Pi (aprox 3.141593) entre mas puntos se usen. Con n=100.000 coinciden hasta la cuarta cifra decimal.

# Convertir Sample_Size a numérico para plot
statistics$Sample_Size <- as.numeric(as.character(statistics$Sample_Size))

# Calcular los límites para el eje y
y_max <- max(statistics$Mean_Estimate + statistics$SD_Estimate)
y_min <- min(statistics$Mean_Estimate - statistics$SD_Estimate)

# Gráfico de la media de las estimaciones vs. tamaño de la muestra con límites ajustados en el eje y
plot(statistics$Sample_Size, statistics$Mean_Estimate, 
     xlab = "Tamaño de la muestra", ylab = "Media de la Estimación de π", 
     main = "Media de la Estimación de π vs. Tamaño de la Muestra", 
     pch = 19, type = "b", col = "blue", log = "x", ylim = c(y_min, y_max))

# Añadir barras de error
arrows(statistics$Sample_Size, statistics$Mean_Estimate - statistics$SD_Estimate, 
       statistics$Sample_Size, statistics$Mean_Estimate + statistics$SD_Estimate, 
       angle = 90, code = 3, length = 0.05, col = "red")

Se puede observar que la desviación estándar del promedio de las estimaciones disminuye a medida que incrementa el tamaño de la muestra. El error promedio es menor a 0.0001 desde aproximadamente 100 mil puntos (1e+05).

# Gráfico del error promedio vs. tamaño de la muestra
plot(statistics$Sample_Size, statistics$Error_Promedio, 
     xlab = "Tamaño de la muestra", ylab = "Error Promedio", 
     main = "Error Promedio vs. Tamaño de la Muestra", 
     pch = 19, type = "b", col = "red", log = "x")

# Añadir líneas de cuadrícula para facilitar la lectura
grid(nx = NULL, ny = NULL, col = "gray", lty = "dotted")

Se puede observar que el error promedio disminuye con e tamaño de muestra, con variaciones mas pequeñas a partir de 100.000 puntos.

# Vector vacío para almacenar los tiempos de ejecución de cada ejecución
e_times <- numeric(length = 8)
# Tamaños de muestra
sample_sizes <- c(10, 100, 1000, 10000, 100000, 1000000)

# Ejecución de la estimación de π para cada tamaño de muestra y registro del tiempo de ejecución
for (i in 1:length(sample_sizes)) {
  start_time <- Sys.time()
  
  # Llamada a la función estimar_pi con el tamaño de muestra actual
  estimar_pi(sample_sizes[i])
  
  end_time <- Sys.time()
  e_times[i] <- as.numeric(end_time - start_time, units="secs")
}

sample_sizes <- sample_sizes[1:length(e_times)]

# Crear el gráfico
plot(sample_sizes, e_times, type = "b", pch = 19, col = "blue",
     xlab = "Tamaño de la muestra", ylab = "Tiempo de ejecución (segundos)",
     main = "Tiempo de Ejecución vs. Tamaño de la Muestra", log = "x")

# Añadir líneas de cuadrícula para facilitar la lectura
grid(nx = NULL, ny = NULL, col = "gray", lty = "dotted")

En el gráfico, el tiempo de ejecución presenta un incremento exponencial con el incremento en potencias de 10 del tamaño de la muestra, que generan una pequeña disminución en el error de estimación. Este tiempo puede variar con la capacidad computacional usada, pero da una idea del incremento exponencial del incremento en potencias de 10 del tamaño de muestra. Esta dificultad se puede mitigar en parte, gracias a que la función estimar_pi(n) es altamente paralelizable, ya que el ciclo for de generación de puntos puede realizarse en paralelo al no requerír datos de un cálculo anterior para el siguiente. Al ser paralelizable, se puede aprovechar hardware como las GPU actuales como la NVidia RTX4090 con 16.384 núcleos CUDA, para reducir en varios ordenes de magnitud los tiempos de ejecución.

Actividad 2, Problema 1 de Métodos y Simulación Estadística

1. Estimación del valor de π

R:/