87.1. Convertir una cantidad en una esperanza

C86 terminó con una idea que ahora se vuelve método computacional. Si podemos escribir una cantidad como

θ=E[h(X)],\theta=\mathbb E[h(X)],

y podemos generar realizaciones

X1,…,XNX_1,\ldots,X_N

de la distribución de XX, entonces podemos aproximar θ\theta mediante

θ^N=1N∑i=1Nh(Xi).\widehat\theta_N = \frac1N \sum_{i=1}^N h(X_i).

Esta es la estructura elemental de Monte Carlo.

La ley de grandes números da la justificación básica:

θ^N→θ\widehat\theta_N \to \theta

bajo condiciones apropiadas.

El teorema central del límite añade una descripción aproximada del error cuando la varianza es finita:

Nθ^N−θσh⟶dN(0,1),\sqrt N \frac{ \widehat\theta_N-\theta }{ \sigma_h } \overset d\longrightarrow N(0,1),

donde

σh2=Var⁡(h(X)).\sigma_h^2 = \operatorname{Var}(h(X)).

Monte Carlo transforma un problema determinista o probabilístico difícil en un problema de promediar realizaciones.

La potencia del método proviene de que ciertas cantidades pueden representarse como expectativas y de que sabemos cómo se comportan los promedios, sin que el azar «encuentre» nada por sí mismo.

87.2. El estimador Monte Carlo elemental

Sea

Y=h(X)Y=h(X)

con

E[∣Y∣]<∞.\mathbb E[|Y|]<\infty.

Generamos realizaciones independientes e idénticamente distribuidas

Y1,…,YN.Y_1,\ldots,Y_N.

El estimador es

θ^N=YˉN.\widehat\theta_N = \bar Y_N.

Si

θ=E[Y],\theta=\mathbb E[Y],

entonces

E[θ^N]=θ.\mathbb E[\widehat\theta_N] = \theta.

El estimador es insesgado.

Si además

Var⁡(Y)=σY2<∞,\operatorname{Var}(Y)=\sigma_Y^2<\infty,
Var⁡(θ^N)=σY2N.\operatorname{Var}(\widehat\theta_N) = \frac{\sigma_Y^2}{N}.

Por tanto, su error estándar teórico es

SE⁡(θ^N)=σYN.\operatorname{SE}(\widehat\theta_N) = \frac{\sigma_Y}{\sqrt N}.

Esta tasa N−1/2N^{-1/2} es el sello del Monte Carlo ordinario independiente.

87.3. Duplicar precisión es caro

La tasa

SE⁡∝N−1/2\operatorname{SE} \propto N^{-1/2}

implica rendimientos decrecientes.

Para dividir el error estándar por dos necesitamos aproximadamente 4N4N realizaciones.

Para dividirlo por diez necesitamos aproximadamente 100N.100N.

Monte Carlo es robusto, pero no converge rápido.

Esto explica el papel de los métodos de reducción de varianza en problemas costosos. No los desarrollaremos aquí. Aparecerán más adelante en el currículo, cuando el objetivo sea mejorar la eficiencia sin alterar la cantidad que se estima.

C87 se concentrará en el Monte Carlo directo y en cómo auditarlo correctamente.

87.4. Integración como esperanza

Consideremos la integral

I=∫01e−x2 dx.I = \int_0^1 e^{-x^2}\,dx.

Si

U∼Unif⁡(0,1),U\sim\operatorname{Unif}(0,1),

su densidad es uno sobre el intervalo. Entonces

I=E[e−U2].I = \mathbb E[ e^{-U^2} ].

Definamos H=e−U2.H=e^{-U^2}.

Podemos aproximar la integral generando

U1,…,UNU_1,\ldots,U_N

uniformes independientes y calculando

I^N=1N∑i=1Ne−Ui2.\widehat I_N = \frac1N \sum_{i=1}^N e^{-U_i^2}.

El problema de integración se ha convertido en un promedio.

87.5. Valor exacto para poder auditar el método

Esta integral tiene una expresión conocida:

I=π2erf⁡(1)≈0.746824133.I = \frac{\sqrt\pi}{2} \operatorname{erf}(1) \approx 0.746824133.

El paquete reproducible usa este valor como referencia y lo contrasta además con cuadratura numérica de alta precisión.

La concordancia entre ambas es mejor que 10−12.10^{-12}.

Elegir un ejemplo con solución conocida permite separar dos preguntas, la de si el algoritmo Monte Carlo está implementado correctamente y la de qué comportamiento finito exhibe el estimador. Así evitamos confiar en el propio Monte Carlo para decidir si el Monte Carlo funciona.

