Conceptos básicos de probabilidad

Recordemos algunos conceptos de probabilidad. Sea \(X\) una variable aleatoria continua, si su función de densidad es \(f(x)\), en un intervalo \([\alpha,\beta]\), la probabilidad de que la variable aleatoria tome el valor a \(a\in[\alpha,\beta]\) está dado por: \[ P[X\leq a] = \int_\alpha^a f(x)\,dx \] La probabilidad es el área bajo la curva \(y = f(x)\) del punto \(\alpha\) al punto \(a\).

El cálculo de la esperanza de una variable aleatoria \(X\) continua con función de densidad \(f(x)\) se define por \[ E(x) = \int_\alpha^\beta xf(x)\,dx \]

De igual forma se puede definir la esperanza de una función continua de una variable aleatoria. Sea \(g : \Re \longrightarrow \Re\) continua en el intervalo \([\alpha, \beta]\) entonces \[ E[g(x)] = \int_\alpha^\beta g(x)\,f(x)\,dx \]

Definamos la función de densidad de una distribución uniforme en el intervalo \((0, 1)\) dada por: \[ f(x) =\begin{cases} 1 & 0<x<1\\ 0 & \text{e.o.c.} \end{cases} \] Por lo cual definimos la esperanza de una función continua de una variable aleatoria continua uniforme de la siguiente manera \[ E[g(x)] = \int_0^1 g(x)\,dx \] Como la esperanza de una variable aleatoria continua uniforme es una integral, la media muestral se puede usar para estimar el valor de una integral.

Integración por aproximación de Monte Carlo

Si \(U_1,U_2,\ldots,U_k\) son una muestra de v.a. independientes uniformes \((0,1)\), entonces las v.a. \(g(U_1),g(U_2),\ldots,g(U_k)\) son i.i.d. con media \(\theta\). Entonces, por la ley de los grandes números \[ \sum_{i=1}^k\frac{g(U_i)}{k} \longrightarrow E[g(x)]=\theta \quad\text{si}\quad k\longrightarrow \infty \] Está aproximación de las integrales es llamada aproximación de Monte Carlo.

Ejemplo 1

Realice la siguiente integral médiante el método montecarlo \[ \int_0^1e^{-x^2}\,dx \] Para ello debemos generar una función

g <- function(x){
  exp(-x^2)
}

Gráficamente se desea integar la siguiente función

curve(g,xlab="Function",ylab="",lwd=2)

Para ello entonces generamos 10000 números aleatorios entre \((0,1)\)

set.seed(20150907)
U<-runif(10000)
head(U)

## [1] 0.04589137 0.92775662 0.04946757 0.36739956 0.49676768 0.67351615

Finalmente evaluamos cada valor de U entre g y promedia todos los resultados

mean(g(U))

## [1] 0.7461881

A continuación se describe gráficamente como se aproxima la integral mediante los números aleatorios generados.

n<-1:10000
I<-cumsum(g(U))/n
plot(n,I,type="l")

Utilizando R se realiza de la siguiente manera

integrate(g,lower = 0,upper = 1)

## 0.7468241 with absolute error < 8.3e-15

Si deseamos computar \[ \int_a^b g(x)\,dx \] Sustituyendo \(u=\frac{x-a}{b-a}\) cuando \(x\longrightarrow a\) entonces \(u\longrightarrow 0\); cuando \(x\longrightarrow b\) entonces \(u\longrightarrow 1\).

Derivando tenemos que \(du=\frac{dx}{b-a}\), por lo cual \(x=(b-a)\,u+a\), y \(dx=(b-a)\,du\).

Reemplazando tenemos que \[ \int_0^1 g((b-a)u+a)\,(b-a)\,du=\int_0^1 h(u)\,du \] donde \(h(u)=(b-a)\,g(a+(b-a)u)\)

Ejemplo 2

Realice la siguiente integral médiante el método montecarlo \[ \int_{-2}^2[\cos(50x) + \sin(20x)]^2\,dx \] Para ello debemos generar una función

g <- function(x){
  (cos(50*x)+sin(20*x))^2
}

Gráficamente se desea integar la siguiente función

curve(g,from= -2,to = 2,xlab="Function",ylab="",lwd=2)

Para ello entonces generamos 10000 números aleatorios entre \((0,1)\)

