Quanto a estimativa pode errar?
Uma estimativa Monte Carlo é a média de n sorteios aleatórios, então ela própria é aleatória. Seu erro típico é o erro padrão: o desvio padrão dos sorteios dividido por sqrt(n), a mesma fórmula do teste t. Duas consequências importam. Primeira: o erro não depende da dificuldade aparente do problema, só da dispersão da quantidade que você tira a média. Segunda: ele diminui devagar: para ganhar mais um dígito de precisão são necessários 100 vezes mais sorteios.
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 erro e custo da precisão
Informe um número simulado com seu intervalo de confiança, média -/+ 1,96 * erro padrão: a faixa que contém o valor verdadeiro cerca de 95% das vezes. Se precisar da resposta dentro de uma tolerância, inverta a fórmula do erro padrão: n = (std / erro_alvo)², arredondado para cima. Com dispersão de cerca de 2,2 e alvo de 0,001 são necessários quase cinco milhões de sorteios; para 0,01, só 47 000. Rode primeiro uma pequena simulação piloto para conhecer std e depois decida quantos sorteios a execução real precisa.
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 drawsPrecisão mais barata: variáveis antitéticas
Além de sortear mais, é possível sortear com mais inteligência. Para uma normal padrão Z, a imagem espelhada -Z é um sorteio igualmente válido. Se f é monótona, f(Z) e f(-Z) se movem em sentidos opostos, então a média deles é muito menos ruidosa que cada uma. O estimador antitético tira a média de (f(z) + f(-z)) / 2 sobre os sorteios. No exemplo o erro padrão cai de 0,0217 para 0,0119 com o mesmo número de números aleatórios, e para uma função ímpar como 3z + 1 o ruído some por completo. É um de vários truques de redução de variância.
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]- Erro padrão =
std / sqrt(n): 100 vezes mais sorteios dão 10 vezes menos erro; informemean -/+ 1,96 * se. - Sorteios necessários para um erro-alvo:
ceil((std / alvo) ** 2); rode antes um piloto pequeno para conhecerstd. - Variáveis antitéticas: tire a média de
(f(z) + f(-z)) / 2para reduzir o ruído sem números aleatórios extras; sempre verifique a convergência com outra semente.
A seguir: valor presente e desconto, o primeiro passo na renda fixa: um euro amanhã vale menos que um euro hoje.
Lab: o que é avaliado
- Medir o ruído de uma estimativa com
standard_error25 pts - Pôr barras de erro na resposta com
mc_confidence_interval25 pts - Orçar a simulação com
samples_needed25 pts - Obter a mesma precisão por menos com
antithetic_mean25 pts