1.2 Aproximación Monte Carlo
Cuando trabajamos con variables aleatorias, es común el problema de querer estimar el valor esperado de cualquier cantidad \(h(X)\), siendo \(X\) una variable aleatoria. Utilizando la simulación es posible obtener estimaciones de dichos valores de un modo relativamente sencillo: a través de la integración Monte Carlo.
Ante un problema genérico relativo a calcular el valor esperado de cierta función \(h(X)\) para una variable aleatoria con fdp \(X \sim f(x)\), \[\begin{equation} E[h(X)]=\int_S h(x) f(x) dx (\#eq:problemaestimacionMC) \end{equation}\] donde \(S\) denota el conjunto en el que la variable \(X\) toma valores,
Definición 1.8 Ante el problema de estimar (??), el procedimiento de estimación Monte Carlo propone simular una muestra aleatoria de la distribución de \(X\), \(x_1,\ldots,x_n\) y con ella obtener una aproximación empírica a través del promedio de las cantidades \(h(x_1),\ldots,h(x_n)\), \[\begin{equation} \bar{h_n}=\frac{\sum_{i=1}^n h(x_i)}{n}. (\#eq:estimacionMC) \end{equation}\]
Por la Ley de los Grandes Números, \(\bar{h_n}\) converge casi seguro a la cantidad de interés \(E[h(X)]\). Además, cuando \(h^2(X)\) tiene un valor esperado finito, la velocidad de convergencia de \(\bar{h_n}\) se puede calcular y la varianza asintótica de la aproximación es \[Var(\bar{h_n})=\frac{1}{n} \int_S (h(x)-E[h(X)])^2 f(x) dx,\] que se puede estimar con la muestra \(x_1,\ldots,x_n\) a través de \[\begin{equation} v_n=\frac{1}{n^2} \sum_{i=1}^n [h(x_i)-\bar{h_n}]^2. (\#eq:varianzaMC) \end{equation}\]
Más específicamente, por el Teorema Central del Límite, para \(n\) grande tendremos que \[\begin{equation} \frac{\bar{h_n}-E[h(X)]}{\sqrt{v_n}} \sim N(0,1), (\#eq:distMC) \end{equation}\]
lo que permitirá además, construir bandas de confianza para la aproximación Monte Carlo, a la que en adelante nos referiremos por “aproximación MC.”
\[\begin{equation} IC_{1-\alpha}[\bar{h_n}]=[\bar{h_n}-z_{1-\alpha/2} \sqrt{v_n}, \bar{h_n}+z_{1-\alpha/2} \sqrt{v_n}] \tag{1.1} \end{equation}\] donde \(z_{1-\alpha/2}\) es el cuantil \(1-\alpha/2\) de una normal estándar, asociado al nivel de confianza \(1-\alpha\).
1.2.1 MC y probabilidad
Aunque en las situaciones más sencillas se puede evaluar cualquier probabilidad mediante la correspondiente función de distribución, en muchas ocasiones resulta sencillo obtener una muestra simulada de la variable aleatoria y aproximar dicha probabilidad de interés mediante el denominado estimador Monte-Carlo.
Se obtiene la estimación Monte-Carlo (MC) de la probabilidad de que una variable aleatoria \(X\) tome algún valor en un conjunto \(A\) de valores dentro de su espacio probabilístico, \(Pr(X \in A)\), a partir de un conjunto de observaciones o simulaciones \(x_1, x_2,...,x_N\) de la variable \(X\), mediante el cociente entre el número de simulaciones que están en dicha región \(A\), y el número de datos disponibles,
\[\begin{equation} Pr(X \in A) \approx \frac{\#\{x_1, x_2,...,x_N\} \in A}{N}. \end{equation}\]Básicamente esta definición se fundamenta en la interpretación empírica de la probabilidad, a través del ratio de casos favorables por casos posibles: \[probabilidad=\frac{\mbox{casos favorables}}{\mbox{casos posibles}}.\]
La fórmula anterior la podemos expresar en notación similar a la utilizada en (??) definiendo una variable dicotómica \[I_A(X)=1 \text{, si } X \in A\] y 0 en otro caso, de manera que el problema de calcular una probabilidad según una distribución de probabilidad para \(X\) se convierte en el valor esperado de una distribución Bernouilli para una distribución de \[Pr(X \in A)=\int_A f(x)dx=\int_S I_A(X)f(x)dx=E[I_A(X)]\] y la estimación MC (??) se puede expresar, para un conjunto de \(n\) simulaciones \(x_1,\ldots,x_n\), como \[\begin{equation} Pr(X \in A)\approx \hat{h_n}= \frac{\sum_{i=1}^n I_A(x_i)}{n}, (\#eq:estimacionMCprob) \end{equation}\] y su varianza con \[\begin{equation} v_n=\frac{1}{n^2} \sum_{i=1}^n [I_A(x_i)-\bar{h_n}]^2. (\#eq:varMCprob) \end{equation}\]
Ejemplo 1.3 En la situación del ejemplo 1.1, supongamos que disponemos de un conjunto de 200 simulaciones del número de huevos defectuosos en 200 cajas distribuídas. Queremos descubrir, utilizando exclusivamente esas simulaciones:
- La probabilidad de que una caja tenga más de 3 huevos defectuosos, \(Pr(X>3)\).
- La probabilidad de que una caja tenga a lo sumo 3 huevos defectuosos, \(Pr(X\leq 3)\).
- La probabilidad de que una caja no tenga ningún huevo defectuoso, \(Pr(X=0)\).
- La probabilidad de que la proporción de huevos defectuosos en una caja sea inferior al 1%, \(Pr(X/144 < 0.01)=Pr(X<1.44)\).
# Simulaciones disponibles
defectos <- c(2, 2, 0, 0, 0, 2, 1, 2, 1, 4, 1, 0, 0, 2, 4,
0, 0, 0, 0, 1, 1, 1, 2, 2, 3, 1, 0, 4, 3, 1, 0,
2, 2, 2, 3, 1, 0, 2, 2, 2, 3, 1, 0, 1, 0, 1, 2,
0, 0, 2, 3, 2, 3, 2, 4, 4, 0, 1, 1, 3, 0, 0, 3,
2, 0, 0, 0, 3, 0, 1, 4, 1, 1, 2, 1, 1, 4, 1, 1,
1, 0, 1, 0, 1, 2, 2, 1, 3, 1, 2, 1, 2, 3, 1, 2,
5, 1, 1, 1, 1, 0, 1, 1, 1, 2, 1, 0, 0, 1, 2, 2,
1, 1, 1, 1, 0, 3, 1, 1, 1, 1, 4, 4, 0, 6, 6, 1,
1, 1, 0, 2, 3, 1, 0, 0, 2, 0, 2, 1, 1, 1, 2, 1,
1, 1, 1, 2, 5, 0, 1, 3, 1, 1, 4, 1, 2, 1, 1, 0,
2, 1, 2, 1, 3, 3, 2, 0, 3, 0, 1, 3, 0, 1, 2, 0,
1, 0, 0, 2, 2, 1, 2, 0, 0, 0, 1, 1, 2, 3, 1, 0,
1, 0, 1, 1, 1, 1, 1, 5, 3)
# Número de simulaciones/observaciones
nsim=length(defectos)
# Tamaño de la caja
tamaño <- rep(144,200)
# Conjunto de datos
huevos <- data.frame(tamaño, defectos)Y vamos aproximando por Monte-Carlo las probabilidades de interés, estableciendo las condiciones lógicas en cada situación.
# Pr(X > 3)
sel <- dplyr::filter(huevos, defectos >3)
prob <- nrow(sel)/nsim
cat("Probabilidad estimada [Pr(X > 3)]: ", prob)## Probabilidad estimada [Pr(X > 3)]: 0.075# Pr(X <= 3)
sel <- dplyr::filter(huevos, defectos <= 3)
prob <- nrow(sel)/nsim
cat("Probabilidad estimada [Pr(X <= 3)]: ", prob)## Probabilidad estimada [Pr(X <= 3)]: 0.925# Pr(X = 0)
sel <- dplyr::filter(huevos, defectos == 0)
prob <- nrow(sel)/nsim
cat("Probabilidad estimada [Pr(X = 0): ", prob)## Probabilidad estimada [Pr(X = 0): 0.235# Pr(X/144 <=0.01)
sel <- dplyr::filter(huevos, defectos <= (0.01*144))
prob <- nrow(sel)/nsim
cat("Probabilidad estimada [Pr(X/144 <=0.01): ", prob)## Probabilidad estimada [Pr(X/144 <=0.01): 0.62El error asociado a la estimación de la probabilidad anterior, lo calculamos aplicando (??), y también el intervalo de confianza al 95% con (1.1).
I.a=(huevos$defectos <= 1.44)*1
prob=mean(I.a)
cat("prob.estim=",prob)## prob.estim= 0.62vn=sum((I.a-prob)^2)/(nsim^2)
error=sqrt(vn)
cat("Error Estimado=",round(error,3))## Error Estimado= 0.034# límites del IC redondeados a 3 cifras decimales
ic.low=round(prob-qnorm(0.975)*error,3)
ic.up=round(prob+qnorm(0.975)*error,3)
cat("IC(95%)[AproxMC(prob.estim)]=[",ic.low,",",ic.up,"]")## IC(95%)[AproxMC(prob.estim)]=[ 0.553 , 0.687 ]De esta forma resultará posible estimar cualquier probabilidad asociada a una variable aleatoria, aun desconociendo su función de distribución o de probabilidad, siempre que dispongamos de una muestra simulada, y esta estimación será más precisa cuanto mayor sea el tamaño de dicha muestra.
1.2.2 MC y momentos
Si disponemos de una muestra (de observaciones o simulaciones) lo suficientemente grande de la variable aleatoria \(X\) podemos aproximar -de modo razonablemente preciso- el valor esperado y la varianza (o desviación típica) sin más que calcular la media y varianza de los datos que componen la muestra.
La precisión de estas estimaciones estará directamente relacionada con el tamaño de la muestra utilizada.
Sean \(x_1, x_2,...,x_N\) simulaciones disponibles para \(X\). Entonces \[\begin{eqnarray*} E(X) &\approx& \bar{x}_N=\sum_{i=1}^N x_i /N \\ Var(X) &\approx& \sum_{i=1}^N (x_i-\bar{x}_N)^2/n = \bar{x}_N^2-\sum_{i=1}^N x_i^2 /N. \end{eqnarray*}\]
Ejemplo 1.4 Con los datos del ejemplo anterior, queremos saber cuál es el número aproximado de huevos que se rompen en cada caja, conocer su dispersión y tener así mismo un intervalo de confianza para la media.
# media
media=mean(huevos$defectos)
# dispersión
varianza=var(huevos$defectos)
desvtip=sd(huevos$defectos)
# ic para la media
error=sqrt(sum((huevos$defectos-media)^2)/(nsim^2))
cat("\n Error Estimado (media)=",round(error,3))##
## Error Estimado (media)= 0.089# límites del IC redondeados a 3 cifras decimales
ic.low=round(media-qnorm(0.975)*error,3)
ic.up=round(media+qnorm(0.975)*error,3)
cat("IC(95%)[AproxMC(media)]=[",ic.low,",",ic.up,"]")## IC(95%)[AproxMC(media)]=[ 1.255 , 1.605 ]