87.6. Varianza exacta del integrando

El segundo momento es

E[H2]=∫01e−2x2 dx.\mathbb E[H^2] = \int_0^1 e^{-2x^2}\,dx.

Utilizando la función error,

E[H2]=π22erf⁡(2).\mathbb E[H^2] = \frac{\sqrt\pi}{2\sqrt2} \operatorname{erf}(\sqrt2).

Por tanto,

σH2=E[H2]−I2≈0.0403977.\sigma_H^2 = \mathbb E[H^2]-I^2 \approx 0.0403977.

La desviación estándar es aproximadamente σH≈0.200992.\sigma_H \approx 0.200992.

Así, antes de simular sabemos que

SE⁡(I^N)=0.200992N.\operatorname{SE}(\widehat I_N) = \frac{0.200992}{\sqrt N}.

La simulación tiene una predicción cuantitativa contra la cual compararse.

87.7. Una trayectoria Monte Carlo no es la convergencia teórica

La primera figura muestra la estimación acumulada

I^1,I^2,…,I^N\widehat I_1, \widehat I_2,\ldots,\widehat I_N

para una única secuencia pseudoaleatoria.

Convergencia acumulada de la integral.
Convergencia acumulada de la integral Ampliar figura ↗

La curva oscila alrededor del valor exacto y se estabiliza progresivamente.

Con

N=106,N=10^6,

la ejecución reproducible obtiene aproximadamente

I^N=0.746547314.\widehat I_N = 0.746547314.

El error de esta realización concreta es cercano a −2.77×10−4.-2.77\times10^{-4}.

Ese error es el resultado aleatorio de una ejecución. Leerlo como «el error del método para un millón de muestras» sería incorrecto.

Para estudiar la distribución del error necesitamos replicaciones independientes del experimento completo.

87.8. Repetir el propio Monte Carlo

El paquete ejecuta 2500 estimaciones independientes para cada uno de varios tamaños:

N=100, 300, 1000, 3000, 10000, 30000.N= 100,\ 300,\ 1000,\ 3000,\ 10000,\ 30000.

Para cada NN calcula la media de los estimadores, el sesgo empírico, la desviación estándar entre replicaciones, el RMSE, el error estándar teórico y la cobertura de un intervalo normal aproximado.

El estimador es insesgado en teoría. Los sesgos empíricos observados fluctúan alrededor de cero y disminuyen en importancia respecto de su desviación.

La desviación estándar entre replicaciones coincide estrechamente con

σH/N.\sigma_H/\sqrt N.

Así se verifica operacionalmente la teoría de C86.

87.9. RMSE y tasa de convergencia

Para un estimador

θ^,\widehat\theta,

el error cuadrático medio es

MSE⁡=E[(θ^−θ)2].\operatorname{MSE} = \mathbb E[ (\widehat\theta-\theta)^2 ].

Se descompone como

MSE⁡=Var⁡(θ^)+Bias⁡(θ^)2.\operatorname{MSE} = \operatorname{Var}(\widehat\theta) + \operatorname{Bias}(\widehat\theta)^2.

El RMSE es

RMSE⁡=MSE⁡.\operatorname{RMSE} = \sqrt{\operatorname{MSE}}.

Como el Monte Carlo elemental del ejemplo es insesgado,

RMSE⁡=σHN.\operatorname{RMSE} = \frac{\sigma_H}{\sqrt N}.

La regresión log-log de los RMSE empíricos produce una pendiente aproximada de

−0.496,-0.496,

muy próxima al valor teórico

−12.-\frac12.
RMSE Monte Carlo.
RMSE Monte Carlo Ampliar figura ↗

La simulación verifica que este experimento reproduce la consecuencia esperada de la ley general, sin demostrarla.

87.10. Estimar el error desde la propia muestra

En un problema real no conoceremos σY.\sigma_Y.

Podemos estimarla mediante la desviación estándar muestral

sY=1N−1∑i=1N(Yi−Yˉ)2.s_Y = \sqrt{ \frac1{N-1} \sum_{i=1}^N (Y_i-\bar Y)^2 }.

Entonces utilizamos

SE⁡^=sYN.\widehat{\operatorname{SE}} = \frac{s_Y}{\sqrt N}.

Bajo las hipótesis del TCL, un intervalo aproximado del 95 % es

θ^N±1.96SE⁡^.\widehat\theta_N \pm 1.96 \widehat{\operatorname{SE}}.

