Integración de Monte Carlo

Técnica numérica

Una ilustración de la integración de Monte Carlo. En este ejemplo, el dominio D es el círculo interior y el dominio E es el cuadrado. Como el área del cuadrado (4) se puede calcular fácilmente, el área del círculo (π*1,0 ² ) se puede estimar mediante la relación (0,8) de los puntos dentro del círculo (40) con el número total de puntos (50), lo que da una aproximación para el área del círculo de 4*0,8 = 3,2 ≈ π.

En matemáticas , la integración de Monte Carlo es una técnica de integración numérica que utiliza números aleatorios . Es un método particular de Monte Carlo que calcula numéricamente una integral definida . Mientras que otros algoritmos suelen evaluar el integrando en una cuadrícula regular, ^[1] Monte Carlo elige aleatoriamente puntos en los que se evalúa el integrando. ^[2] Este método es particularmente útil para integrales de dimensiones superiores. ^[3]

Existen diferentes métodos para realizar una integración de Monte Carlo, como el muestreo uniforme , el muestreo estratificado , el muestreo de importancia , el Monte Carlo secuencial (también conocido como filtro de partículas) y los métodos de partículas de campo medio .

Descripción general

En la integración numérica, los métodos como la regla trapezoidal utilizan un enfoque determinista . La integración de Monte Carlo, por otro lado, emplea un enfoque no determinista : cada realización proporciona un resultado diferente. En Monte Carlo, el resultado final es una aproximación del valor correcto con las barras de error respectivas, y es probable que el valor correcto se encuentre dentro de esas barras de error.

El problema que aborda la integración de Monte Carlo es el cálculo de una integral definida multidimensional donde Ω, un subconjunto de , tiene volumen $${\displaystyle I=\int _{\Omega }f({\overline {\mathbf {x} }})\,d{\overline {\mathbf {x} }}}$$ ${\displaystyle \mathbb {R} ^{m}}$ $${\displaystyle V=\int _{\Omega }d{\overline {\mathbf {x} }}}$$

El enfoque ingenuo de Monte Carlo consiste en muestrear puntos uniformemente en Ω: ^[4] dadas N muestras uniformes, $${\displaystyle {\overline {\mathbf {x} }}_{1},\cdots ,{\overline {\mathbf {x} }}_{N}\in \Omega ,}$$

Me puedo aproximar por $${\displaystyle I\approx Q_{N}\equiv V{\frac {1}{N}}\sum _{i=1}^{N}f({\overline {\mathbf {x} }}_{i})=V\langle f\rangle .}$$

Esto se debe a que la ley de los grandes números garantiza que $${\displaystyle \lim _{N\to \infty }Q_{N}=I.}$$

Dada la estimación de I a partir de Q _N , las barras de error de Q _N se pueden estimar mediante la varianza de la muestra utilizando la estimación imparcial de la varianza .

$${\displaystyle \mathrm {Var} (f)\equiv \sigma _{N}^{2}={\frac {1}{N-1}}\sum _{i=1}^{N}\left(f({\overline {\mathbf {x} }}_{i})-\langle f\rangle \right)^{2}.}$$Lo que conduce a $${\displaystyle \mathrm {Var} (Q_{N})={\frac {V^{2}}{N^{2}}}\sum _{i=1}^{N}\mathrm {Var} (f)=V^{2}{\frac {\mathrm {Var} (f)}{N}}=V^{2}{\frac {\sigma _{N}^{2}}{N}}.}$$

Mientras la secuencia esté acotada, esta varianza disminuye asintóticamente a cero como 1/ N . La estimación del error de Q _N es entonces que disminuye como . Este es el error estándar de la media multiplicado por . Este resultado no depende del número de dimensiones de la integral, que es la ventaja prometida de la integración de Monte Carlo frente a la mayoría de los métodos deterministas que dependen exponencialmente de la dimensión. ^[5] Es importante notar que, a diferencia de los métodos deterministas, la estimación del error no es un límite de error estricto; el muestreo aleatorio puede no descubrir todas las características importantes del integrando que pueden resultar en una subestimación del error. $${\displaystyle \left\{\sigma _{1}^{2},\sigma _{2}^{2},\sigma _{3}^{2},\ldots \right\}}$$ $${\displaystyle \delta Q_{N}\approx {\sqrt {\mathrm {Var} (Q_{N})}}=V{\frac {\sigma _{N}}{\sqrt {N}}},}$$ ${\displaystyle {\tfrac {1}{\sqrt {N}}}}$ ${\displaystyle V}$

