خطأ مونت كارلو وتقاربها

كم يمكن أن يخطئ التقدير؟

تقدير مونت كارلو متوسط n سحبة عشوائية، فهو نفسه عشوائي. وخطؤه المعتاد هو الخطأ المعياري: الانحراف المعياري للسحوبات مقسومًا على sqrt(n)، وهي الصيغة نفسها في اختبار t. ونتيجتان مهمتان. الأولى: لا يعتمد الخطأ على صعوبة المسألة الظاهرة بل على تشتت الكمية التي تأخذ متوسطها. الثانية: أنه يتناقص ببطء: لكسب خانة دقة إضافية تحتاج إلى 100 ضعف من السحوبات.

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

أشرطة الخطأ وكلفة الدقة

اذكر العدد المحاكى مع فاصل الثقة المتوسط -/+ 1.96 * الخطأ المعياري: النطاق الذي يحوي القيمة الحقيقية نحو 95% من المرات. وإن أردت الجواب ضمن تفاوت معين فاعكس صيغة الخطأ المعياري: n = (std / الخطأ_المستهدف)² مع التقريب لأعلى. ومع تشتت نحو 2.2 وهدف 0.001 تحتاج إلى قرابة خمسة ملايين سحبة؛ ولهدف 0.01 تكفي 47000. شغّل محاكاة تجريبية صغيرة أولًا لمعرفة std ثم قرّر عدد سحوبات التشغيل الحقيقي.

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

دقة أرخص: المتغيرات المتعاكسة

إلى جانب السحب أكثر يمكنك السحب بذكاء أكبر. فبالنسبة لمتغير طبيعي معياري Z تكون الصورة المرآتية -Z سحبة صالحة تمامًا. وإذا كانت f رتيبة تتحرك f(Z) وf(-Z) في اتجاهين متعاكسين فيكون متوسطهما أقل ضجيجًا بكثير من كلٍّ منهما. والمقدّر المتعاكس يأخذ متوسط (f(z) + f(-z)) / 2 على السحوبات. وفي المثال ينخفض الخطأ المعياري من 0.0217 إلى 0.0119 بالعدد نفسه من الأعداد العشوائية، وفي دالة فردية مثل 3z + 1 يختفي الضجيج تمامًا. وهذه إحدى حيل تقليل التباين.

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]

• الخطأ المعياري = std / sqrt(n): 100 ضعف السحوبات يعطي خطأ أقل بعشرة أضعاف؛ اذكر mean -/+ 1.96 * se.
• السحوبات اللازمة لخطأ مستهدف: ceil((std / target) ** 2)؛ شغّل تجربة صغيرة أولًا لمعرفة std.
• المتغيرات المتعاكسة: خذ متوسط (f(z) + f(-z)) / 2 لتقليل الضجيج دون أعداد عشوائية إضافية؛ وتحقق دائمًا من التقارب ببذرة أخرى.

التالي: القيمة الحالية والخصم، الخطوة الأولى نحو الدخل الثابت: يورو غدًا يساوي أقل من يورو اليوم.

المختبر: ما الذي يُقيَّم

• قِس ضجيج التقدير باستخدام standard_error25 نقطة
• ضع أشرطة خطأ على الجواب باستخدام mc_confidence_interval25 نقطة
• ضع ميزانية المحاكاة باستخدام samples_needed25 نقطة
• احصل على الدقة نفسها بكلفة أقل باستخدام antithetic_mean25 نقطة