En las replicaciones del ejemplo, la cobertura observada queda aproximadamente entre 0.9460.946 y 0.953.0.953.

Eso es compatible con la aproximación nominal del 95 %.

87.11. Qué significa esa cobertura

Un intervalo Monte Carlo no expresa incertidumbre sobre todo el modelo.

Si repitiéramos el procedimiento de simulación muchas veces bajo la misma distribución y el mismo NN, aproximadamente el 95 % de esos intervalos contendría el valor objetivo en el régimen donde la aproximación normal funciona adecuadamente.

Quedan fuera de ese cálculo el error en el modelo probabilístico, la incertidumbre sobre parámetros externos, el sesgo de discretización, el error de un solucionador determinista interno y los errores de programación.

El error Monte Carlo es solo una capa de la incertidumbre total. Una simulación puede tener un error estándar minúsculo y representar con enorme precisión un modelo equivocado.

87.12. Error de modelo frente a error Monte Carlo

Supongamos que realmente queremos

θ=EP[h(X)],\theta = \mathbb E_P[h(X)],

pero simulamos desde una distribución incorrecta QQ.

Al aumentar NN,

θ^N→EQ[h(X)].\widehat\theta_N \to \mathbb E_Q[h(X)].

El error Monte Carlo respecto de la expectativa bajo QQ puede tender a cero mientras persiste

EQ[h(X)]−EP[h(X)].\mathbb E_Q[h(X)] - \mathbb E_P[h(X)].

Más simulaciones no corrigen ese sesgo estructural, que pertenece a la especificación y no al muestreo. Esta distinción es análoga a la de C80 entre precisión de una pendiente y validez del régimen de escala elegido.

87.13. Una semilla no genera aleatoriedad física

Los computadores ordinarios utilizan generadores pseudoaleatorios, algoritmos deterministas que producen secuencias diseñadas para comportarse estadísticamente como realizaciones de ciertas distribuciones.

Una semilla fija produce una secuencia reproducible dentro de la implementación especificada.

El paquete demuestra que dos generadores inicializados con la misma semilla producen los mismos primeros valores, mientras una semilla diferente produce otra secuencia.

La reproducibilidad es una ventaja científica.

Esa ventaja, sin embargo, exige conservar información suficiente sobre el algoritmo o biblioteca utilizados, su versión, la semilla o estado inicial y los parámetros de simulación.

La palabra «aleatorio» no debe utilizarse como excusa para resultados irreproducibles.

87.14. El error de reiniciar la misma semilla

Supongamos que queremos 100 replicaciones independientes.

Si dentro de cada repetición ejecutamos

rng = generator(seed=123)\texttt{rng = generator(seed=123)}

antes de generar datos, produciremos la misma realización 100 veces.

El promedio entre «replicaciones» duplica exactamente la misma información y no estima la variabilidad entre ejecuciones.

Una estrategia correcta consiste en mantener un flujo generador bien definido o crear subflujos independientes mediante mecanismos apropiados del generador.

El principio es el mismo que en C86, porque repetir datos idénticos no aumenta el tamaño efectivo de muestra.

87.15. Monte Carlo para probabilidades

Una probabilidad puede escribirse como esperanza de una variable indicadora.

Si

AA

es un evento,

P(A)=E[1A].P(A) = \mathbb E[ \mathbf 1_A ].

Por tanto, el estimador Monte Carlo es

p^=1N∑i=1N1A(Xi).\widehat p = \frac1N \sum_{i=1}^N \mathbf 1_A(X_i).

Es simplemente la fracción de realizaciones donde ocurre el evento.

Su varianza es

Var⁡(p^)=p(1−p)N.\operatorname{Var}(\widehat p) = \frac{p(1-p)}{N}.

El error estándar es

p(1−p)N.\sqrt{ \frac{p(1-p)}{N} }.

Esto funciona bien para probabilidades moderadas. Los eventos raros muestran una limitación seria.

87.16. El caso P(Z>3)P(Z>3)

Sea

Z∼N(0,1).Z\sim N(0,1).

Queremos estimar p=P(Z>3).p=P(Z>3).

El valor exacto es aproximadamente p=0.00134990.p=0.00134990.

El Monte Carlo ingenuo genera normales estándar y cuenta cuántas superan tres.

Con

N=10 000,N=10\,000,

esperamos solo Np≈13.5Np\approx13.5 aciertos.

En la ejecución reproducible aparecen ocho.

