<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Integración numérica Montecarlo
Referencia
Step1:  Integración Montecarlo tipo 1
Se basa en la definición de valor promedio de una función y en el valor esperado de una variable aleatoria uniforme.
Presentamos esto mediante un ejemplo.
Ejemplo. Aproxime el área bajo la curva $y=x^2$ en el intervalo $\left[0,1\right]$.
Veamos primero cómo luce dicha área.
Step2:  Entonces, lo que queremos es aproximar el área de la región $\mathcal{D}$. Llamaremos esta área $A(\mathcal{D})$.
De cálculo integral, sabemos que
$$A(\mathcal{D})=\int_{0}^{1}y\text{d}x=\int_{0}^{1}x^2\text{d}x$$.
Por definición, el valor promedio de una función $f
Step3:  En este caso, la integral se puede hacer fácilmente. Comparemos el resultado con el valor real
Step4:  Ver que los resultados son distintos cada vez (¿porqué?). Sin embargo, se aproximan más o menos en la misma medida.
Aproximación de integrales en intervalos distintos a $\left[0,1\right]$.
Sin embargo, no todas las integrales que hacemos son en el intervalo $\left[0,1\right]$. En general, podemos integrar cualquier función continua en el intervalo $\left[a,b\right]$, donde $a,b\in\mathbb{R}$ con $a<b$.
Sea $f
Step5:  Actividad. Utilizar la anterior función para realizar las siguientes integrales. Poner los resultados en una tabla cuyas filas correspondan a la cantidad de términos utilizados en la aproximación (usar 10, 100, 1000, 10000 y 100000 términos) y cuyas columnas correspondan a las funciones.
- $\int_{4}^{5} e^{x^2}\text{d}x$.
- $\int_{4}^{5} \frac{1}{log(x)}\text{d}x$.
- $\int_{4}^{5} \frac{sin(x)}{x}\text{d}x$.
Step6:  Integración Montecarlo tipo 2
Con la integración montecarlo tipo 1 pudimos aproximar integrales de funciones continuas de una variable en un intervalo dado. En realidad este mismo análisis se puede ampliar para aproximar integrales definidas de funciones continuas de varias variables (integrales sobre áreas, volúmenes e hipervolúmenes) dado que la noción de valor promedio de una función se extiende a cualquier dimensión.
Este es en realidad el caso interesante, pues las integrales de funciones complicadas también se pueden aproximar por métodos numéricos clásicos, pero cuando la dimensión aumenta es cuando montecarlo se vuelve una herramienta relevante. Dado que no lo veremos en clase por la limitación de que la mayoría no han visto cálculo en varias variables, este tema puede ser elegido como proyecto de módulo, donde se exploraría también como mejorar la aproximación de integrales montecarlo.
Como vimos en el ejemplo (y como debe ser claro de su curso de cálculo integral) una de las aplicaciones más importantes de la integración es hallar áreas. Y no solo el área bajo una curva, sino áreas entre curvas y áreas de regiones más complicadas.
Antes de ver la integración montecarlo tipo 2, ¿cómo podemos usar la integración montecarlo tipo 1 para aproximar el área entre curvas?
Ejemplo. Aproxime el área entre las curvas $y=x$, y $y=x^2$ en el intervalo $\left[0,1\right]$.
Veamos primero cómo luce dicha área.
Step7:  De cálculo integral, sabemos que
$$A(\mathcal{D})=\int_{0}^{1}x-x^2\text{d}x.$$
Entonces...
Step8:  De modo que si la región se puede describir fácilmente, diría el ferras 'no hay pedo, lo pago' (podemos usar montecarlo tipo 1).
Step9:  Pero, ¿qué pasa si la geometría de la región no se puede describir fácilmente?
Como en el caso anterior, motivaremos el método con un caso conocido. Vamos a aproximar el valor de $\pi$ usando el área de un círculo unitario.
Dibujemos el círculo unitario en la región $\mathcal{R}=\left[-1,1\right]\times\left[-1,1\right]$.
Step10:  Si aproximamos $A(\mathcal{D})$ aproximamos el valor de $\pi$, pues el área del círculo unitario es
Step11:  La probabilidad de que el punto $(X,Y)$ esté en el círculo unitario $\mathcal{D}$ es
$$P((X,Y)\in\mathcal{D})=\frac{A(\mathcal{D})}{A(\mathcal{R})}=\frac{\pi}{4}.$$
Luego, definimos una variable aleatoria de Bernoulli $B$ de manera que
$$B=\left\lbrace\begin{array}{ccc}0 & \text{si} & (X,Y)\notin\mathcal{D}\1 & \text{si} & (X,Y)\in\mathcal{D} \end{array}\right.=\left\lbrace\begin{array}{ccc}0 & \text{si} & X^2+Y^2>1\1 & \text{si} & X^2+Y^2\leq 1 \end{array}\right..$$
Entonces, el valor esperado de la variable aleatoria $B$ es
$$E\left[B\right]=\theta=P((X,Y)\in\mathcal{D})=\frac{A(\mathcal{D})}{A(\mathcal{R})}.$$
De lo anterior, una estimación de theta se puede obtener como
$$\theta=\frac{A(\mathcal{D})}{A(\mathcal{R})}\approx \frac{1}{N}\sum_{i=1}^{N}b_i,$$
donde 
$$b_i=\left\lbrace\begin{array}{ccc}0 & \text{si} & x_i^2+y_i^2>1\1 & \text{si} & x_i^2+y_i^2\leq 1 \end{array}\right.$$
son realizaciones de la variable aleatoria $B$, que a su vez es producto de las realizaciones $x_i$ e $y_i$ de las variables aleatorias $X$ e $Y$, respectivamente.
Finalmente, la aproximación montecarlo tipo 2 con $N$ términos es 
$$A(\mathcal{D})\approx \frac{A(\mathcal{R})}{N}\sum_{i=1}^{N}b_i.$$
Step12:  De nuevo, comparemos con el valor real
Step13:  Escribamos una función que tenga como entradas
Step14:  Actividad. Utilizar la anterior función para aproximar el área de la región descrita por
$$4(2x-1)^4+8(2y-1)^8<1+2(2y-1)^3(3x-2)^2$$
Poner los resultados en una tabla cuyas filas correspondan a la cantidad de términos utilizados en la aproximación (usar 10, 100, 1000, 10000 y 100000 términos).
<ASSISTANT_TASK:>
Python Code:
from IPython.display import YouTubeVideo
YouTubeVideo('Ti5zUD08w5s')
YouTubeVideo('jmsFC0mNayM')
Explanation: Integración numérica Montecarlo
Referencia:
- https://ocw.mit.edu/courses/mechanical-engineering/2-086-numerical-computation-for-mechanical-engineers-fall-2014/nutshells-guis/MIT2_086F14_Monte_Carlo.pdf
- http://ta.twi.tudelft.nl/mf/users/oosterle/oosterlee/lec8-hit-2009.pdf
- Sauer, Timothy. Análisis Numérico, 2da. Edición, ISBN: 978-607-32-2059-0.
<img style="float: center; margin: 0px 0px 15px 15px;" src="https://upload.wikimedia.org/wikipedia/commons/f/f2/Integral_as_region_under_curve.svg" width="300px" height="100px" />
Motivación
En análisis de ingeniería, normalmente debemos evaluar integrales definidas sobre un dominio complejo o en un espacio de dimensión alta.
Por ejemplo, podríamos querer calcular:
- la deflexión en una viga de geometría complicada, 
- el volumen de una parte tridimensional de una aeronave,
- o evaluar alguna medida de rendimiento (rentabilidad) en algún proceso que sea expresada como una integral de alguna función sin antiderivada primitiva (que se pueda expresar en términos de funciones elementales).
A la mano tenemos herramientas de integración analítica cuando tanto el espacio de integración como la función a integrar son simples. Cuando la función a integrar es difícil (incluso, imposible) de integrar podemos aún recurrir a métodos numéricos de integración. 
Desafortunadamente, los métodos determinísiticos de integración fallan cuando:
- la región es demasiado compleja para discretizarla,
- o la función a integrar es demasiado irregular,
- o la convergencia es demasiado lenta debido a la alta dimensionalidad del espacio de integración (ver Maldición de la dimensionalidad).
Por eso en esta clase veremos una técnica alternativa de integración numérica: Integración Montecarlo.
Ejemplos de funciones sin antiderivada primitiva.
De su curso de cálculo integral seguro recordarán (o estarán viendo) que existen funciones cuya integral no tiene primitiva. Es decir, que no podemos encontrar una función que se pueda expresar en forma de funciones elementales cuya derivada sea tal función.
Esto no significa que dicha función no se pueda integrar, ya que sabemos que cualquier función continua es integrable (y la mayoría de funciones que vemos a ese nivel, lo son). Lo que ocurre es que no podemos expresar dicha integral de una forma sencilla (por ejemplo, en función de exponenciales, senos, cosenos, logaritmos...).
Algunas integrales que no son elementales son:
- $\int e^{p(x)}\text{d}x$, donde $p(x)$ es un polinomio de grado mayor o igual a dos.
- $\int \frac{1}{log(x)}\text{d}x$.
- $\int \frac{sin(x)}{x}\text{d}x$
Referencia:
- https://www.gaussianos.com/funciones-sin-primitiva-elemental/
Ejemplos de regiones difíciles de discretizar.
End of explanation
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
def parab(x):
    return x**2
