Erreur et convergence de Monte-Carlo

De combien l'estimation peut-elle se tromper ?

Une estimation de Monte Carlo est la moyenne de n tirages aléatoires, donc elle est elle-même aléatoire. Son erreur typique est l'erreur standard : l'écart type des tirages divisé par sqrt(n), la même formule que dans le test t. Deux conséquences comptent. Première : l'erreur ne dépend pas de la difficulté apparente du problème, seulement de la dispersion de la quantité que vous moyennez. Seconde : elle diminue lentement : pour gagner un chiffre de précision de plus, il faut 100 fois plus de tirages.

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.00216

Barres d'erreur et coût de la précision

Donnez un nombre simulé avec son intervalle de confiance, moyenne -/+ 1,96 * erreur standard : la plage qui contient la vraie valeur environ 95 % du temps. Si vous avez besoin de la réponse à une tolérance donnée, inversez la formule de l'erreur standard : n = (std / erreur_cible)², arrondi au supérieur. Avec une dispersion d'environ 2,2 et une cible de 0,001, il faut presque cinq millions de tirages ; pour 0,01, seulement 47 000. Lancez d'abord une petite simulation pilote pour connaître std, puis décidez du nombre de tirages de la vraie exécution.

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 draws

Une précision moins chère : les variables antithétiques

Outre tirer davantage, on peut tirer plus malin. Pour une loi normale centrée réduite Z, l'image miroir -Z est un tirage tout aussi valable. Si f est monotone, f(Z) et f(-Z) varient en sens opposés, donc leur moyenne est bien moins bruitée que chacune. L'estimateur antithétique fait la moyenne de (f(z) + f(-z)) / 2 sur les tirages. Dans l'exemple, l'erreur standard passe de 0,0217 à 0,0119 pour le même nombre de nombres aléatoires, et pour une fonction impaire comme 3z + 1 le bruit disparaît complètement. C'est l'une des astuces de réduction de variance.

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]

• Erreur standard = std / sqrt(n) : 100 fois plus de tirages donnent 10 fois moins d'erreur ; donnez mean -/+ 1,96 * se.
• Tirages nécessaires pour une erreur cible : ceil((std / cible) ** 2) ; lancez d'abord un petit pilote pour connaître std.
• Variables antithétiques : moyennez (f(z) + f(-z)) / 2 pour réduire le bruit sans nombre aléatoire supplémentaire ; vérifiez toujours la convergence avec une autre graine.

À suivre : la valeur actuelle et l'actualisation, premier pas vers les taux : un euro demain vaut moins qu'un euro aujourd'hui.

Lab : ce qui est noté

• Mesurer le bruit d'une estimation avec standard_error25 pts
• Mettre des barres d'erreur sur la réponse avec mc_confidence_interval25 pts
• Budgétiser la simulation avec samples_needed25 pts
• Obtenir la même précision pour moins cher avec antithetic_mean25 pts