La estimación es 0.0008,0.0008, muy distinta en términos relativos del valor verdadero, aunque el algoritmo esté funcionando correctamente.

87.17. El error relativo de un evento raro

El error estándar relativo es

p(1−p)/Np.\frac{ \sqrt{p(1-p)/N} }{p}.

Puede escribirse como

1−pNp.\sqrt{ \frac{1-p}{Np} }.

Cuando

p≪1,p\ll1,
RSE⁡≈1Np.\operatorname{RSE} \approx \frac1{\sqrt{Np}}.

La cantidad importante es el número esperado de eventos Np.Np.

Si queremos un error estándar relativo del 10 %, necesitamos aproximadamente N≈7.4×104.N \approx 7.4\times10^4.

Para 3 %, N≈8.2×105.N \approx 8.2\times10^5.

Para 1 %, N≈7.4×106.N \approx 7.4\times10^6.

Costo de un evento raro.
Costo de un evento raro Ampliar figura ↗

El método funciona correctamente. Lo desfavorable es la varianza relativa de la cantidad que intentamos estimar.

87.18. Cero observaciones no significa probabilidad cero

Para eventos aún más raros puede ocurrir que ninguna de NN simulaciones alcance el evento.

Entonces

p^=0.\widehat p=0.

Eso solo indica que la muestra no observó éxitos, sin demostrar p=0.p=0.

La probabilidad de cero aciertos es (1−p)N.(1-p)^N.

Para pp pequeño, (1−p)N≈e−Np.(1-p)^N \approx e^{-Np}.

Si

Np≪1,Np\ll1,

observar cero es lo habitual incluso cuando el evento tiene probabilidad positiva.

Este punto será central en C88, donde estudiaremos colas y riesgo de extrapolación.

87.19. Monte Carlo multidimensional

Consideremos

Id=∫[0,1]dexp⁡(−∑j=1dxj)dx.I_d = \int_{[0,1]^d} \exp\left( -\sum_{j=1}^d x_j \right) dx.

La integral es separable:

Id=∏j=1d∫01e−xj dxj.I_d = \prod_{j=1}^d \int_0^1e^{-x_j}\,dx_j.

Por tanto, Id=(1−e−1)d.I_d = (1-e^{-1})^d.

El valor exacto permite auditar Monte Carlo en varias dimensiones.

Si

U∼Unif⁡([0,1]d),U\sim\operatorname{Unif}([0,1]^d),

definimos

Hd(U)=e−∑jUj.H_d(U) = e^{-\sum_jU_j}.

Entonces

Id=E[Hd(U)].I_d = \mathbb E[H_d(U)].

El mismo estimador medio funciona para cualquier dd.

87.20. La tasa formal no depende de la dimensión

Para Monte Carlo independiente,

SE⁡=Var⁡(Hd)N.\operatorname{SE} = \sqrt{ \frac{ \operatorname{Var}(H_d) }{N} }.

El exponente de NN sigue siendo

−12-\frac12

independientemente de dd.

Esta propiedad hace atractivo a Monte Carlo en problemas de alta dimensión.

Una regla de cuadratura construida como una malla con mm puntos por eje requiere mdm^d evaluaciones.

Con cinco puntos por eje y d=8,d=8, ya obtenemos 58=3906255^8=390625 puntos.

Con diez puntos por eje,

108=100 000 000.10^8=100\,000\,000.

El crecimiento tensorial puede ser prohibitivo.

87.21. Dimensión independiente no significa costo independiente

La frase «la tasa Monte Carlo no depende de la dimensión» debe interpretarse con precisión.

El factor N−1/2N^{-1/2} no cambia.

La varianza

Var⁡(Hd)\operatorname{Var}(H_d)

sí puede depender de dd, en cambio. También puede aumentar el costo de evaluar una sola realización.

En el ejemplo reproducible, con N=200000,N=200000, el error estándar relativo crece desde aproximadamente 6.4×10−46.4\times10^{-4} en una dimensión hasta 4.4×10−34.4\times10^{-3} en veinte dimensiones.

Efecto de la dimensión.
Efecto de la dimensión Ampliar figura ↗

Lo que Monte Carlo evita es la explosión combinatoria directa de una malla tensorial, sin conceder inmunidad universal frente a la dimensionalidad.

87.22. Este ejemplo multidimensional es deliberadamente fácil

La integral anterior es separable.

Eso significa que el mejor método no sería Monte Carlo ni una malla de cien millones de puntos. Podemos calcularla mediante una integral unidimensional y elevar el resultado a la potencia dd.

