Monte-Carlo-Fehler und Konvergenz

Wie falsch kann die Schätzung sein?

Eine Monte-Carlo-Schätzung ist der Mittelwert von n Zufallsziehungen, also selbst zufällig. Ihr typischer Fehler ist der Standardfehler: die Standardabweichung der Ziehungen geteilt durch sqrt(n), dieselbe Formel wie beim t-Test. Zwei Folgen zählen. Erstens hängt der Fehler nicht davon ab, wie schwierig das Problem aussieht, sondern nur von der Streuung der gemittelten Größe. Zweitens schrumpft er langsam: Für eine Stelle mehr Genauigkeit braucht man 100-mal mehr Ziehungen.

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

Fehlerbalken und Kosten der Genauigkeit

Geben Sie eine simulierte Zahl mit ihrem Konfidenzintervall an, Mittel -/+ 1,96 * Standardfehler: dem Bereich, der den wahren Wert in etwa 95 % der Fälle enthält. Brauchen Sie die Antwort innerhalb einer Toleranz, kehren Sie die Standardfehlerformel um: n = (std / Zielfehler)², aufgerundet. Bei einer Streuung von etwa 2,2 und einem Ziel von 0,001 sind fast fünf Millionen Ziehungen nötig; für 0,01 nur 47.000. Laufen Sie zuerst eine kleine Pilotsimulation, um std kennenzulernen, und entscheiden Sie dann, wie viele Ziehungen der eigentliche Lauf braucht.

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

Billigere Genauigkeit: antithetische Variablen

Statt mehr zu ziehen, kann man klüger ziehen. Bei einer Standardnormalen Z ist das Spiegelbild -Z eine ebenso gültige Ziehung. Ist f monoton, bewegen sich f(Z) und f(-Z) gegenläufig, ihr Mittel ist also viel weniger verrauscht als jedes einzelne. Der antithetische Schätzer mittelt (f(z) + f(-z)) / 2 über die Ziehungen. Im Beispiel sinkt der Standardfehler bei gleicher Zahl von Zufallszahlen von 0,0217 auf 0,0119, und bei einer ungeraden Funktion wie 3z + 1 verschwindet das Rauschen ganz. Das ist einer von mehreren Tricks zur Varianzreduktion.

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]

• Standardfehler = std / sqrt(n): 100-mal mehr Ziehungen ergeben 10-mal weniger Fehler; geben Sie mean -/+ 1,96 * se an.
• Nötige Ziehungen für einen Zielfehler: ceil((std / Ziel) ** 2); laufen Sie zuerst einen kleinen Pilot, um std zu kennen.
• Antithetische Variablen: Mitteln Sie (f(z) + f(-z)) / 2, um das Rauschen ohne zusätzliche Zufallszahlen zu senken; prüfen Sie die Konvergenz stets mit einem anderen Seed.

Als Nächstes: Barwert und Diskontierung, der erste Schritt in die festverzinslichen Papiere: Ein Euro morgen ist weniger wert als ein Euro heute.

Lab: Was bewertet wird

• Mit standard_error das Rauschen einer Schätzung messen25 Pkt.
• Mit mc_confidence_interval Fehlerbalken an die Antwort setzen25 Pkt.
• Mit samples_needed die Simulation budgetieren25 Pkt.
• Mit antithetic_mean dieselbe Genauigkeit billiger erreichen25 Pkt.