Integración de Monte Carlo

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 por 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, métodos como la regla trapezoidal emplean un enfoque determinista . La integración de Monte Carlo, en cambio, utiliza 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 sus respectivos márgenes de error, y es probable que el valor correcto se encuentre dentro de dichos márgenes.

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

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

Puedo ser aproximado por $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 $\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 pueden estimarse mediante la varianza de la muestra utilizando la estimación insesgada de la varianza .

$\mathrm {Var} (f)=\mathrm {E} (\sigma _{N}^{2})\equiv {\frac {1}{N-1}}\sum _{i=1}^{N}\mathrm {E} \left[\left(f({\overline {\mathbf {x} }}_{i})-\langle f\rangle \right)^{2}\right].$ lo que conduce a $\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 {\mathrm {E} (\sigma {N}^{2})}{N}}.$

Dado que la secuencia $\left\{\mathrm {E} (\sigma _{1}^{2}),\mathrm {E} (\sigma _{2}^{2}),\mathrm {E} (\sigma _{3}^{2}),\ldots \right\}$ está acotada debido a que es idénticamente igual a Var(f) , siempre que se suponga que es finito, esta varianza disminuye asintóticamente a cero como 1/ N. La estimación del error de Q _N es, por lo tanto, $\delta Q_{N}\approx {\sqrt {\mathrm {Var} (Q_{N})}}=V{\frac {\sqrt {\mathrm {Var} (f)}}{\sqrt {N}}},$ que disminuye a medida que ${\tfrac {1}{\sqrt {N}}}$ . Este es el error estándar de la media multiplicado por $V$ Este resultado no depende del número de dimensiones de la integral, lo cual 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 una cota de error estricta; el muestreo aleatorio puede no revelar todas las características importantes del integrando, lo que puede resultar en una subestimación del error.

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

Ejemplo

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

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

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

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

Ejemplo en C/C++

#include <stdio.h> #include <stdlib.h> #include <time.h> int main () { // Inicializa el número de conteos a 0, estableciendo el número total en 100000 en el bucle. int throws = 99999 , insideCircle = 0 ; double randX , randY , pi ;srand ( tiempo ( NULL ));// Comprueba para cada par aleatorio de x e y si están dentro de un círculo de radio 1. for ( int i = 0 ; i < throws ; i ++ ) { randX = rand () / ( double ) RAND_MAX ; randY = rand () / ( double ) RAND_MAX ; if ( randX * randX + randY * randY < 1 ) { insideCircle ++ ; } }// Calculando pi e imprimiendo. pi = 4.0 * insideCircle / throws ; printf ( "%lf \n " , pi ); }

Ejemplo de Python

Hecho en Python .

import numpy as nprng = np.random.default_rng ( 0 )lanzamientos = 2000 radio = 1# Elige datos X e Y aleatorios centrados en 0,0 x = rng . uniform ( - radio , radio , lanzamientos ) y = rng . uniform ( - radio , radio , lanzamientos )# Cuenta las veces que (x, y) está dentro del círculo, # lo cual sucede cuando sqrt(x^2 + y^2) <= radio. inside_circle = np . count_nonzero ( np . hypot ( x , y ) <= radio )# Calcular el área e imprimirla; 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 que aparece a continuación describe un proceso de integración de la función. $f(x)={\frac {1}{1+\sinh(2x)\log(x)^{2}}}$ de $0.8<x<3$ Utilizando el método de Montecarlo en Mathematica :

func [ x_ ] := 1 / ( 1 + Sinh [ 2 * x ] * ( Iniciar sesión [ x ]) ^ 2 );(*Muestra de la distribución normal truncada para acelerar la convergencia*) Distrib [ x_ , average_ , var_ ] := PDF [ NormalDistribution [ average , var ], 1.1 * x - 0.1 ]; n = 10 ; RV = RandomVariate [ TruncatedDistribution [{ 0.8 , 3 }, NormalDistribution [ 1 , 0.399 ]], n ]; Int = 1 / n Total [ func [ RV ] / Distrib [ RV , 1 , 0.399 ]] * Integrate [ Distrib [ x , 1 , 0.399 ], { x , 0.8 , 3 }]NIntegrate [ func [ x ], { x , 0.8 , 3 }] (*Comparar con la respuesta real*)

Muestreo estratificado recursivo

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

El muestreo estratificado recursivo es una generalización de las cuadraturas adaptativas unidimensionales a integrales multidimensionales. En cada paso de la recursión, la integral y el error se estiman mediante un algoritmo de 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 estos subvolúmenes.