El ejemplo aísla dos propiedades, sin pretender que Monte Carlo sea superior. Una malla tensorial ingenua crece como mdm^d, mientras Monte Carlo conserva su tasa N−1/2N^{-1/2} aunque su constante de error pueda cambiar.

La estructura del problema debe explotarse siempre que exista. Utilizar Monte Carlo donde una derivación exacta resuelve el problema sería desperdiciar información.

87.23. Monte Carlo frente a cuadratura determinista

En una dimensión, para una función suave como e−x2,e^{-x^2}, los métodos deterministas de cuadratura pueden converger muchísimo más rápido que Monte Carlo.

La cuadratura utilizada como control en este capítulo reproduce el valor de referencia prácticamente a precisión de máquina con un costo minúsculo frente a un millón de muestras.

Ambos hechos son compatibles. Monte Carlo resulta atractivo cuando la dimensión es alta, cuando el dominio es complicado, cuando la expectativa proviene de un simulador, cuando la distribución permite muestrear pero no integrar con facilidad o cuando el integrando no se presta a cuadratura estructurada. No existe un ganador universal entre integración determinista y estocástica.

87.24. Una simulación puede ser más costosa que el muestreo

Hasta ahora evaluar h(X)h(X) ha sido casi gratuito.

En aplicaciones reales una realización puede requerir resolver un sistema de ecuaciones diferenciales, ejecutar un modelo basado en agentes o simular una red grande.

Entonces el costo total es aproximadamente

N×costo por realizacioˊn.N\times \text{costo por realización}.

Reducir el error por un factor dos exige cuatro veces más muestras y, en primera aproximación, cuatro veces más tiempo computacional.

Esta relación convierte la varianza en una cantidad operacional. Menor varianza del estimador significa menos simulaciones para una precisión dada.

87.25. Error de discretización dentro de Monte Carlo

Supongamos que cada realización necesita resolver numéricamente una ecuación diferencial con paso

Δt.\Delta t.

Entonces podemos estar estimando

θΔt=E[h(XΔt)],\theta_{\Delta t} = \mathbb E[ h(X_{\Delta t}) ],

no exactamente

θ=E[h(X)].\theta = \mathbb E[ h(X) ].

El error total puede descomponerse conceptualmente en

θ^N,Δt−θ=θ^N,Δt−θΔt⏟error Monte Carlo+θΔt−θ⏟error de discretizacioˊn.\widehat\theta_{N,\Delta t} - \theta = \underbrace{ \widehat\theta_{N,\Delta t} - \theta_{\Delta t} }_{\text{error Monte Carlo}} + \underbrace{ \theta_{\Delta t}-\theta }_{\text{error de discretización}}.

Aumentar NN reduce solo el primer término.

Disminuir Δt\Delta t puede reducir el segundo, pero encarece cada realización.

Una simulación reproducible debe distinguir ambas fuentes de error.

87.26. Detener la simulación según un criterio de precisión

Un procedimiento razonable puede definir de antemano una tolerancia Monte Carlo.

Por ejemplo, queremos

SE⁡^≤τ.\widehat{\operatorname{SE}} \leq \tau.

Como

SE⁡^=sN,\widehat{\operatorname{SE}} = \frac{s}{\sqrt N},

podemos aumentar NN hasta alcanzar la tolerancia.

Sin embargo, el criterio se vuelve aleatorio porque ss también se estima.

En problemas ordinarios esto puede ser suficiente como regla operacional. En problemas más delicados se requieren procedimientos secuenciales con garantías explícitas.

La idea general es que el número de simulaciones debería fijarse a partir de la precisión requerida, antes que como una cifra redonda sin diagnóstico de error.

87.27. No reutilizar el mismo error estándar para todas las cantidades

Una simulación puede estimar varias salidas:

θ1,θ2,…\theta_1,\theta_2,\ldots

Cada una puede tener varianza diferente.

Una cantidad central puede estimarse con gran precisión mientras una probabilidad de cola permanece extremadamente ruidosa.

El hecho de haber ejecutado N=106N=10^6 realizaciones no define por sí solo la calidad de todas las estimaciones producidas.

Debe calcularse el error Monte Carlo de cada cantidad relevante.

Para probabilidades raras, cuantiles extremos o razones de cantidades pequeñas, el comportamiento puede ser mucho peor que para una media central.

87.28. Reproducibilidad y separación entre experimento y análisis