x = np.linspace(0,1)
y = parab(x)
plt.fill_between(x,y)
plt.text(0.8,0.2,'$\mathcal{D}$',fontsize=20)
plt.show()
Explanation: Integración Montecarlo tipo 1
Se basa en la definición de valor promedio de una función y en el valor esperado de una variable aleatoria uniforme.
Presentamos esto mediante un ejemplo.
Ejemplo. Aproxime el área bajo la curva $y=x^2$ en el intervalo $\left[0,1\right]$.
Veamos primero cómo luce dicha área.
End of explanation
help(np.random.uniform)
N = 100000
x = np.random.uniform(0, 1, N)
A_Dapprox = np.sum(parab(x))/N
A_Dapprox
Explanation: Entonces, lo que queremos es aproximar el área de la región $\mathcal{D}$. Llamaremos esta área $A(\mathcal{D})$.
De cálculo integral, sabemos que
$$A(\mathcal{D})=\int_{0}^{1}y\text{d}x=\int_{0}^{1}x^2\text{d}x$$.
Por definición, el valor promedio de una función $f:\left[a,b\right]\to\mathbb{R}$ en un intervalo $\left[a,b\right]$ es
$$\frac{1}{b-a}\int_{a}^{b}f(x)\text{d}x.$$
Entonces, el área bajo la curva $y=x^2$ es exactamente el valor promedio de $f(x)=x^2$ en $\left[0,1\right]$. Este valor promedio puede aproximarse mediante el promedio de los valores de la función en puntos aleatorios uniformemente distribuidos en el intervalo $\left[0,1\right]$. Es decir,
$$A(\mathcal{D})=\int_{0}^{1}x^2\text{d}x=\int_{0}^{1}f(x)\text{d}x\approx \frac{1}{N}\sum_{i=1}^{N}f(u_i)=\frac{1}{N}\sum_{i=1}^{N}u_i^2$$,
donde $u_i$ son realizaciones de la variable aleatoria $U\sim\mathcal{U}\left[0,1\right]$ ($U$ distribuye uniformemente en el intervalo $\left[0,1\right]$).
¿Cómo construit vectores de números aleatorios?
- Ver numpy.random.
En este caso necesitamos $N$ números aleatorios uniformemente distribuidos...
End of explanation
import pandas as pd
A_D = 1/3
N = np.logspace(1,7,7)
df = pd.DataFrame(index=N,columns=['Valor_aproximacion', 'Error_relativo'], dtype='float')
df.index.name = "Cantidad_terminos"
for n in N:
    x = np.random.uniform(0, 1, n.astype(int))
    df.loc[n,"Valor_aproximacion"] = np.sum(parab(x))/n
    df.loc[n,"Error_relativo"] = np.abs(df.loc[n,"Valor_aproximacion"]-A_D)/A_D
    
