Aplicaciones de simulación con números aleatorios

Alex J. Zambrano
alexzambrano@usantotomas.edu.co

Cálculo numérico y Simulación

Visualización del Teorema Central del Límite

Sea \(X_1,X_2,\ldots,X_n\) una muestra aleatoria procedente de una población con media \(\mu\) y desviación estándar \(\sigma\). Si \(n\) es suficientemente grande, la v.a. \[ \bar{X}=\frac{1}{n}\sum_{i=1}^nX_i \] sigue aproximadamente una distribución normal con media \(\mu_{\bar{X}}=\mu\) y desviación estándar \(\sigma_{\bar{X}}=\sigma/\sqrt{n}\)

Veamos esté ejemplo en R, utilizando números aleatorios uniforme entre 0 y 1.

aleatorios <- matrix(runif(12000),ncol=12)
medias <- apply(aleatorios,1,mean)
hist(medias,freq = F)
lines(seq(from = 0,to = 1,by = 1/1000),dnorm(seq(from = 0,to = 1,by = 1/1000),mean = 1/2,sd = 1/12),lwd=2,col="sienna")

Ahora vemos si los datos se distribuyen exponencial con \(\lambda=1/5\)

aleatorios <- matrix(rexp(n = 12000,rate = 1/5),ncol=12)
medias <- apply(aleatorios,1,mean)
hist(medias,freq = F)
lines(seq(from = 0,to = 12,by = 1/1000),dnorm(seq(from = 0,to = 12,by = 1/1000),mean = 5,sd = 5/sqrt(12)),lwd=2,col="sienna")

Problema de la secretaria

Famoso problema en Estadística que se puede resumir de la siguiente manera: Una persona debe seleccionar una secretaria para llenar una vacante. Asumamos que en la población hay \(n\) aspirantes que pueden ordenarse desde la peor, hasta la mejor. Si asumimos que las secretarias llegan a la cita aleatoriamente y que una vez una secretaria sea entrevistada la desición debe tomarse inmediatamente: se le da el puesto o se le dice que no cumple con los requisitos. Para hacer un poco más realista el problema, supongamos que cada entrevisita tiene un costo (en tiempo y dinero) por lo tanto no es aconsejable entrevistar muchas aspirantes.

Consideremos el caso en el que se entrevista 3 aspirantes. La probabilidad de escoger la mejor aspirante es de 1/3. Sin embargo, no se puede realizar las 3 entrevistas, ya que cada una de ellas requiere un costo. Por lo cual se opta realizar la siguiente estrategia:

¿En cuantos casos hemos acertado con esta estrategia?

Via simulación veamos que tan acertada es la estrategia

Supongamos que la calificación se rankea de 1 a 3 siendo 3 la calificación más alta. Simulemos la estrategia anterior 10000 veces miremos que tan acertada fue la elección

ranking <- 1:3
resultados <- numeric(10000)
for(i in 1:10000){
  entrevista <- sample(x = ranking,size = 3,replace = F) # Se estable un orden para realizar las entrevistas
  sel <- ifelse(entrevista[1] < entrevista[2],2,3)  # Candidata seleccionada
  resultados[i] <- ifelse(sel == 3,1,0)  # resultados acertados
}
mean(resultados)
## [1] 0.4934

Estimación de \(\pi\)

Suponga que el vector aleatorio \((X,Y)\) está distribuido uniformemente en el cuadrado de área 4 centrado en el origen.

Considere la probabilidad de que un punto aleatorio del cuadrado esté contenido en el círculo inscrito de radio 1.

Note que \((X,Y)\) se distribuye uniformemente en el cuadrado, se tiene que: \[ \begin{align*} P\{(X,Y) \text{ esta en el circulo}\} &= P\{X^2+Y^2\leq 1\}\\ &=\frac{\text{Area del circulo}}{\text{Area del cuadrado}}\\ &=\frac{\pi}{4} \end{align*} \]