Una implementación robusta debería separar:

  1. generación de realizaciones;
  2. cálculo de la cantidad h(X)h(X);
  3. agregación del estimador;
  4. diagnóstico de incertidumbre;
  5. producción de figuras y tablas.

Esta separación facilita verificar cada etapa.

El paquete reproducible del capítulo conserva datos fuente resumidos, semilla, resultados de replicaciones, cálculos de error y verificaciones analíticas.

La reproducibilidad exige que sea posible reconstruir qué se estimó, con qué distribución y con qué tolerancia, más allá de que el código «corra».

87.29. Pseudorandomness no sustituye validación estadística

Una secuencia pseudoaleatoria puede ser perfectamente reproducible y aun así utilizarse mal.

Entre los errores frecuentes están correlacionar involuntariamente flujos que deberían ser independientes, reiniciar semillas, transformar mal uniformes a otra distribución, reutilizar realizaciones sin reconocer dependencia y descartar valores extremos porque parecen errores.

La calidad del generador no corrige un algoritmo probabilístico equivocado.

Del mismo modo, obtener cifras distintas al cambiar de semilla solo muestra que existe variabilidad entre ejecuciones, sin acreditar que el simulador sea «realista».

87.30. Qué significa aprender mediante realizaciones repetidas

Monte Carlo puede interpretarse como aprendizaje computacional sobre una distribución especificada.

Cada realización aporta información sobre

E[h(X)].\mathbb E[h(X)].

La media acumulada reduce el ruido de muestreo.

El método, sin embargo, aprende dentro del modelo. No descubre por sí solo qué distribución debería usarse, qué variables se omitieron o si el simulador representa el fenómeno real.

La cadena lógica es

modelo probabilıˊstico→realizaciones→estimador→error Monte Carlo.\text{modelo probabilístico} \to \text{realizaciones} \to \text{estimador} \to \text{error Monte Carlo}.

La validez del primer eslabón sigue requiriendo argumentos externos a la simulación.

87.31. La distribución de muestreo forma parte del estimador

La integral unidimensional utilizó una uniforme porque el dominio era [0,1][0,1]:

I=∫01h(x) dx=E[h(U)].I = \int_0^1 h(x)\,dx = \mathbb E[h(U)].

Si el dominio D⊂RdD\subset\mathbb R^d tiene volumen finito ∣D∣|D| y

U∼Unif⁡(D),U\sim\operatorname{Unif}(D),

entonces

∫Dh(x) dx=∣D∣ E[h(U)].\int_D h(x)\,dx = |D|\, \mathbb E[h(U)].

El factor ∣D∣|D| no puede omitirse salvo que el volumen sea uno.

Más generalmente, si queremos calcular

θ=∫h(x)p(x) dx,\theta = \int h(x)p(x)\,dx,

donde pp es una densidad normalizada, podemos escribir

θ=Ep[h(X)].\theta = \mathbb E_p[h(X)].

El subíndice recuerda respecto de qué distribución se toma la esperanza.

La elección de la distribución de muestreo y la definición de la función promediada forman una pareja. Cambiar una sin modificar la otra cambia la cantidad estimada.

Existe una identidad más general. Si generamos desde otra densidad q(x)q(x), positiva donde h(x)p(x)h(x)p(x) es relevante,

θ=∫h(x)p(x)q(x)q(x) dx=Eq[h(X)p(X)q(X)].\theta = \int h(x) \frac{p(x)}{q(x)} q(x)\,dx = \mathbb E_q \left[ h(X)\frac{p(X)}{q(X)} \right].

Esta expresión es la base del muestreo de importancia. No desarrollaremos aquí cómo elegir qq ni sus diagnósticos, porque el currículo lo retomará en C235. El punto de C87 es conceptual. Monte Carlo consiste en construir una expectativa bajo una distribución de muestreo conocida y promediar el integrando correcto, lo que rara vez equivale a «muestrear uniformemente».

Una implementación puede generar millones de puntos sin estimar la cantidad deseada si los pesos, el volumen del dominio o la distribución de referencia están mal especificados.

87.32. Insesgado no se conserva bajo cualquier posprocesamiento

El estimador de la integral satisface

E[I^N]=I.\mathbb E[\widehat I_N]=I.

Eso no significa que cualquier función posterior del estimador siga siendo insesgada.

Supongamos que queremos informar

g(I)g(I)

pero calculamos

g(I^N).g(\widehat I_N).

En general,

E[g(I^N)]≠g(I).\mathbb E[g(\widehat I_N)] \neq g(I).

Por ejemplo, si

g(x)=x2,g(x)=x^2,

