Integración de Monte Carlo - Monte Carlo integration

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. Debido a que el área del cuadrado (4) se puede calcular fácilmente, el área del círculo (π * 1.0 ² ) se puede estimar por la razón (0.8) de los puntos dentro del círculo (40) al número total de puntos (50) , lo que arroja 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 para la integración numérica utilizando números aleatorios . Es un método particular de Monte Carlo que calcula numéricamente una integral definida . Mientras que otros algoritmos generalmente evalúan el integrando en una cuadrícula regular, Monte Carlo elige al azar puntos en los que se evalúa el integrando. Este método es particularmente útil para integrales de dimensiones superiores.

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

Visión general

En la integración numérica, 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 respectivas barras de error, y es probable que el valor correcto esté 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

${\ Displaystyle I = \ int _ {\ Omega} f ({\ overline {\ mathbf {x}}}) \, d {\ overline {\ mathbf {x}}}}$

donde Ω, un subconjunto de R ^m , tiene volumen

${\ Displaystyle V = \ int _ {\ Omega} d {\ overline {\ mathbf {x}}}}$

El enfoque ingenuo de Monte Carlo es muestrear puntos uniformemente en Ω: dadas N muestras uniformes,

${\ Displaystyle {\ overline {\ mathbf {x}}} _ {1}, \ cdots, {\ overline {\ mathbf {x}}} _ {N} \ in \ Omega,}$

Puedo ser aproximado 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 asegura 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 pueden estimarse mediante la varianza de la muestra utilizando la estimación insesgada 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 lleva 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

${\ Displaystyle \ left \ {\ sigma _ {1} ^ {2}, \ sigma _ {2} ^ {2}, \ sigma _ {3} ^ {2}, \ ldots \ right \}}$

está limitada, esta variación disminuye asintóticamente a cero como 1 / N . Por tanto, la estimación del error de Q _N es

${\ Displaystyle \ delta Q_ {N} \ approx {\ sqrt {\ mathrm {Var} (Q_ {N})}} = V {\ frac {\ sigma _ {N}} {\ sqrt {N}}}, }$

que disminuye a medida que . Este es el error estándar de la media multiplicada 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. 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 {\ tfrac {1} {\ sqrt {N}}}}$${\ Displaystyle V}$

Si bien el ingenuo Monte Carlo funciona con ejemplos simples, una mejora con respecto a los algoritmos deterministas solo se puede lograr con algoritmos que utilizan distribuciones de muestreo específicas del problema. Con una distribución de muestra adecuada, 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. Una gran parte de la literatura de Monte Carlo está dedicada 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

${\ Displaystyle H \ left (x, y \ right) = {\ begin {cases} 1 & {\ text {if}} x ^ {2} + y ^ {2} \ leq 1 \\ 0 & {\ text {else }} \ end {cases}}}$

y el conjunto Ω = [−1,1] × [−1,1] con V = 4. Observa que

${\ Displaystyle I _ {\ pi} = \ int _ {\ Omega} H (x, y) dxdy = \ pi.}$

Por lo tanto, una forma burda 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}}}}$

C ejemplo

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

int i, throws = 99999, insideCircle = 0;
double randX, randY, pi;

srand(time(NULL));

for (i = 0; i < throws; ++i) {
  randX = rand() / (double) RAND_MAX;
  randY = rand() / (double) RAND_MAX;
  if (randX * randX + randY * randY < 1) ++insideCircle;
}

pi = 4.0 * insideCircle / throws;

Ejemplo de Wolfram Mathematica

El siguiente código describe un proceso de integración de la función

${\ Displaystyle f (x) = {\ frac {1} {1+ \ sinh (2x) \ log (x) ^ {2}}}}$

de usar el método Monte-Carlo en Mathematica : ${\ Displaystyle 0.8 <x <3}$

func[x_] := 1/(1 + Sinh[2*x]*(Log[x])^2);

(*Sample from truncated normal distribution to speed up convergence*)
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}] (*Compare with real answer*)

Muestreo estratificado recursivo

Una 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 fueron registrados y graficados. El algoritmo de muestreo claramente 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 cuadraturas adaptativas unidimensionales a integrales multidimensionales. En cada paso de recursividad, la integral y el error se estiman utilizando 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 de forma recursiva a los subvolúmenes.

La estrategia ordinaria de 'dividir entre dos' no funciona para múltiples dimensiones, ya que el número de subvolúmenes crece demasiado rápido para realizar un seguimiento. En su lugar, se estima en qué dimensión una subdivisión debería generar la mayor cantidad de dividendos y solo subdivide el volumen a lo largo de esta 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 gran varianza 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 un muestreo estratificado recursivo . Esta técnica tiene como objetivo reducir el error de integración general al concentrar los puntos de integración en las regiones de mayor varianza.