df
Explanation: En este caso, la integral se puede hacer fácilmente. Comparemos el resultado con el valor real:
$$A(\mathcal{D})=\int_{0}^{1}x^2\text{d}x=\left.\frac{x^3}{3}\right|_{x=0}^{x=1}=\frac{1}{3}$$
Hagamos una tabla viendo:
- cantidad de terminos
- valor de la aproximacion
- error relativo
End of explanation
# Escribir la función acá
def int_montecarlo1(f, a, b, N):
    return (b-a)/N*np.sum(f(np.random.uniform(a,b,N)))
Explanation: Ver que los resultados son distintos cada vez (¿porqué?). Sin embargo, se aproximan más o menos en la misma medida.
Aproximación de integrales en intervalos distintos a $\left[0,1\right]$.
Sin embargo, no todas las integrales que hacemos son en el intervalo $\left[0,1\right]$. En general, podemos integrar cualquier función continua en el intervalo $\left[a,b\right]$, donde $a,b\in\mathbb{R}$ con $a<b$.
Sea $f:\left[a,b\right]\to\mathbb{R}$ una función continua en el intervalo $\left(a,b\right)$ (por lo tanto es integrable endicho intervalo). Queremos resolver:
$$\int_{a}^{b}f(x)\text{d}x.$$
¿Cómo podemos usar la idea del valor promedio para resolver esto?
El valor promedio de $f$ en $\left[a,b\right]$ es:
$$\frac{1}{b-a}\int_{a}^{b}f(x)\text{d}x.$$
Este valor promedio puede aproximarse mediante el promedio de $N$ valores de la función en puntos aleatorios uniformemente distribuidos en el intervalo $\left[a,b\right]$. Es decir,
$$\frac{1}{b-a}\int_{a}^{b}f(x)\text{d}x\approx \frac{1}{N}\sum_{i=1}^{N}f(u_i)$$,
donde $u_i$ son realizaciones de la variable aleatoria $U\sim\mathcal{U}\left[a,b\right]$ ($U$ distribuye uniformemente en el intervalo $\left[a,b\right]$).
Finalmente, la aproximación montecarlo tipo 1 con $N$ términos es 
$$\int_{a}^{b}f(x)\text{d}x\approx \frac{b-a}{N}\sum_{i=1}^{N}f(u_i)$$,
Escribamos una función que tenga como entradas:
- la función a integrar $f$,
- los límites de integración $a$ y $b$, y
- el número de términos que se usará en la aproximación $N$,
y que devuelva la aproximación montecarlo tipo 1 de la integral $\int_{a}^{b}f(x)\text{d}x$.
End of explanation
# Resolver
def func1(x):
    return np.exp(x**2)