entonces

E[I^N2]=Var⁡(I^N)+I2.\mathbb E[\widehat I_N^2] = \operatorname{Var}(\widehat I_N) + I^2.

Por tanto,

E[I^N2]>I2\mathbb E[\widehat I_N^2] > I^2

si la varianza es positiva.

El sesgo es exactamente

σH2N.\frac{\sigma_H^2}{N}.

Desaparece cuando N→∞N\to\infty, pero existe para tamaño finito.

Esta observación cuenta en simulaciones donde la cantidad final es una razón, un logaritmo, un máximo o alguna transformación no lineal de varias expectativas estimadas.

La ley de grandes números puede garantizar consistencia mediante continuidad de gg, pero consistencia e insesgamiento finito no son sinónimos.

También recuerda una idea de C85:

E[g(X)]≠g(E[X]).\mathbb E[g(X)] \neq g(\mathbb E[X]).

Monte Carlo aproxima expectativas. Aplicar después una transformación cambia el objeto algebraico y puede introducir otro componente del error.

87.33. Planificar el tamaño Monte Carlo

Si conocemos aproximadamente la desviación estándar

σY\sigma_Y

y queremos que un semiancho normal aproximado sea como máximo δ\delta,

z1−α/2σYN≤δ.z_{1-\alpha/2} \frac{\sigma_Y}{\sqrt N} \leq \delta.

Despejando,

N≥(z1−α/2σYδ)2.N \geq \left( \frac{ z_{1-\alpha/2}\sigma_Y }{ \delta } \right)^2.

Para un intervalo aproximado del 95 %, z0.975≈1.96.z_{0.975}\approx1.96.

Si σY\sigma_Y es desconocida, puede estimarse mediante una ejecución piloto.

Esta planificación contiene incertidumbre porque el propio sYs_Y es aleatorio. Aun así, es preferible a elegir NN solo porque una potencia de diez resulta cómoda.

También hace visible el costo cuadrático de la precisión. Si reducimos δ\delta a la mitad, NN se multiplica aproximadamente por cuatro.

Para eventos raros conviene utilizar precisión relativa. Como

RSE⁡(p^)=1−pNp,\operatorname{RSE}(\widehat p) = \sqrt{\frac{1-p}{Np}},

resolver

RSE⁡≤r\operatorname{RSE}\leq r

produce

N≥1−ppr2.N \geq \frac{1-p}{pr^2}.

El cálculo del capítulo para P(Z>3)P(Z>3) proviene exactamente de esta identidad.

87.34. Replicaciones, lotes y diagnóstico de estabilidad

Una única media acumulada puede ocultar problemas.

Dividir las realizaciones en lotes independientes permite observar si las estimaciones parciales son compatibles con el error esperado. Si cuatro lotes enormes producen resultados separados por muchas desviaciones estándar, puede existir dependencia inadvertida, un cambio de régimen, un error de implementación o una subestimación de la varianza.

En el Monte Carlo i.i.d. elemental, si cada lote contiene mm realizaciones,

θ^(1),…,θ^(B)\widehat\theta^{(1)},\ldots, \widehat\theta^{(B)}

deberían comportarse como estimadores aproximadamente independientes con varianza σY2/m.\sigma_Y^2/m.

La variabilidad entre lotes aporta una comprobación externa al cálculo de sYs_Y dentro de una única secuencia.

No deben confundirse lotes con repetir la misma secuencia debido a una semilla reiniciada. Los lotes necesitan información pseudoaleatoria distinta.

El diseño por lotes gana peso cuando las muestras no son independientes, como ocurrirá posteriormente con MCMC. Allí las observaciones consecutivas pueden estar autocorrelacionadas y la fórmula

s/Ns/\sqrt N

subestimaría la incertidumbre si se aplica ingenuamente.

C87 no desarrolla ese caso, pero fija un principio que sobrevivirá. El tamaño nominal de una simulación no determina el tamaño informativo sin conocer la dependencia entre realizaciones.

87.35. Verificar un simulador con problemas que sí sabemos resolver

Monte Carlo suele utilizarse cuando la respuesta exacta es desconocida. Eso no elimina la posibilidad de verificar partes del procedimiento.

Una estrategia robusta consiste en construir primero casos de referencia donde la respuesta sí sea conocida. El ejemplo de la integral

∫01e−x2 dx\int_0^1 e^{-x^2}\,dx

cumple esa función. También lo hace la integral multidimensional separable.

Si el código falla en esos problemas, no existe fundamento para confiar en una aplicación más compleja donde el valor verdadero está oculto.