Si bien el método Monte Carlo ingenuo funciona para ejemplos simples, solo se puede lograr una mejora con respecto a los algoritmos deterministas con algoritmos que utilicen distribuciones de muestreo específicas del problema. Con una distribución de muestreo adecuada, es posible aprovechar el hecho de que casi todos los integrandos de dimensiones superiores están muy localizados y solo un subespacio pequeño contribuye notablemente a la integral. ^[6] Una gran parte de la literatura de Monte Carlo se dedica al desarrollo de estrategias para mejorar las estimaciones de error. En particular, el muestreo estratificado (dividir la región en subdominios) y el muestreo de importancia (muestreo de distribuciones no uniformes) son dos ejemplos de tales técnicas.

Ejemplo

Error relativo en función del número de muestras, mostrando la escala ${\displaystyle {\tfrac {1}{\sqrt {N}}}}$

Un ejemplo paradigmático de una integración de Monte Carlo es la estimación de π. Considere la función y el conjunto Ω = [−1,1] × [−1,1] con V = 4. Observe que $${\displaystyle H\left(x,y\right)={\begin{cases}1&{\text{if }}x^{2}+y^{2}\leq 1\\0&{\text{else}}\end{cases}}}$$ $${\displaystyle I_{\pi }=\int _{\Omega }H(x,y)dxdy=\pi .}$$

Por lo tanto, una forma básica de calcular el valor de π con la integración de Monte Carlo es elegir N números aleatorios en Ω y calcular $${\displaystyle Q_{N}=4{\frac {1}{N}}\sum _{i=1}^{N}H(x_{i},y_{i})}$$

En la figura de la derecha, el error relativo se mide en función de N , lo que confirma el . ${\displaystyle {\tfrac {Q_{N}-\pi }{\pi }}}$ ${\displaystyle {\tfrac {1}{\sqrt {N}}}}$

Ejemplo C

Tenga en cuenta que se debe utilizar un verdadero generador de números aleatorios.

int i , lanza = 99999 , circuloInterior = 0 ; doble randX , randY , pi ;          srand ( tiempo ( NULL ));para ( i = 0 ; i < lanza ; ++ i ) { randX = rand () / ( doble ) RAND_MAX ; randY = rand () / ( doble ) RAND_MAX ; si ( randX * randX + randY * randY < 1 ) ++ circuloInterior ; }                               pi = 4.0 * circulointerior / lanza ;      

Ejemplo de Python

Hecho en Python .

Desde  numpy,  importar  al azarlanza  =  2000 círculo_interior  =  0 i  =  0 radio  =  1 mientras  i  <  lanza : # Elige X  e Y aleatorios centrados alrededor de 0,0  x = aleatorio.uniforme  (  -radio , radio ) y = aleatorio.uniforme ( -radio , radio ) # Si el punto está dentro del círculo, aumenta la variable si x ** 2 + y ** 2 < = radio ** 2 : círculo_interior + = 1 i += 1                  # Calcular el área e imprimir; debería estar más cerca de Pi con un número creciente de lanzamientos área  =  ((( 2  *  radio )  **  2 )  *  círculo_interior )  /  lanzamientos imprimir ( área )

Ejemplo de Wolfram Mathematica

El código siguiente describe un proceso de integración de la función mediante el método de Monte Carlo en Mathematica : $${\displaystyle f(x)={\frac {1}{1+\sinh(2x)\log(x)^{2}}}}$$ ${\displaystyle 0.8<x<3}$

func [ x_ ] := 1 / ( 1 + Sinh [ 2 * x ] * ( Iniciar sesión [ x ]) ^ 2 );    (*Muestra de una distribución normal truncada para acelerar la convergencia*) Distrib [ x_ , promedio_ , var_ ] := PDF [ DistribuciónNormal [ promedio , var ], 1.1 * x - 0.1 ]; n = 10 ; RV = VariableAleatoria [ DistribuciónTruncada [{ 0.8 , 3 }, DistribuciónNormal [ 1 , 0.399 ]], n ]; Int = 1 / n Total [ func [ RV ] / Distrib [ RV , 1 , 0.399 ]] * Integrar [ Distrib [ x , 1 , 0.399 ], { x , 0.8 , 3 }]                          NIntegrate [ func [ x ], { x , 0.8 , 3 }] (*Comparar con la respuesta real*)    

Muestreo estratificado recursivo

Ilustración de muestreo estratificado recursivo. En este ejemplo, la función de la ilustración anterior se integró dentro de un cuadrado unitario utilizando el algoritmo sugerido. Los puntos muestreados se registraron y se graficaron. Es claro que el algoritmo de muestreo estratificado concentra los puntos en las regiones donde la variación de la función es mayor. ${\displaystyle f(x,y)={\begin{cases}1&x^{2}+y^{2}<1\\0&x^{2}+y^{2}\geq 1\end{cases}}}$

El muestreo estratificado recursivo es una generalización de las cuadraturas adaptativas unidimensionales a integrales multidimensionales. En cada paso de recursión, la integral y el error se estiman utilizando un algoritmo Monte Carlo simple. Si la estimación del error es mayor que la precisión requerida, el volumen de integración se divide en subvolúmenes y el procedimiento se aplica recursivamente a los subvolúmenes.