La idea de muestreo estratificado comienza con la observación de que para dos disjuntos regiones una y b con Monte Carlo estimaciones de la integral y y varianzas y , la varianza Var ( f ) de la estimación combinado ${\ 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)}$

es dado por,

${\ 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 manera que,

${\ Displaystyle {\ frac {N_ {a}} {N_ {a} + N_ {b}}} = {\ frac {\ sigma _ {a}} {\ sigma _ {a} + \ sigma _ {b} }}}$

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

El algoritmo MISER procede al bisecar 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 posibles bisecciones y seleccionando la que minimizará 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. A continuación, se repite el mismo procedimiento de forma recursiva para cada uno de los dos medios espacios 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 simple de Monte Carlo. Estos valores individuales y sus estimaciones de error se combinan luego hacia arriba para dar un resultado general y una estimación de su error.

Muestreo de importancia

Existe una variedad de algoritmos de muestreo de importancia, como

Algoritmo de muestreo de importancia

El muestreo de importancia proporciona una herramienta muy importante para realizar la integración Monte-Carlo. El principal resultado del muestreo de 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 pueden elegir para reducir la varianza de la medición Q _N . ${\ Displaystyle {\ overline {\ mathbf {x}}}}$${\ Displaystyle p ({\ overline {\ mathbf {x}}})}$${\ Displaystyle p ({\ overline {\ mathbf {x}}})}$

Considere el siguiente ejemplo donde a uno le gustaría integrar numéricamente una función gaussiana, centrada en 0, con σ = 1, de −1000 a 1000. Naturalmente, si las muestras se extraen uniformemente en el intervalo [−1000, 1000], solo un una parte muy pequeña de ellos sería significativa para la integral. Esto se puede mejorar eligiendo una distribución diferente de la que se eligen 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

${\ Displaystyle p ({\ overline {\ mathbf {x}}}): \ qquad {\ overline {\ mathbf {x}}} _ {1}, \ cdots, {\ overline {\ mathbf {x}}} _ {N} \ en V,}$

el estimador de I viene dado por

${\ 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 en particular el doble que otras muestras, la ponderamos la mitad que las otras muestras. Este estimador es naturalmente válido para un muestreo uniforme, el caso donde es constante. ${\ Displaystyle p ({\ overline {\ mathbf {x}}})}$

El algoritmo Metropolis-Hastings es uno de los algoritmos más utilizados para generar a partir de , proporcionando así una forma eficiente de las integrales de computación. ${\ Displaystyle {\ overline {\ mathbf {x}}}}$${\ Displaystyle p ({\ overline {\ mathbf {x}}})}$

VEGAS Montecarlo

El algoritmo de 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. Para evitar que el número de intervalos de histograma crezca como K ^d , la distribución de probabilidad se aproxima mediante una función separable:

${\ Displaystyle g (x_ {1}, x_ {2}, \ ldots) = g_ {1} (x_ {1}) g_ {2} (x_ {2}) \ ldots}$

de modo que el número de bins necesarios es solo Kd . Esto equivale a ubicar 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 por importancia.

Ver también

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

Notas

Referencias

• Caflisch, RE (1998). "Métodos de Montecarlo y cuasi-Montecarlo". Acta Numerica . 7 : 1-49. Bibcode : 1998AcNum ... 7 .... 1C . doi : 10.1017 / S0962492900002804 .
• Weinzierl, S. (2000). "Introducción a los métodos de Monte Carlo". arXiv : hep-ph / 0006269 .
• Presione, WH; Farrar, GR (1990). "Muestreo estratificado recursivo para la integración multidimensional de Monte Carlo" . Computadoras en Física . 4 (2): 190. Código Bibliográfico : 1990ComPh ... 4..190P . doi : 10.1063 / 1.4822899 .
• Lepage, GP (1978). "Un nuevo algoritmo para la integración multidimensional adaptativa". Revista de Física Computacional . 27 (2): 192-203. Código Bibliográfico : 1978JCoPh..27..192L . doi : 10.1016 / 0021-9991 (78) 90004-9 .
• Lepage, GP (1980). "VEGAS: un programa de integración multidimensional adaptativo". Cornell Preprint CLNS 80-447 .
• Hammersley, JM; Handscomb, DC (1964). Métodos de Monte Carlo . Methuen. ISBN 978-0-416-52340-9.
• Presione, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007). Recetas numéricas: el arte de la informática científica (3ª ed.). Nueva York: Cambridge University Press. ISBN 978-0-521-88068-8.
• Newman, MEJ; Barkema, GT (1999). Métodos de Monte Carlo en Física Estadística . Prensa de Clarendon.
• Robert, CP; Casella, G (2004). Métodos estadísticos de Monte Carlo (2ª ed.). Saltador. 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 Naive Monte-Carlo: Documentación para las rutinas Naive Monte-Carlo de C ++
• Applet de Monte Carlo aplicado en problemas de física estadística