La verificación puede realizarse por capas. El generador se contrasta mediante distribuciones conocidas. Las transformaciones se comprueban contra identidades analíticas. El estimador se audita usando medias o probabilidades exactas. La tasa de convergencia se compara con

N−1/2.N^{-1/2}.

Cuando el simulador interno incluye dinámica, pueden utilizarse invariantes, soluciones especiales o límites simplificados, siguiendo el método de C29.

Esto separa verificación de validación.

Verificar pregunta si resolvemos correctamente el modelo que escribimos. Validar pregunta si ese modelo representa suficientemente el fenómeno de interés.

Monte Carlo puede aprobar todas las pruebas de verificación y seguir siendo científicamente inadecuado porque las distribuciones, parámetros o mecanismos elegidos no corresponden al sistema real. Por eso conviene conservar, junto con el resultado final, una jerarquía explícita de controles que cubra la exactitud del generador, la corrección de las transformaciones, la consistencia del estimador, la estabilidad frente a semillas y el acuerdo con casos límite. Cuanto más opaco sea el simulador, más importante resulta disponer de problemas reducidos donde una parte de la respuesta pueda calcularse por otra vía independiente.

Esta distinción pesa más cuando la simulación produce imágenes o distribuciones visualmente convincentes. La apariencia de variabilidad realista no sustituye una prueba de que cada capa del algoritmo implementa la estructura declarada. En simulación científica, la trazabilidad de esas capas es parte del resultado y no una tarea administrativa posterior.

87.36. Verificación computacional

El paquete reproducible ejecuta controles en cuatro niveles.

Para la integral

I=∫01e−x2 dx,I=\int_0^1e^{-x^2}\,dx,

comprueba que la expresión con erf⁡\operatorname{erf} coincide con cuadratura de alta precisión. Deriva el segundo momento y la varianza del integrando. La estimación con un millón de muestras permanece dentro de cinco errores estándar teóricos.

Después ejecuta 2500 replicaciones para seis tamaños muestrales. La desviación estándar empírica de los estimadores coincide con

σH/N,\sigma_H/\sqrt N,

la cobertura de intervalos normales aproximados permanece cerca de 95 % y la pendiente log-log del RMSE es −0.496,-0.496, compatible con la tasa teórica −1/2-1/2.

Para el evento P(Z>3),P(Z>3), verifica la fórmula del error estándar relativo y el crecimiento cuadrático del tamaño requerido al exigir mayor precisión relativa. Con un millón de realizaciones aparecen 1370 aciertos y la estimación queda dentro de cuatro errores estándar del valor exacto.

En el ejemplo multidimensional se calculan integrales exactas para dimensiones entre 1 y 20 y se verifica que todas las estimaciones queden dentro de cinco errores estándar. También se documenta el crecimiento mdm^d de una malla tensorial.

Finalmente, dos generadores inicializados con la misma semilla reproducen exactamente la misma secuencia, mientras una semilla distinta produce otra.

Se ejecutaron 25 comprobaciones automáticas y las 25 resultaron satisfactorias.

87.37. Qué queda establecido

Monte Carlo convierte expectativas en promedios de realizaciones.

La ley de grandes números justifica convergencia. El TCL clásico explica la escala aproximada del error cuando la varianza es finita.

El estimador elemental tiene error estándar

σ/N,\sigma/\sqrt N,

por lo que mejorar precisión es costoso.

El método puede ser competitivo en dimensiones altas porque el exponente de convergencia respecto de NN no depende directamente de la dimensión, pero la constante de error y el costo por realización sí pueden hacerlo.

Los eventos raros revelan otra limitación, porque si solo una fracción pp de las simulaciones aporta un acierto la precisión relativa depende aproximadamente de NpNp.

Finalmente, el error Monte Carlo debe separarse de error del modelo, discretización y programación.

El siguiente capítulo, C88, estudiará las colas amplias, los cuantiles y el riesgo de extrapolación. Allí veremos por qué las regiones menos visitadas por Monte Carlo ordinario pueden ser las decisivas para ciertas preguntas y por qué resumir incertidumbre mediante media y varianza puede resultar insuficiente.

Referencias

Recursos reproducibles

Descargar recursos del capítulo (ZIP, 287 kB).

6 archivos CSV · 4 figuras · 1 programa Python.

El programa se ha ejecutado en una carpeta vacía y genera los archivos de datos y figuras incluidos en el paquete.

Ver instrucciones y contenido del paquete.