def func2(x):
    return 1/np.log(x)
def func3(x):
    return np.sin(x)/x
a, b = 4, 5
N = np.logspace(1,5,5)
df = pd.DataFrame(index=N,columns=['Funcion1', 'Funcion2', 'Funcion3'], dtype='float')
df.index.name = "Cantidad_terminos"
for n in N:
    df.loc[n,"Funcion1"] = int_montecarlo1(func1, a, b, n.astype(int))
    df.loc[n,"Funcion2"] = int_montecarlo1(func2, a, b, n.astype(int))
    df.loc[n,"Funcion3"] = int_montecarlo1(func3, a, b, n.astype(int))
    
df
Explanation: Actividad. Utilizar la anterior función para realizar las siguientes integrales. Poner los resultados en una tabla cuyas filas correspondan a la cantidad de términos utilizados en la aproximación (usar 10, 100, 1000, 10000 y 100000 términos) y cuyas columnas correspondan a las funciones.
- $\int_{4}^{5} e^{x^2}\text{d}x$.
- $\int_{4}^{5} \frac{1}{log(x)}\text{d}x$.
- $\int_{4}^{5} \frac{sin(x)}{x}\text{d}x$.
End of explanation
x = np.linspace(-0.1,1.1)
y = parab(x)
plt.plot(x,x,'k--',label='$y=x$')
plt.plot(x,y,'k',label='$y=x^2$')
plt.fill_between(x,x,y)
plt.text(0.5,0.4,'$\mathcal{D}$',fontsize=20)
plt.legend(loc='best')
plt.show()
Explanation: Integración Montecarlo tipo 2
Con la integración montecarlo tipo 1 pudimos aproximar integrales de funciones continuas de una variable en un intervalo dado. En realidad este mismo análisis se puede ampliar para aproximar integrales definidas de funciones continuas de varias variables (integrales sobre áreas, volúmenes e hipervolúmenes) dado que la noción de valor promedio de una función se extiende a cualquier dimensión.
Este es en realidad el caso interesante, pues las integrales de funciones complicadas también se pueden aproximar por métodos numéricos clásicos, pero cuando la dimensión aumenta es cuando montecarlo se vuelve una herramienta relevante. Dado que no lo veremos en clase por la limitación de que la mayoría no han visto cálculo en varias variables, este tema puede ser elegido como proyecto de módulo, donde se exploraría también como mejorar la aproximación de integrales montecarlo.
Como vimos en el ejemplo (y como debe ser claro de su curso de cálculo integral) una de las aplicaciones más importantes de la integración es hallar áreas. Y no solo el área bajo una curva, sino áreas entre curvas y áreas de regiones más complicadas.
Antes de ver la integración montecarlo tipo 2, ¿cómo podemos usar la integración montecarlo tipo 1 para aproximar el área entre curvas?
Ejemplo. Aproxime el área entre las curvas $y=x$, y $y=x^2$ en el intervalo $\left[0,1\right]$.
Veamos primero cómo luce dicha área.
End of explanation
# Usar la funcion int_montecarlo1
def f(x):
    return x-x**2
