Problema 1

Actividad 2. Métodos y Simulación

Problema 1

Estimación del valor de \(\pi\)

La siguiente figura sugiere como estimar el valor de \(\pi\) con una simulación.

En la figura, un circuito con un área igual a \(\frac{\pi}{4}\), está inscrito en un cuadrado cuya área es igual a 1. Se elige de forma aleatoria n puntos dentro del cuadrado. La probabilidad de que un punto esté dentro del círculo es igual a la fracción del área del cuadrado que abarca a este, la cual es \(\frac{\pi}{4}\). Por tanto, se puede estimar el valor de \(\frac{\pi}{4}\) al contar el número de puntos dentro del círculo, para obtener la estimación de \(\frac{\pi}{4}\). De este último resultado se encontrar una aproximación para el valor de \(\pi\).

Pasos sugeridos:

1. Genere n coordenadas \(x: x_1, \dots, x_n\). Utilice la distribución uniforme con valor mínimo de 0 y valor máximo de 1. La distribución uniforme genera variables aleatorias que tienen la misma probabilidad de venir de cualquier parte del intervalo \((0,1)\).
2. Genere 1000 coordenadas \(y: y_1,\dots, y_n\), utilizando nuevamente la distribución uniforme con valor mínimo de 0 y valor máximo de 1.
3. Cada punto \((x_i\),\(y_i)\) se encuentra dentro del círculo si su distancia desde el centro \((0.5,0.5)\) es menor a 0.5. Para cada par \((x_i\),\(y_i)\) determine si la distancia desde el centro es menor a 0.5. Esto último se puede realizar al calcular el valor \((x_i−0.5)^2+(y_i−0.5)^2\), que es el cuadrado de la distancia, y al determinar si es menor que 0.25.
4. ¿Cuántos de los puntos están dentro del círculo? ¿Cuál es su estimación de \(\pi\)?

Nota

• Con sólo 1000 puntos, es probable que la estimación presente un error de 0.05 o más. Una simulación con 10000 y 100000 puntos tiene mayores probabilidades de dar como resultado una estimación muy cercana al valor verdadero.
• Funciones recomendadas : runif(), funtion(){}
• Entregable : enlace en RPubs con informe 1
• Problema tomado de Navidi(2006)

Desarrollo

Se nos presentan 2 figuras geométricas concéntricas, un cuadrado rectángulo y un círculo. El área del cuadrado es base x altura, y la del círculo es \(\pi r^2\). El radio mide 0,5, luego cada lado del cuadrado es 2 x 0,5, es decir 1. Las áreas serían: * Área Círculo: \(\pi r^2 = \pi(0,5)^2 = 0,7854\) * Área Cuadrado: base x altura \(= (2r)*(2r) = 4r^2 = 4 * 0,5^2 = 4 * 0,25 = 1\)

Como se quiere calcular la probabilidad de que un punto esté dentro del circulo (siendo la superficie del cuadrado la superficie en donde se encuentran todos los puntos), esta sería la razón del área del círculo y del área del cuadrado. Entonces,

• Área círculo / Área cuadrado = \(\frac{\pi r^2}{4r^2} = \frac{\pi}{4}\)

Se realizarán los pasos mencionados en el siguiente código:

# Function to estimate the value of π
estimation_of_pi <- function(iterations = 1000, seed = 689){
  set.seed(seed)
  # [iterations] values to each coordinate (x,y) in a uniform distribution for numbers between 0 and 1 (Inside the area of the square)
  # step 1 and 2
  x <- runif(iterations, min = 0, max = 1) 
  y <- runif(iterations, min = 0, max = 1)
  # step 3. Separation from the center of the circle (0.5,0.5)
  point <- (x - 0.5)^2 + (y - 0.5)^2
  # The point is inside the circle if _point_ is smaller than 0.25
  reference = 0.25
  inside_the_circle <- which(point <= reference) 
  # step 4. Number of points within the circle
  number_of_points <- length(inside_the_circle)
  # Estimation of π is 4 * inside_the_circle
  estimated_pi <- 4 * number_of_points / iterations 
  
  #cat("Percentage of points inside the circle:", number_of_points/iterations*100,"%")
  
  return(estimated_pi)
}

Revisión del método

Se sabe que \(\pi\) = 3.1415926535. Probemos la lógica utilizada con mil iteraciones (El valor establecido por defecto) para ver cuál es el valor estimado.

estimation_of_pi()

## [1] 3.172

Solo hay un decimal correcto. Todavía estamos lejos. Se sabe que para conjuntos de datos más numerosos, mejorará la estimación, entonces ejecutamos la función para valores más grandes:

number_of_iterations <- c(123456, 5123456, 99999999, 55999999)
res <- sapply(number_of_iterations, function(n)estimation_of_pi(iterations = n))
names(res) <- number_of_iterations
res

##   123456  5123456 99999999 55999999 
## 3.138381 3.141848 3.141654 3.141591

Tal como se puede ver en la tabla, el valor estimado de π se acerca al valor real entre más iteraciones se realicen.

Gráficamente se puede representar el cambio del valor de π pintándolo con cada número de iteraciones. Para esto, al ejecutar la función, debemos guardar el valor estimado de \(\pi\), para después poder usar todos sus valores estimados en una sola gráfica.

# Function to estimate the value of π
estimation_of_pi <- function(iterations = 1000, seed = 689){
  set.seed(seed)
  # store the values in a data frame
  df <- data.frame(x <- runif(n = iterations, min = 0, max = 1),
                   y <- runif(n = iterations, min = 0, max = 1))
  # store the number of iterations in the data frame
  df$iteration <- 1:iterations
  
  # count the number of points inside the circle
  reference = 0.25
  df$points_inside_circle <- ifelse (((x - 0.5)^2 + (y - 0.5)^2) <= reference, 1, 0)
  
  # Estimation of π is 4 * inside_the_circle
  df$estimated_pi <- 4 * cumsum(df$points_inside_circle) / df$iteration 
  
  return(df)
}

Con esta función, podemos generar los valores estimados de \(\pi\) de los primeros 100 mil puntos

pi_data <- estimation_of_pi(iterations = 100000)

# plot showing estimated pi as number of points increases
ggplot(pi_data, aes(x = iteration, y = estimated_pi)) +
            geom_line(col = "blue") +
            geom_hline(yintercept = pi, col = "darkgreen") +
            ylim(c(3, 3.3)) +
            labs(title = expression(paste("Aproximación de ", pi)),
                 x = "Cantidad de puntos",
                 y = expression(paste("Valor estimado de ",pi)))

## Warning: Removed 7 rows containing missing values (`geom_line()`).

Podemos ver que en la medida en que hay mas puntos, el valor estimado de π se acerca a su valor real!