set.seed(20150907)
U<-runif(10000)
head(U)

## [1] 0.04589137 0.92775662 0.04946757 0.36739956 0.49676768 0.67351615

Finalmente evaluamos cada valor de \(U\) entre \(h=g((b-a)U+a)(b-a)\) y se promedia todos los resultados

a <- -2
b <- 2
h <- g((b-a)*U+a)*(b-a)
mean(h)

## [1] 4.025422

A continuación se describe gráficamente como se aproxima la integral mediante los números aleatorios generados.

n<-1:10000
I<-cumsum(g((b-a)*U+a)*(b-a))/n
plot(n,I,type="l")

Utilizando R se realiza de la siguiente manera

integrate(g,lower = -2,upper = 2)

## 4.016114 with absolute error < 5.8e-10

Integración multivarida por aproximación de Monte Carlo

Supogase que \(g\) es una función con argumento \(n\)-dimensionaly estamos interesados en el computo en el computo de \[ \theta=\int_0^1\int_0^1\cdots\int_0^1g(x_1,x_2,\ldots,x_n)\,dx_1\,dx_2\cdots\,dx_n \] La clave de la aproximación de Monte Carlo para estimar \(\theta\) puede expresarse como la siguiente esperanza \[ \theta=E[g(U_1,U_2,\ldots,U_n)] \] donde \(U_1,U_2,\ldots,U_n\) son v.a. uniformes independientes en \((0,1)\). Por lo tanto, generando \(k\) conjuntos independientes, cada uno de \(n\) v.a. independientes uniformes \((0,1)\) \[ \begin{matrix} U_{11}& U_{12}& \ldots & U_{1n}\\ U_{21}& U_{22}& \ldots & U_{2n}\\ \vdots& \vdots & \ddots & \vdots\\ U_{k1}& U_{k2}& \ldots & U_{kn}\\ \end{matrix} \] Entonces, des que las v.a. \(g(U_{i1},U_{i2},\ldots,U_{in})\) para \(i=1,\ldots,k\), sean todas v.a. i.i.d. con media \(\theta\), podemos estimar \(\theta\) como \[ \frac{\sum_{i=1}^kg(U_{i1},U_{i2},\ldots,U_{in})}{k} \]

Ejemplo 3

Realice la siguiente integral médiante el método montecarlo \[ \int_{0}^1\int_{0}^1(x_1^2+x_2^2)\,dx_2\,dx_1 \] Supongamos que \(X=[x_1,x_2]\), por lo cual la función a tener en cuenta es la siguiente

g <- function(X){
  X[1]^2+X[2]^2
}

Para ello entonces generamos 10000 valores de \(U=[u_1,u_2]\), t.q. \(U_1,U_2\) se distribuyen uniformemente entre \((0,1)\).

set.seed(20150907)
U1<-runif(10000)
U2<-runif(10000)
U<-cbind(U1,U2)
head(U)

##              U1        U2
## [1,] 0.04589137 0.4499325
## [2,] 0.92775662 0.8019573
## [3,] 0.04946757 0.3215803
## [4,] 0.36739956 0.2402872
## [5,] 0.49676768 0.3376965
## [6,] 0.67351615 0.7855377

Evaluamos cada valor de \(U\) en la función g y promediamos los resultados

I <- apply(U,1,g)
mean(I)

## [1] 0.6660147

A continuación se describe gráficamente como se aproxima la integral mediante los números aleatorios generados.

n<-1:10000
theta<-cumsum(I)/n
plot(n,theta,type="l")

Ejercicios

Emplee la simulación para aproximar las siguientes integrales. Compare su estimación con la respuesta exacta, si ésta se conoce.

1. \(\int_0^1\exp(\exp(x))\,dx\).
2. \(\int_0^1(1-x^2)^{3/2}\,dx\).
3. \(\int_{-2}^2\exp{(x+x^2)}\,dx\).
4. \(\int_0^\infty x(1+x^2)^{-2}\,dx\).
5. \(\int_{-\infty}^\infty\exp(-x^2)\,dx\).
6. \(\int_0^1\int_0^1\exp(-(x^2+y^2))\,dy\,dx\).
7. \(\int_0^\infty\int_0^x\exp(-(x+y))\,dy\,dx\).