A_Daprox = int_montecarlo1(f, 0, 1, 100000000)
A_Daprox
Explanation: De cálculo integral, sabemos que
$$A(\mathcal{D})=\int_{0}^{1}x-x^2\text{d}x.$$
Entonces...
End of explanation
YouTubeVideo('G8fOTMYDPEA')
Explanation: De modo que si la región se puede describir fácilmente, diría el ferras 'no hay pedo, lo pago' (podemos usar montecarlo tipo 1).
End of explanation
def circ_arriba(x, r):
    return np.sqrt(r**2-x**2)
def circ_abajo(x, r):
    return -np.sqrt(r**2-x**2)
x = np.linspace(-1,1,100)
y1 = circ_arriba(x, 1)
y2 = circ_abajo(x, 1)
plt.figure(figsize=(5,5))
plt.plot(x,y1,'k')
plt.plot(x,y2,'k')
plt.fill_between(x,y1,y2)
plt.text(0,0,'$\mathcal{D}$',fontsize=20)
plt.text(0.8,0.8,'$\mathcal{R}$',fontsize=20)
plt.show()
Explanation: Pero, ¿qué pasa si la geometría de la región no se puede describir fácilmente?
Como en el caso anterior, motivaremos el método con un caso conocido. Vamos a aproximar el valor de $\pi$ usando el área de un círculo unitario.
Dibujemos el círculo unitario en la región $\mathcal{R}=\left[-1,1\right]\times\left[-1,1\right]$.
End of explanation
N = 1000000
x = np.random.uniform(-1, 1, N)
y = np.random.uniform(-1, 1, N)
X, Y = np.meshgrid(x,y)
plt.figure(figsize=(5,5))
plt.scatter(X,Y)
plt.show()
Explanation: Si aproximamos $A(\mathcal{D})$ aproximamos el valor de $\pi$, pues el área del círculo unitario es:
$$A(\mathcal{D})=\pi(1)^2=\pi.$$
Por otra parte es claro que el área de la región $\mathcal{R}=\left[-1,1\right]\times\left[-1,1\right]$ es
$$A(\mathcal{R})=4.$$
Ahora, haremos uso de nuestro generador de números aleatorios. Supongamos que escogemos un punto aleatorio en la región $\mathcal{R}=\left[-1,1\right]\times\left[-1,1\right]$. Describimos este punto como $(X,Y)$ para $X$ e $Y$ variables aleatorias uniformes sobre el intervalo $\left[-1,1\right]$.
¿Cómo generamos puntos aleatorios en un rectángulo?
End of explanation
def reg_circ(x,y):
    return x**2+y**2<=1