La estrategia habitual de "dividir por dos" no funciona para volúmenes multidimensionales, ya que el número de subvolúmenes crece demasiado rápido como para llevar un registro. En lugar de ello, se calcula en qué dimensión una subdivisión debería generar más dividendos y solo se subdivide el volumen en esa dimensión.

El algoritmo de muestreo estratificado concentra los puntos de muestreo en las regiones donde la varianza de la función es mayor, reduciendo así la varianza general y haciendo que el muestreo sea más efectivo, como se muestra en la ilustración.

La popular rutina MISER implementa un algoritmo similar.

AVIADOR Monte Carlo

El algoritmo MISER se basa en un muestreo estratificado recursivo . Esta técnica tiene como objetivo reducir el error de integración global al concentrar los puntos de integración en las regiones de mayor varianza. ^[7]

La idea del muestreo estratificado comienza con la observación de que para dos regiones disjuntas a y b con estimaciones de Monte Carlo de las varianzas integrales y y y , la varianza Var( f ) de la estimación combinada está dada por, ${\displaystyle E_{a}(f)}$ ${\displaystyle E_{b}(f)}$ ${\displaystyle \sigma _{a}^{2}(f)}$ ${\displaystyle \sigma _{b}^{2}(f)}$ $${\displaystyle E(f)={\tfrac {1}{2}}\left(E_{a}(f)+E_{b}(f)\right)}$$ $${\displaystyle \mathrm {Var} (f)={\frac {\sigma _{a}^{2}(f)}{4N_{a}}}+{\frac {\sigma _{b}^{2}(f)}{4N_{b}}}}$$

Se puede demostrar que esta varianza se minimiza distribuyendo los puntos de tal manera que, $${\displaystyle {\frac {N_{a}}{N_{a}+N_{b}}}={\frac {\sigma _{a}}{\sigma _{a}+\sigma _{b}}}}$$

Por lo tanto, la estimación del error más pequeño se obtiene asignando puntos de muestra en proporción a la desviación estándar de la función en cada subregión.

El algoritmo MISER procede dividiendo en dos la región de integración a lo largo de un eje de coordenadas para dar dos subregiones en cada paso. La dirección se elige examinando todas las d bisecciones posibles y seleccionando la que minimice la varianza combinada de las dos subregiones. La varianza en las subregiones se estima mediante un muestreo con una fracción del número total de puntos disponibles para el paso actual. Luego se repite el mismo procedimiento de forma recursiva para cada uno de los dos semiespacios a partir de la mejor bisección. Los puntos de muestra restantes se asignan a las subregiones utilizando la fórmula para N _a y N _b . Esta asignación recursiva de puntos de integración continúa hasta una profundidad especificada por el usuario donde cada subregión se integra utilizando una estimación de Monte Carlo simple. Luego, estos valores individuales y sus estimaciones de error se combinan hacia arriba para dar un resultado general y una estimación de su error.

Muestreo de importancia

Hay una variedad de algoritmos de muestreo de importancia, como

Algoritmo de muestreo por importancia

El muestreo por importancia proporciona una herramienta muy importante para realizar la integración de Montecarlo. ^{[3] [8]} El principal resultado del muestreo por importancia para este método es que el muestreo uniforme de es un caso particular de una elección más genérica, en la que las muestras se extraen de cualquier distribución . La idea es que se puede elegir para disminuir la varianza de la medición Q _N . ${\displaystyle {\overline {\mathbf {x} }}}$ ${\displaystyle p({\overline {\mathbf {x} }})}$ ${\displaystyle p({\overline {\mathbf {x} }})}$

Consideremos el siguiente ejemplo, en el que se desea integrar numéricamente una función gaussiana centrada en 0, con σ = 1, desde −1000 hasta 1000. Naturalmente, si las muestras se extraen uniformemente en el intervalo [−1000, 1000], solo una parte muy pequeña de ellas sería significativa para la integral. Esto se puede mejorar eligiendo una distribución diferente de la que se utiliza para las muestras, por ejemplo, muestreando según una distribución gaussiana centrada en 0, con σ = 1. Por supuesto, la elección "correcta" depende en gran medida del integrando.

Formalmente, dado un conjunto de muestras elegidas de una distribución, el estimador para I viene dado por ^[3] $${\displaystyle p({\overline {\mathbf {x} }}):\qquad {\overline {\mathbf {x} }}_{1},\cdots ,{\overline {\mathbf {x} }}_{N}\in V,}$$ $${\displaystyle Q_{N}\equiv {\frac {1}{N}}\sum _{i=1}^{N}{\frac {f({\overline {\mathbf {x} }}_{i})}{p({\overline {\mathbf {x} }}_{i})}}}$$

