8.2: Simulación Montecarlo

Supongamos que hacemos\(N\) simulaciones, cada vez dibujando los parámetros aleatorios necesarios\(x_i\) a partir de un número aleatorio” caja negra” (sobre el cual daremos más detalles en la siguiente sección). Definimos la salida de alto nivel de nuestro sistema\(S\) para ser\(g(\vec{x})\). Por simplicidad, decimos que\(g(\vec{x})\) es un escalar. \(g(\vec{x})\)puede ser prácticamente cualquier salida de interés, por ejemplo: el valor de un estado en un momento dado después de una entrada impulsiva, o la integral en el tiempo de la trayectoria de una de las salidas, con una entrada dada. En lo que sigue, dejaremos caer la notación vectorial\(x\) para mayor claridad.

Dejar que el estimador\(G\) de\(g(x)\) ser definido como

\[ G = \dfrac{1}{N} \sum_{j=1}^{N} g(x_j). \]

Es posible que reconozcas esto como un promedio directo. En efecto, tomando la expectativa de ambas partes,

\[ E(G) = \dfrac{1}{N} \sum_{j=1}^{N} E(g(x_i)), \]

es claro que\(E(G) = E(g)\). Al mismo tiempo, sin embargo, no lo sabemos\(E(g)\); calculamos\(G\) entendiendo que con un número muy grande de ensayos,\(G\) deben acercarse\(E(g)\). Ahora veamos la varianza del estimador. Esto conceptualmente resulta de un número infinito de ensayos con estimadores, cada uno de los cuales implica\(N\) evaluaciones de\(g\) acuerdo con la definición anterior. Es importante tener en cuenta que tal varianza involucra muestras del estimador (cada una involucrando\(N\) evaluaciones), no de la función subyacente\(g(x)\). Tenemos

\ begin {align}\ sigma^2 (G)\, &=\,\ sigma^2\ izquierda [\ dfrac {1} {N}\ sum_ {j=1} ^N g (x_j)\ derecha]\\ [4pt] &=\,\ dfrac {1} {N^2}\ sigma^2\ izquierda [\ sum_ {j=1} ^N g (x_j)\ derecha]\\ [4pt] &=\,\ dfrac {1} {N^2}\ suma_ {j=1} ^N\ sigma^2 (g)\\ [4pt] &=\,\ dfrac {1} {N}\ sigma^2 (g). \ end {align}

Esta relación es clave. La primera igualdad se deriva del hecho de que\(\sigma^2 (nx) = n^2 \sigma^2 (x)\), si\(n\) es una constante. La segunda igualdad es verdadera porque\(\sigma^2 (x+y) = \sigma^2 (x) + \sigma^2 (y)\), donde\(x\) y\(y\) son variables aleatorias. El resultado principal es que\(\sigma^2 (G) = \sigma^2 (g)\) si sólo se consideran ensayos de una muestra, pero eso\(\sigma^2 (G) \rightarrow 0\) como\(N \rightarrow \infty\). De ahí que con un tamaño suficientemente grande\(N\) podemos esperar que nuestro\(G\) va a estar muy cerca de\(E(g)\).

Llevemos esto un poco más lejos, para obtener una estimación explícita del error en a\(G\) medida que vamos a lo grande\(N\). Definir un error de estimador no dimensional

\ begin {align} q\, &=\,\ dfrac {G - E (g)} {\ sigma (G)}\\ [4pt] &=\,\ dfrac {(G - E (g))\ sqrt {N}} {\ sigma (g)},\ end {align}

donde la segunda igualdad proviene del resultado anterior. Llamamos al factor\(\sigma(g) / \sqrt{N}\) el error estándar. Invocando el teorema del límite central, que garantiza que la distribución de G se vuelve gaussiana para N suficientemente grande, tenemos

\ begin {align}\ lim_ {N\ a\ infty}\ textrm {prob} (a<q<b)\, &=\,\ int\ limits_ {a} ^ {b}\ dfrac {1} {\ sqrt {2\ pi}} e^ {-t^2/2}\, dt\\ [4pt] &=\, F (a) - F (b),\ end {align}

donde\(F(x)\) es la función de probabilidad acumulativa de la variable gaussiana estándar:

\[ F(a) = \int\limits_{- \infty}^{a} \dfrac{1}{\sqrt{2 \pi}} e^{-t^2 / 2} \, dt \]

Al buscar algunos valores para\(F(x)\), vemos que el error no dimensional es menor de uno en 68.3% de los ensayos; es menos de dos en 95.4% de los ensayos y menos de tres en 99.7% de los ensayos. El intervalo de confianza del 99.7% corresponde con

\ begin {align} -3\, &\ leq\, (G - E (g))\ sqrt {N}/\ sigma (g)\,\ leq\, 3\,\ fila derecha\\ [4pt] -3\ sigma (g)/\ sqrt {N}\, &\ leq\,\ quad\,\, G - E (g)\,\,\ quad\ quad\,\ leq\, 3\ sigma (g)/\ sqrt {N}. \ end {align}

En general, cuadruplicar el número de ensayos mejora el error en un factor de dos.

Hasta el momento hemos estado describiendo un solo estimador\(G\), que recupera la media. La media, sin embargo, es de hecho una integral sobre el dominio aleatorio:

\[ E(g) \, = \int\limits_{x \epsilon X} p(x) g(x) \, dx, \]

donde\(p(x)\) esta el pdf de la variable aleatoria\(x\). Entonces el estimador de Monte Carlo\(G\) es de hecho un integrador:

\[ G \, \simeq \int\limits_{x \epsilon X} p(x) g(x) \, dx. \]

Podemos definir fácilmente estimadores de momentos estadísticos:

\[ G_n \, = \, \dfrac{1}{N} \sum_{j=1}^{N} x_j^n g(x_j) \, \simeq \int\limits_{x \exists X} x^n p(x) g(x) \, dx, \]que seguirá las mismas tendencias básicas de convergencia del estimador de medias\(G\). Estos momentos se pueden calcular todos utilizando las mismas\(N\) evaluaciones de\(g(x)\).

La ecuación anterior da otro punto de vista para entender cómo funciona el enfoque de Monte Carlo: el efecto de la función de densidad de probabilidad en la integral es reemplazado por el hecho de que las variables aleatorias en MC se extraen de la misma distribución. En otras palabras,\(g(x)\) ahí se\(X\) amplifica un alto\(p(x)\) en un área determinada del dominio. MC hace lo mismo, porque de hecho hay más\(x\) sacados de esta área, en hacer las\(N\) evaluaciones.