A_R = 4
A_Dapprox = A_R*np.sum(reg_circ(x,y))/N
A_Dapprox
Explanation: La probabilidad de que el punto $(X,Y)$ esté en el círculo unitario $\mathcal{D}$ es
$$P((X,Y)\in\mathcal{D})=\frac{A(\mathcal{D})}{A(\mathcal{R})}=\frac{\pi}{4}.$$
Luego, definimos una variable aleatoria de Bernoulli $B$ de manera que
$$B=\left\lbrace\begin{array}{ccc}0 & \text{si} & (X,Y)\notin\mathcal{D}\1 & \text{si} & (X,Y)\in\mathcal{D} \end{array}\right.=\left\lbrace\begin{array}{ccc}0 & \text{si} & X^2+Y^2>1\1 & \text{si} & X^2+Y^2\leq 1 \end{array}\right..$$
Entonces, el valor esperado de la variable aleatoria $B$ es
$$E\left[B\right]=\theta=P((X,Y)\in\mathcal{D})=\frac{A(\mathcal{D})}{A(\mathcal{R})}.$$
De lo anterior, una estimación de theta se puede obtener como
$$\theta=\frac{A(\mathcal{D})}{A(\mathcal{R})}\approx \frac{1}{N}\sum_{i=1}^{N}b_i,$$
donde 
$$b_i=\left\lbrace\begin{array}{ccc}0 & \text{si} & x_i^2+y_i^2>1\1 & \text{si} & x_i^2+y_i^2\leq 1 \end{array}\right.$$
son realizaciones de la variable aleatoria $B$, que a su vez es producto de las realizaciones $x_i$ e $y_i$ de las variables aleatorias $X$ e $Y$, respectivamente.
Finalmente, la aproximación montecarlo tipo 2 con $N$ términos es 
$$A(\mathcal{D})\approx \frac{A(\mathcal{R})}{N}\sum_{i=1}^{N}b_i.$$
End of explanation
A_D = np.pi
N = np.logspace(1,7,7)
df = pd.DataFrame(index=N,columns=['Valor_aproximacion', 'Error_relativo'], dtype='float')
df.index.name = "Cantidad_terminos"
for n in N:
    x = np.random.uniform(-1, 1, n.astype(int))
    y = np.random.uniform(-1, 1, n.astype(int))
    df.loc[n,"Valor_aproximacion"] = A_R*np.sum(reg_circ(x,y))/n
    df.loc[n,"Error_relativo"] = np.abs(df.loc[n,"Valor_aproximacion"]-A_D)/A_D
    
df
Explanation: De nuevo, comparemos con el valor real
End of explanation
# Escribir la función acá
def int_montecarlo2(region, a1, b1, a2, b2, N):
    A_R = (b1-a1)*(b2-a2)
    x = np.random.uniform(a1, b1, N.astype(int))
    y = np.random.uniform(a2, b2, N.astype(int))
    return A_R*np.sum(region(x,y))/N
Explanation: Escribamos una función que tenga como entradas:
- la función que describe la region $region$,
- los límites de la region $a_1$, $b_1$, $a_2$ y $b_2$, con $R=\left[a_1,b_1\right]\times\left[a_2,b_2\right]$ y
- el número de términos que se usará en la aproximación $N$,
y que devuelva la aproximación montecarlo tipo 2 del area de la region.
End of explanation
N = 100
x = np.linspace(0, 1, N)
y = np.linspace(0, 1, N)
def region(x,y):
    return 4*(2*x-1)**4+8*(2*y-1)**8 < 1+2*(2*y-1)**3*(3*x-2)**2
X, Y = np.meshgrid(x,y)
plt.figure(figsize=(5,5))
plt.scatter(X,Y,c=~region(X,Y),cmap='bone')
plt.show()
# Resolver
a1, a2, b1, b2 = 0, 0, 1, 1
N = np.logspace(1,5,5)
df = pd.DataFrame(index=N,columns=['Valor_aproximacion'], dtype='float')
df.index.name = "Cantidad_terminos"
for n in N:
    df.loc[n,"Valor_aproximacion"] = int_montecarlo2(region, a1, b1, a2, b2, n)
        
df
Explanation: Actividad. Utilizar la anterior función para aproximar el área de la región descrita por
$$4(2x-1)^4+8(2y-1)^8<1+2(2y-1)^3(3x-2)^2$$
Poner los resultados en una tabla cuyas filas correspondan a la cantidad de términos utilizados en la aproximación (usar 10, 100, 1000, 10000 y 100000 términos).
End of explanation
<END_TASK>