Ahora si \(X\) e \(Y\) son independientes entonces \(x\in(-1,1)\) e \(y\in(-1,1)\). Ahora si \(U\sim U(0,1)\), entonces \(2U-1\sim U(-1,1)\). Por lo tanto, si generamos números aleatorios \(U_1\) y \(U_2\), fijamos \(X=2U_1-1\) e \(Y=2U_2-1\) y definimos lo siguiente: \[ I=\begin{cases} 1 & \text{ si } X^2+Y^2\leq 1\\ 0 & \text{e. o. c.} \end{cases} \] entonces \[ E[I]=P\{X^2+Y^2\leq 1\}=\frac{\pi}{4} \] Veamos lo siguiente en R

U1 <- runif(10000)
U2 <- runif(10000)
X <- 2*U1-1
Y <- 2*U2-1
I <- ifelse(X^2+Y^2<=1,1,0)
mean(I) # estimación de pi/4
## [1] 0.7901
4*mean(I) # estimación de pi
## [1] 3.1604

Juego de dados

Suponga que lo invitan a participar en este juego: Ud. lanza un dado, si sale 1 o 6, ud. repite el lanzamiento. Si sale un impar ud. gana en pesos el cuadrado del número resultante, pero si sale par, ud. le toca pagar en pesos el cubo del resultado. ¿Cuál es la ganancia esperada de este juego?

Realizando en simulación tenemos lo siguiente

n <- 10000 # Número de experimentos a realizar
ganancia <- numeric(n)
for(i in 1:n){
  dados <- sample(1:6,size = 1,replace = F) # Seleccionamos una cara
  while(dados == 1 || dados == 6){ # Si cae 1 o 6 repetimos lanzamiento
    dados <- sample(1:6,size = 1,replace = F)
  }
  ganancia[i] <- ifelse(dados == 2 || dados ==4, -dados^3, dados^2)
}
mean(ganancia)
## [1] -9.602

Juego del dardo

Este juego consiste en lanzar un dardo sobre el tablero de tiro al blanco con 5 regiones enumeradas. Si la probabilidad de que el dardo caíga en la región \(i\) (\(i=1,2,\ldots,5\)) es proporcional al área de la región, ¿Cuál es el valor esperado del juego, si al caer el dardo en la región \(i\) se le asigna \(i\) puntos al jugador?

La probabilidad de obtener \(i\) puntos es \[ p(i)=\frac{(6-i)^2-(5-i)^2}{5^2}\quad i=1,2,\ldots,5 \] de está forma \[ \mu=\sum_{i=1}^5 i \cdot p(i)=2.2 \] Realizando un codigo de simulación en R, se tiene que:

x <- 1:5
p <- ((6-x)^2-(5-x)^2)/5^2
i <- sample(1:5,100000,replace=T,prob=p)
mean(i)
## [1] 2.20192

Ejercicios

  1. Se lanza un dado normal y luego se lanza tantas veces como haya mostrado el primer lanzamiento. Se suman los resultados de estos lanzamientos. Defina una variable aleatoria \(X\) que modele este experimento. También defina una variable aleatoria \(Y\) que cuente el número de lanzamientos. Defina y programe un algoritmo que simule este experimento.

  2. Se lanza de manera continua un par de dados normales, hasta que todos los posibles resultados \(2,3,\ldots,12\) hayan aparecido al menos una vez. Desarrolle un estudio de simulación para estimar el número esperado de lanzamientos necesarios.

  3. Utilice la simulación para aproximar \(cov(U,e^U)\), donde \(U\) es uniforme en \((0, 1)\). Compare su aproximaeión con la respuesta exacta.

  4. Sea \(U\) uniforme en \((0, 1)\). Utilice simulación para aproximar lo siguiente \(cov(U,\sqrt{1 - U^2})\).

  5. Para \(U_1,U_2,\ldots,\) variables aleatorias uniformes en \((0, 1)\), definimos \[ N= \min\left\{n\mid \sum_{i=1}^n U_i> l\right\} \] Es decir, N es igual a la cantidad de números aleatorios que deben sumarse para exceder a l. ¿Cuál cree que sea el valor de \(E[N]\)?

References

Correa, Juan Carlos. 2001. “Elementos de Simulación.”

Jalon, David. 2014. “El Problema de La Secretaria.” http://jalon.org/MAES/SECRETARIA/SECRETARYv082.pdf.

Ross, Sheldon M. 1999. Simulación. Segunda. Prentice Hall.