El cálculo de integrales es un cálculo de áreas. La idea de la integración por Monte Carlo es lanzar un gran número de puntos uniformemente distribuidos en un rectángulo $R$ y contar los que cayeron debajo de la curva. Si lanzo $N$ puntos y n cayeron por debajo, entonces:
$$ \frac{n}{N} \simeq \frac{\int f(x) dx}{A_{R}} $$y la aproximación es cada vez mejor si N es grande.
¿Cómo decidimos si los puntos que lanzamos están por debajo de la curva? Si $\vec{x}=(x,y) \in R$ es un punto lanzado al azar, $\vec{x}$ está por debajo de la curva si:
$$ y \le f(x) $$(estamos suponiendo $x,y \ge 0$).
Vamos a usar la integración Monte Carlo para calcular una aproximación de $\pi$. Si tomo el arco del círculo unitario $f(x) = \sqrt{1-x^2}$ con $x \in [0,1]$, calcular el área bajo la curva me da $\pi /4$.
Una ilustración gráfica:
In [2]:
using PyPlot
In [4]:
f(x) = sqrt(1-x^2)
Out[4]:
In [15]:
# Graficamos la función f(x)
x = linspace(0,1,256)
figure(figsize=(4,4))
plot(x, [f(a) for a in x],"k-")
# Lanzamos un número grande de puntos en el rectángulo unitario
N = 10_000
puntos = rand(N,2)
# Los clasificamos según si están abajo o arriba de la curva
abj_x = Float64[]
abj_y = Float64[]
arr_x = Float64[]
arr_y = Float64[]
for i in 1:N
x, y = puntos[i,:]
if y <= f(x)
push!(abj_x, x) ; push!(abj_y, y)
else
push!(arr_x, x) ; push!(arr_y, y)
end
end
# Los graficamos
plot(abj_x,abj_y, "r.", markersize=3.)
plot(arr_x,arr_y, "k.", markersize=3.) ;
Ahora, como estamos trabajando con el círculo unitario, hay una prueba más sencilla (usa menos recursos computacionalmente) para saber sí el punto está dentro del círculo:
$$\vec{x} \in C \,\, si \,\, x^2 + y^2 \le 1$$Vamos a escribir las funciones que calculan $\pi$ con ésta condición.
In [16]:
function norma_cuadrado(v)
x, y = v
x^2 + y^2
end
Out[16]:
In [17]:
norma_cuadrado([1,2])
Out[17]:
Incidentalmente, Julia tiene implementada una función para calcular la norma de un vector pero si no sabemos eso usamos nuestra función arriba.
In [21]:
norm([1,2])
Out[21]:
In [24]:
function aproxima_pi(N) # N es el número de puntos
puntos = rand(N,2)
contador = 0
for i in 1:N
v = puntos[i,:]
if norma_cuadrado(v) <= 1
contador += 1
end
end
4*contador/N
end
Out[24]:
En el notebook, Ctrl + Enter evalúa la celda y la mantiene seleccionada, es útil para correr muchas veces el mismo código.
In [31]:
aproxima_pi(100)
Out[31]:
Ahora nos interesa saber qué tan buena es nuestra aproximación de $\pi$ cuando N aumenta:
$$ e_{rel} = \frac{|\pi - Aprox|}{\pi}$$
In [36]:
Ns = [10:2:100] # Va de 10 hasta 100 con pasos de 2
error_relativo = Float64[]
for N in Ns
aprox = aproxima_pi(N)
push!(error_relativo, abs(aprox-pi)/pi )
end
figure(figsize=(5,3))
ylabel(L"$e_{rel}$")
xlabel(L"$N$")
plot(Ns, error_relativo, "o-") ;
No se ve gran cosa, necesitamos $N$ mucho más grande. Vamos a tomar N en una distribución logarítmica:
In [40]:
Ns = [int(10^i) for i in 1:0.1:5]
error_relativo = Float64[]
for N in Ns
aprox = aproxima_pi(N)
push!(error_relativo, abs(aprox-pi)/pi )
end
figure(figsize=(5,3))
ylabel(L"$e_{rel}$")
xlabel(L"$N$")
loglog(Ns, error_relativo, "o-") ;
Ahora sí podemos ver una tendencia. A pesar de las fluctuaciones, podemos ver que la tendencia parece ser una recta decreciente. Vamos a tomar varias corridas para limpiar la curva. Vamos a requerir bastantes cálculos y en Julia ésto es más rápido si lo metemos dentro de una función:
In [54]:
function error_aproximacion(Ns)
N = length(Ns)
errores = Array(Float64,N) # Esto construye una lista de tamaño N, que guarda flotantes
for i in 1:N
aprox = aproxima_pi(Ns[i])
errores[i] = abs(aprox-pi)/pi
end
errores
end
Out[54]:
Hace lo mismo que lo que teníamos arriba:
In [55]:
Ns = [int(10^i) for i in 1:0.1:5]
figure(figsize=(5,3))
ylabel(L"$e_{rel}$")
xlabel(L"$N$")
loglog(Ns, error_aproximacion(Ns), "o-") ;
Ahora tomamos el promedio:
In [61]:
function error_promedio(Ns, n_promedio)
N = length(Ns)
# Vamos a almacenar cada cálculo en las columnas de una matriz
corridas = Array(Float64,(N, n_promedio))
for i in 1:n_promedio
corridas[:,i] = error_aproximacion(Ns)
end
#Calculamos el valor promedio: Julia tiene una función para ésto
out = Array(Float64, N) # Va contener los promedios
for j in 1:N
out[j] = mean(corridas[j,:])
end
out
end
Out[61]:
In [70]:
Ns = [int(10^i) for i in 1:0.1:5]
n_promedio = 50
figure(figsize=(5,3))
ylabel(L"$\langle e \rangle$")
xlabel(L"$N$")
loglog(Ns, error_promedio(Ns,n_promedio), "o-") ;
# Ajustamos una recta
loglog(Ns, [0.5/sqrt(N) for N in Ns], "r--")
Out[70]:
Ahora sí vemos que el comportamiento en la gráfica log-log es una recta con pendiente aproximadamente igual a $1/2$. Una recta en log-log se llama una ley de potencias, y el comportamiento:
$$ e \sim \frac{1}{\sqrt{N}} $$es típico de la integración Monte Carlo.
La gran utilidad de la integración Monte Carlo, es que mantiene la misma convergencia aún para integrales en dimensiones mayores.
In [ ]: