El proceso de simulación constituye una herramienta poderosa para la estadística que se pueden emplear para entender relaciones complejas y estimar valores difíciles de calcular directamente. Para entenderlo utilizaremos se plantean los siguientes problemas:

0.1 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:

0.1.1 Paso 1

Genere \(n\) coordenadas \(x:X_1,...,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).

Para este punto se utilizará La función \(runif\) permite obtener \(n\) observaciones aleatorias de una distribución uniforme entre el intervalo (0,1). Adicionalmente, se probarán diferentes valores para \(n\).

par(mfrow = c(1, 3))

x <- seq(-0.5, 1.5, 0.01)

set.seed(1)

# n = 1000
hist(runif(1000), main = "n = 1000", xlim = c(-0.2, 1.25),
     xlab = "", prob = TRUE)
lines(x, dunif(x), col = "red", lwd = 2)

# n = 10000
hist(runif(10000), main = "n = 10000", xlim = c(-0.2, 1.25),
     xlab = "", prob = TRUE)
lines(x, dunif(x), col = "red", lwd = 2)

# n = 100000
hist(runif(100000), main = "n = 100000", xlim = c(-0.2, 1.25),
     xlab = "", prob = TRUE)
lines(x, dunif(x), col = "red", lwd = 2)

0.1.2 Paso 2

Genere 1000 coordenadas \(y:Y_1,...,Y_n\), utilizando nuevamente la distribución uniforme con valor mínimo de 0 y valor máximo de 1

Para este punto se utilizará La función \(runif\) permite obtener \(n\) observaciones aleatorias de una distribución uniforme entre el intervalo (0,1).

par(mfrow = c(1, 3))

set.seed(1)

y=runif(1000, min = 0, max = 1)

# n = 1000
hist(y, main = "n = 1000", xlim = c(-0.2, 1.25),
     xlab = "", prob = TRUE)

0.1.3 Paso 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\).

Para reolver este punto se gráfican cada uno de los puntos de \((X_i, Y_i)\) dentro de un circulo con radio \(0.5\).

y=runif(1000, min = 0, max = 1)
x=runif(1000, min = 0, max = 1)

rad     <- 0.5   # Valor del radio
xcenter <- 0.5  # Coordenada en x del centro
ycenter <- 0.5   # Coordenada en y del centro

theta <- seq(0, 2 * pi, length = 100)
plot(x,y, lines(x=rad * cos(theta) + xcenter,
                   y=rad * sin(theta) + ycenter, type = "l", lwd = 3, col="red"), lwd=1, col="black")

Gráficamente se observa que varios de los puntos generados para \((X_i, Y_i)\) se encuentran fuera de un circulo con radio \(0.5\). Ahora, se procede a calcular el valor \((X_i-0.5)^{2} + (Y_i-0.5)^{2}\), que es el cuadrado de la distancia y determinar si es menor que \(0.25\).

Gráficamente tenemos que:

y=runif(1000, min = 0, max = 1)
x=runif(1000, min = 0, max = 1)

x_1=1:1000

y_1=1:1000
y_1=(y_1/y_1)*0.25


Distancia=function(x,y){((x-0.5)^2)+((y-0.5)^2)}

plot(Distancia(x,y), lwd=1, col="black")
lines(x_1,y_1, type = "l", col="red")

El número de puntos por debajo de 0.25 es 746

p=Distancia(x,y)
sum(p<= 0.25)
## [1] 746

0.1.4 Paso 4

¿Cuántos de los puntos están dentro del círculo? ¿Cuál es su estimación de \(\pi\)?

El número de puntos dentro del círculo, es 746, como se presentó en el punto anterior. Dado que se trabajó con 1000 repeticiones, se divide el número de puntos dentro del círculo entre el total de puntos, obteniendo \(\frac{\pi}{4} = 0.746\). Por tanto, \(\pi = 4(0.746)\) y \(\pi = 2.98\)