Intuitivamente, esto dice que si elegimos una muestra particular con el doble de peso que otras muestras, la ponderamos con la mitad de peso que las otras muestras. Este estimador es naturalmente válido para el muestreo uniforme, caso en el que es constante. ${\displaystyle p({\overline {\mathbf {x} }})}$

El algoritmo Metropolis-Hastings es uno de los algoritmos más utilizados para generar a partir de , ^[3] proporcionando así una forma eficiente de calcular integrales. ${\displaystyle {\overline {\mathbf {x} }}}$ ${\displaystyle p({\overline {\mathbf {x} }})}$

VEGAS Montecarlo

El algoritmo VEGAS aproxima la distribución exacta haciendo una serie de pasadas sobre la región de integración que crea el histograma de la función f . Cada histograma se utiliza para definir una distribución de muestreo para la siguiente pasada. Asintóticamente, este procedimiento converge a la distribución deseada. ^[9] Para evitar que la cantidad de compartimentos del histograma crezca como K ^d , la distribución de probabilidad se aproxima mediante una función separable: de modo que la cantidad de compartimentos necesarios sea solo Kd . Esto es equivalente a localizar los picos de la función a partir de las proyecciones del integrando sobre los ejes de coordenadas. La eficiencia de VEGAS depende de la validez de esta suposición. Es más eficiente cuando los picos del integrando están bien localizados. Si un integrando se puede reescribir en una forma que sea aproximadamente separable, esto aumentará la eficiencia de la integración con VEGAS. VEGAS incorpora una serie de características adicionales y combina tanto el muestreo estratificado como el muestreo de importancia. ^[9] $${\displaystyle g(x_{1},x_{2},\ldots )=g_{1}(x_{1})g_{2}(x_{2})\ldots }$$

Véase también

• Método cuasi-Monte Carlo
• Campo auxiliar Monte Carlo
• El método de Monte Carlo en física estadística
• Método de Monte Carlo
• Reducción de la varianza

Notas

1. Press et al. 2007, cap. 4
2. Press et al. 2007, cap. 7
3. Newman y Barkema 1999, cap. 2
4. Newman y Barkema 1999, cap. 1
5. Prensa y otros, 2007
6. MacKay 2003, págs. 284-292
7. Press & Farrar 1990, págs. 190-195
8. Kroese, Taimre y Botev 2011
9. Lepage 1978

Referencias

• Caflisch, RE (1998). "Métodos de Monte Carlo y cuasi-Monte Carlo". Acta Numerica . 7 : 1–49. Bibcode :1998AcNum...7....1C. doi :10.1017/S0962492900002804. S2CID  5708790.
• Weinzierl, S. (2000). "Introducción a los métodos de Monte Carlo". arXiv : hep-ph/0006269 .
• Press, WH; Farrar, GR (1990). "Muestreo estratificado recursivo para la integración multidimensional de Monte Carlo". Computers in Physics . 4 (2): 190. Bibcode :1990ComPh...4..190P. doi : 10.1063/1.4822899 .
• Lepage, GP (1978). "Un nuevo algoritmo para la integración multidimensional adaptativa". Journal of Computational Physics . 27 (2): 192–203. Bibcode :1978JCoPh..27..192L. doi :10.1016/0021-9991(78)90004-9.
• Lepage, GP (1980). "VEGAS: Un programa de integración multidimensional adaptativo". Preimpresión de Cornell CLNS 80-447 .
• Hammersley, JM; Handscomb, DC (1964). Métodos de Monte Carlo . Methuen. ISBN 978-0-416-52340-9.
• Kroese, DP ; Taimre, T.; Botev, ZI (2011). Manual de métodos de Monte Carlo . John Wiley & Sons.
• Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007). Recetas numéricas: el arte de la computación científica (3.ª ed.). Nueva York: Cambridge University Press. ISBN 978-0-521-88068-8.
• MacKay, David (2003). "capítulo 4.4 Tipicidad y capítulo 29.1" (PDF) . Teoría de la información, inferencia y algoritmos de aprendizaje. Cambridge University Press. ISBN 978-0-521-64298-9. Sr.  2012999.
• Newman, MEJ; Barkema, GT (1999). Métodos de Monte Carlo en Física Estadística . Clarendon Press.
• Robert, CP; Casella, G (2004). Métodos estadísticos de Monte Carlo (2.ª ed.). Springer. ISBN 978-1-4419-1939-7.

Enlaces externos

• Café math: Integración de Monte Carlo: Un artículo de blog que describe la integración de Monte Carlo (principio, hipótesis, intervalo de confianza)
• Boost.Math: Integración de Naive Monte Carlo: Documentación para las rutinas de Naive Monte Carlo de C++
• Aplicación del applet de Monte Carlo en problemas de física estadística