¿Cuánto puede equivocarse la estimación?
Una estimación Monte Carlo es la media de n extracciones aleatorias, así que es aleatoria. Su error típico es el error estándar: la desviación típica de las extracciones dividida por sqrt(n), la misma fórmula que en la prueba t. Dos consecuencias importan. Primera: el error no depende de lo difícil que parezca el problema, solo de la dispersión de la cantidad que promedias. Segunda: se reduce despacio: para ganar un dígito más de precisión hacen falta 100 veces más extracciones.
import math
import numpy as np
# Estimate E[exp(Z)] for a standard normal Z; the exact value is exp(0.5) = 1.6487
for n in (100, 10000, 1000000):
values = np.exp(np.random.default_rng(5).standard_normal(n))
print(n, round(values.mean(), 4), round(values.std(ddof=1) / math.sqrt(n), 5))
# 100 1.2167 0.14583
# 10000 1.6887 0.0217 the error falls roughly like 1 / sqrt(n)
# 1000000 1.6508 0.00216Barras de error y coste de la precisión
Informa un número simulado con su intervalo de confianza, media -/+ 1,96 * error estándar: el rango que contiene el valor verdadero alrededor del 95 % de las veces. Si necesitas la respuesta dentro de una tolerancia, invierte la fórmula del error estándar: n = (std / error_objetivo)², redondeado hacia arriba. Con una dispersión de unos 2,2 y un objetivo de 0,001 hacen falta casi cinco millones de extracciones; para 0,01, solo 47 000. Ejecuta primero una pequeña simulación piloto para conocer std y luego decide cuántas extracciones necesita la ejecución real.
values = np.exp(np.random.default_rng(5).standard_normal(10000))
mean, std = values.mean(), values.std(ddof=1)
std_error = std / math.sqrt(len(values))
print(round(mean - 1.96 * std_error, 4), round(mean + 1.96 * std_error, 4)) # 1.6462 1.7312 contains the exact 1.6487
print(math.ceil((std / 0.01) ** 2)) # 47111 draws for a standard error of 0.01
print(math.ceil((std / 0.001) ** 2)) # 4710946 10 times the precision, 100 times the drawsPrecisión más barata: variables antitéticas
Además de extraer más, se puede extraer con más cabeza. Para una normal estándar Z, la imagen especular -Z es una extracción igual de válida. Si f es monótona, f(Z) y f(-Z) se mueven en sentidos opuestos, así que su promedio es mucho menos ruidoso que cada una. El estimador antitético promedia (f(z) + f(-z)) / 2 sobre las extracciones. En el ejemplo el error estándar baja de 0,0217 a 0,0119 con el mismo número de números aleatorios, y para una función impar como 3z + 1 el ruido desaparece por completo. Es uno de varios trucos de reducción de varianza.
z = np.random.default_rng(5).standard_normal(10000)
plain = np.exp(z)
antithetic = (np.exp(z) + np.exp(-z)) / 2 # each draw is paired with its mirror image
print(round(plain.mean(), 4), round(plain.std(ddof=1) / 100, 5)) # 1.6887 0.0217
print(round(antithetic.mean(), 4), round(antithetic.std(ddof=1) / 100, 5)) # 1.6565 0.01192 half the noise, same 10,000 numbers
print(round(math.exp(0.5), 4)) # 1.6487 the exact value# The running mean of one long simulation settles down as n grows
values = np.exp(np.random.default_rng(5).standard_normal(100000))
running = np.cumsum(values) / np.arange(1, len(values) + 1)
print([round(float(running[k - 1]), 4) for k in (10, 100, 1000, 100000)]) # [1.5522, 1.2167, 1.634, 1.6502]- Error estándar =
std / sqrt(n): 100 veces más extracciones dan 10 veces menos error; informamean -/+ 1,96 * se. - Extracciones necesarias para un error objetivo:
ceil((std / objetivo) ** 2); ejecuta primero un piloto pequeño para conocerstd. - Variables antitéticas: promedia
(f(z) + f(-z)) / 2para reducir el ruido sin números aleatorios extra; comprueba siempre la convergencia con otra semilla.
A continuación: el valor actual y el descuento, el primer paso hacia la renta fija: un euro mañana vale menos que un euro hoy.
Lab: qué se califica
- Medir el ruido de una estimación con
standard_error25 pts - Poner barras de error a la respuesta con
mc_confidence_interval25 pts - Presupuestar la simulación con
samples_needed25 pts - Obtener la misma precisión con menos con
antithetic_mean25 pts