La estrategia habitual de «dividir por dos» no funciona para múltiples dimensiones, ya que el número de subvolúmenes crece demasiado rápido como para poder controlarlo. En su lugar, se estima a lo largo de qué dimensión una subdivisión debería generar mayores beneficios y solo se subdivide el volumen a lo largo de 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.

MISER Monte Carlo

El algoritmo MISER se basa en el muestreo estratificado recursivo . Esta técnica tiene como objetivo reducir el error de integración global concentrando 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 la integral $E_{a}(f)$ y $E_{b}(f)$ y variaciones $\sigma _{a}^{2}(f)$ y $\sigma _{b}^{2}(f)$ , la varianza Var( f ) de la estimación combinada $E(f)={\tfrac {1}{2}}\left(E_{a}(f)+E_{b}(f)\right)$ está dado por, $\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, ${\frac {N_{a}}{N_{a}+N_{b}}}={\frac {\sigma _{a}}{\sigma _{a}+\sigma _{b}}}$

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

El algoritmo MISER procede dividiendo la región de integración por la mitad a lo largo de un eje de coordenadas para obtener dos subregiones en cada paso. La dirección se elige examinando todas las posibles bisecciones y seleccionando la que minimice la varianza combinada de las dos subregiones. La varianza en las subregiones se estima muestreando con una fracción del número total de puntos disponibles en el paso actual. El mismo procedimiento se repite recursivamente para cada uno de los dos semiplanos 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. Estos valores individuales y sus estimaciones de error se combinan para obtener un resultado global y una estimación de su error.

Muestreo de importancia

Existen diversos algoritmos de muestreo por importancia, como por ejemplo:

Algoritmo de muestreo por importancia

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

Consideremos el siguiente ejemplo: 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 pequeña parte de ellas será significativa para la integral. Esto se puede mejorar eligiendo una distribución diferente a la de 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 de la función a integrar.

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

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

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

VEGAS Monte Carlo

El algoritmo VEGAS aproxima la distribución exacta realizando varias pasadas sobre la región de integración, lo 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 el número de intervalos del histograma crezca como K ^d , la distribución de probabilidad se aproxima mediante una función separable: $g(x_{1},x_{2},\ldots )=g_{1}(x_{1})g_{2}(x_{2})\ldots$ de modo que el número de intervalos requeridos sea solo Kd . Esto equivale 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 aproximadamente separable, esto aumentará la eficiencia de la integración con VEGAS. VEGAS incorpora varias características adicionales y combina el muestreo estratificado y el muestreo por importancia. ^{[ 9 ]}

Véase también

Notas

↑ Press et al. 2007 , Cap. 4
↑ Press et al. 2007 , Cap. 7
1 2 3 4 Newman y Barkema 1999 , Cap. 2
↑ Newman y Barkema 1999 , Cap. 1
↑ Prensa y otros, 2007
↑ MacKay 2003 , págs. 284–292
↑ Press & Farrar 1990 , págs. 190–195
^ Kroese, Taimre y Botev 2011
1 2 Lepage 1978

Referencias

Caflisch, RE (1998). "Métodos de Montecarlo y cuasi-Montecarlo". Acta Numérica . 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 de Monte Carlo multidimensional" . 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.{{cite book}}: Incompatibilidad de ISBN/Fecha ( ayuda )
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). Numerical Recipes: The Art of Scientific Computing (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. MR 2012999 .
Newman, MEJ; Barkema, GT (1999). Métodos de Monte Carlo en física estadística . Clarendon Press.
Roberto, CP; Casella, G (2004). Métodos estadísticos de Monte Carlo (2ª ed.). Saltador. ISBN 978-1-4419-1939-7.

Enlaces externos

Matemáticas de café : 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 ingenua de Monte Carlo: Documentación para las rutinas ingenuas de Monte Carlo en C++
Aplicación de la metodología de Monte Carlo en problemas de física estadística.

[1] Press et al. 2007 , Cap. 4

[2] Press et al. 2007 , Cap. 7

[newman1999ch2-3] 1 2 3 4 Newman y Barkema 1999 , Cap. 2

[newman1999ch1-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

[kr11-8] Kroese, Taimre y Botev 2011

[Lepage,_1978-9] 1 2 Lepage 1978

[ 1 ]

[ 2 ]

[ 3 ]

[ 4 ]

[ 5 ]

[ 6 ]

[ 7 ]

[ 8 ]

[ 9 ]