אינטגרציית מונטה קרלו

אנימציה של אינטגרציית מונטה קרלו של הפולינום $x^3-6x^2+9x+1$ בתחום $[0,4]$. הקווים הנוספים מסמלים את כיסוי התנהגות הפונקציה ע״י הוספת נקודות דגימה נוספות. בכותרת מוצגת השגיאה היחסית שקטנה ככל שמספר נקודות הדגימה עולה.

אינטגרציית מונטה קרלו היא שיטת אינטגרציה נומרית מבוססת דגימה אקראית של נקודות בתחום האינטגרל, והיא מקרה פרטי של שיטת מונטה קרלו. השיטה שימושית במיוחד לחישוב אינטגרלים רב ממדיים מאחר שאינה מחלקת את תחום האינטגרל לרשת בדידה שגדלה בצורה מעריכית כתלות בממד.

שיטה זאת זוכה לשימושים נרחבים במתמטיקה פיננסית, פיזיקה סטטיסטית, גרפיקה ממוחשבת ובעוד תחומים רבים.

רקע

נניח כי ברצוננו לחשב את האינטגרל המסוים $I$ של פונקציה חד או רב ממדית $g(𝐱)$ על התחום $\mathrm{\Omega }\subset {ℝ}^n$ כאשר $\mathrm{\Omega }$ היא קבוצה סופית.

$$I=\underset{\mathrm{\Omega }}{\int }g(𝐱)d𝐱$$נניח כי $\mathrm{\Omega }$ בעלת נפח סופי $V=\underset{\mathrm{\Omega }}{\int }d𝐱$.

נגדיר משתנה מקרי $\overrightarrow{X}\sim \text{Unif}(\mathrm{\Omega })$ שמתפלג בצורה אחידה על $\mathrm{\Omega }$ ונשים לב כי פונקציית הצפיפות $f$ של $\overrightarrow{X}$ תוגדר להיות$$f(x)=\{\begin{array}{cc}p=1/V& x\in \mathrm{\Omega }\\ 0& \text{otherwise}\end{array}$$

נשים לב כי, התוחלת של $\overrightarrow{X}$ כאשר עליו מופעלת הפונקציה $g$ היא

$$𝔼[g(\overrightarrow{X})]=\underset{\overline{\mathrm{\Omega }}}{\int }g(𝐱)f(𝐱)d𝐱=\underset{\overline{\mathrm{\Omega }}}{\int }g(𝐱)pd𝐱=pI=I/V$$

מאידך, עבור $N$ דגימות מהמ״מ $\overrightarrow{X}$, ${𝐱}_1,{𝐱}_2,\dots ,{𝐱}_N\in \overline{\mathrm{\Omega }}$ נשערך את $𝔼[g(\overrightarrow{X})]$ ע״י

$$𝔼[g(\overrightarrow{X})]\approx \frac{1}{N}\sum _{i=1}^{N}g({𝐱}_i)$$ולבסוף נקבל

$$I=V\cdot 𝔼[g(\overrightarrow{X})]\approx Q_N\equiv \frac{V}{N}\sum _{i=1}^{N}g({𝐱}_i)$$

מחוק המספרים הגדולים נקבל

$$\underset{N\to \mathrm{\infty }}{\lim }{Q}_N=V\cdot 𝔼[g(\overrightarrow{X})]=I$$

הערכת השגיאה

ראשית, נשים לב כי $Q_N$ משערך חסר הטיה תבנית:אנ

$$𝔼[Q_N]=𝔼[\frac{V}{N}\sum _{i=1}^{N}g({𝐱}_i)]=\frac{V}{N}\sum _{i=1}^N𝔼[g({𝐱}_i)]=V\cdot 𝔼[g(\overrightarrow{X})]=I$$

נחשב את המשערך חסר ההטייה ${\hat{\sigma }}_N^2$ של השונות $𝕍[g(\overrightarrow{X})]$

$${\hat{\sigma }}_N^2\equiv \frac{1}{N-1}\sum _{i=1}^N(g({𝐱}_i)-I{)}^2$$

שגיאה יחסית של אינטגרציית מונטה קרלו כתלות במספר הדגימות $N$ שדועכת בקצב התכנסות $1/\sqrt{N}$

והשונות של המשערך

$Q_N$

$$𝕍[Q_N]=𝕍[\frac{V}{N}\sum _{i=1}^{N}g({𝐱}_i)]=\frac{V^2}{N^2}\sum _{i=1}^N𝕍[g({𝐱}_i)]=\frac{V^2}{N}\cdot 𝕍[g(\overrightarrow{X})]\approx V^2\frac{{\hat{\sigma }}_N^2}{N}$$

ולבסוף נקבל, כי שורש השגיאה הריבועית ממוצעת, ע״פ משפט 6.9 של והיעדר ההטייה של $Q_N$ הוא

$$\text{RMSE}(Q_N)=\sqrt{𝕍[Q_N]+\text{bias}(Q_N)}=\sqrt{𝕍[Q_N]}\approx V\frac{{\sigma }_N^2}{\sqrt{N}}$$

דוגמה: שערוך הערך $\pi$ ע״י חישוב שטח מעגל היחידה

אנימציה של שערוך הערך $\mathrm{\pi }$ ע״י דגימת נקודות

הדגמה נפוצה של השיטה היא שערוך הערך $\pi$ ע״י חישוב שטח מעגל היחידה ${𝕊}^1$. על אף שזו אינה דרך יעילה לחשיוב הערך של $\pi$, היא מדגימה את השיטה.

בשביל לחשב את $\pi$ ניתן לחשב אינטגרל על התחום $[-1,1]\times [-1,1]$ של פונקציית האינדיקטור של מעגל היחידה ${\mathrm{𝟏}}_{S^1}$ המוגדרת להיות

$${\mathrm{𝟏}}_{S^1}(x,y)=\{\begin{array}{cc}1& \text{if }{x}^2+y^2\le 1\\ 0& \text{else}\end{array}$$ומכך מקבלים:

$$\pi =\text{area}(S^1)=\underset{S^1}{\int }d𝐱=\underset{[-1,1{]}^2}{\int }{\mathrm{𝟏}}_{S^1}d𝐱$$

והמשערך יהיה

$$Q_N(\pi )=\frac{V}{N}\sum _{i=1}^{N}g({𝐱}_i)=\frac{\text{vol}([-1,1{]}^2)}{N}\sum _{i=1}^N{\mathrm{𝟏}}_{S^1}=4\cdot \frac{\mathrm{\#}\{\text{samples in }{S}^1\}}{\mathrm{\#}\{\text{total samples}\}}$$

כלומר למעשה, החישוב נעשה ע״י דגימת נקודות בריבוע היחידה $[-1,1{]}^2$, ספירה של הנקודות שנופלות בתוך מעגל היחידה, חלוקה במספר נקודות הדגימה הכולל (תחולת האינטגרנד) והכפלה בנפח ריבוע היחידה.

מימוש ההדגמה בשפת פייתון:

import numpy as np

### parameters
# number of samples
N = 1e6

# generate N uniform random points inside the unit cube [-1, 1] x [-1, 1]
X = np.random.uniform(low=-1, high=1, size=(int(N), 2))
domain_volume = 4  # unit cube area.

# mask for points inside the unit circle - equivalent to applying the indicator function of the unit circle.
distance = np.sqrt(X[:, 0] ** 2 + X[:, 1] ** 2)
inside_mask = distance <= 1  # array of booleans.

# compute the mean of the indicator function over the domain and multiply by the domain volume to get the integral value.
indicator_mean = np.sum(inside_mask) / N
integral_value = indicator_mean * domain_volume

ניתן כמובן לשערך את הערך

$\pi$

בעזרת חישוב שטחם של מעגלים בעלי רדיוס שונה מ-1

יתרונות וחסרונות

תבנית:תמונות מרובות יתרונות

• פשטות - להבנה ומימוש - השיטה אינטואיטיבית ודי פשוטה למימוש כפי שניתן לראות בקוד לעיל.
• היעדר תלות בממדתבנית:הערה - שיטות אינטגרציה נומרית אחרות (תרבוע גאוסיאני ושיטת סימפסוןתבנית:אנ לדוגמה) מחלקות את תחום האינטגרציה לרשת בדידה בעלת $N$ נקודות, בהן הפונקציה מופעלת. כאשר ממד האינטגרל גדול מאחד $d\ge 1$ אנו מקבלים גידול אקספוננציאלי במספר הנקודות מאחר שעלינו לחלק את התחום לרשת בדידה בכל ממד כלומר - $N^d$ נקודות. שיטת מונטה לעומת זאת, דוגמת $N$ נקודות במרחב הרב ממדי כולו ולכן אינה סובלת מ״קללת הממד״תבנית:אנ

חסרונות

• התכנסות איטית - כאמור, השגיאה הריבועית ממוצעת דועכת בקצב $O(\frac{1}{\sqrt{N}})$ כתלות במספר הדגימות $N$.
• אינה דטרמיניסטית - בעלת מרכיב אקראי (דגימת הנקודות) - בשונה משיטות אינטגרציה נומרית רשתיות כגון שיטת הטרפזתבנית:אנ שמחלקות את תחום האינטגרציה לרשת בדידה דטרמיניסטית, שיטת מונטה קרלו דוגמת אקראית נקודות בתחום האינטגרציה. צורת הדגימה של הנקודות חשובה מאוד. לדוגמה, לא נרצה לקבל מעיין ״גוש״ נקודות קרובות אחת לשנייה (איור 1), אלא דגימה אחידה במרחב ע״מ לקבל כיסוי מספיק של התנהגות הפונקציה בתחום. נוסף על כך, אקראיות השיטה תאלץ אותנו לחזור על חישובנו מספר רב של פעמים על מנת ליצור סטטיסטיקה מספקת על התוצאות ולהשתמש ברווח בר-סמך ו/או מבחנים סטטיסטים.

וריאציות

דגימת שכבות רקורסיבית

הכללה של אינטגרציה מסתגלתתבנית:אנ לשיטת מונטה קרלו. כאשר בכל חזרור, אם הדיוק המתקבל אינו עומד בדיוק הרצוי, תחום האינטגרציה מחולק ל-2 תתי תחומים, אשר עליהם נעשית האינטגרציה, בצורה איטרטיבית. החלוקה תעשה לתחומים בהם שונות הפונקציה היא הגדולה ביותר על מנת להביא לכיסוי טוב יותר של התחום באזורים חשובים אלה.

כאשר האינטגרל הוא רב ממדי, חלוקה של כל ממד ל-2 תגדיל את מספר הנקודות בצורה מעריכית, לכן החלוקה תעשה בממד בו חלוקה נוספת תביא להקטנת השגיאה המשמעותית ביותר.

ניתן להרחיב על דגימת שכבות רקורסיבית ואלגוריתם MISERתבנית:אנ המממש אותה.

דגימת חשיבות

שגיאה יחסית של חישוב ערך האינטגרל של פעמון גאוסיאני סטנדרטי בתחום האינטגרציה $[-50,50]$ כתלות במספר הדגימות $N$. מדגיש את עליונות אלגוריתם וגאס תבנית:אנ אשר דוגם פחות נקודות בקצוות התחום (מאחר שחשיבותן נמוכה מאוד) לעומת דגימה אחידה על התחום האינטגרציה.

בדומה לאלגוריתמים דטרמיניסטים כגון שיטת סימפסוןתבנית:אנ, נקודות בעלות ערך שונה בהרבה מסביבתן, אשר משפיעות על הערך הסופי של האינטגרל, עלולות ״להתפספס״ כנקודות דגימה ולהשפיע לרעה על התוצאה הסופית. לכן ישנה חשיבות רבה על נקודות הדגימה עצמן, הדרך בהן הן נבחרות.

השיטה ״הנאבית״ פשוטה להבנה ומימוש, הנקודות נבחרות בצורה שרירותית ע״י דגימה אחידה על תחום האינטגרציה, ללא התייחסות לחשיבותן בתוצאה הסופית. וריאציות של השיטה מסוג "דגימת חשיבות" מציעות צורת דגימה אחרת.

דוגמה טובה לכך היא חישוב האינטגרל המסוים של פעמון גאוסיאני סטנדרטי בעל תוחלת $\mu =0$ וסטיית תקן $\sigma =1$ על התחום $[-50,50]$. כאשר נשתמש בשיטה הנאיבית, מאחר שהדגימה אחידה על התחום, יתקבלו נקודות דגימה רחוקות ממרכז הפעמון ($0$) למרות היותן בעלות חשיבות נמוכה עד זניחה עבור חישוב האינטגרל שהרי עבור $|x|\ge 3$ (3 סטיות תקן) ערך הפונקציה קטן מ- $5\cdot 10^{-3}$. אי לכך, רצוי להתמקד בדגימה באזורים בהם הנקודות בעלות חשיבות גבוהה (סביב מרכז הפעמון בדוגמה לעיל). רעיון זה ממומש, לדוגמה, באלגוריתם וגאס תבנית:אנ.

מימוש בקוד

קוד פייתון כללי כאשר הפונקציה ותחום האינטגרציה ידועים, והתחום הוא מלבן רב ממדי:

import numpy as np
from typing import Callable


def monte_carlo_integration(
    domain: np.ndarray, N: int, f: Callable[[np.ndarray], float]
) -> float:
    """monte_carlo_integration

    Args:
        domain (np.ndarray): d-dimensional domain - [a_i, b_i] x ... x [a_d, b_d] - represented as a 2 x d numpy array.
        N (int): number of random points to generate.
        f (Callable[[np.ndarray], float]): function to integrate that receives a d-dimensional point and returns a float.

    Returns:
        float: integral value of f over the domain.
    """

    # generate N uniform random points - a random matrix of size N x d with each row representing a random point.
    # first column of the domain holds the lower bounds and the second column holds the upper bounds.
    X = np.random.uniform(
        low=domain[:, 0], high=domain[:, 1], size=(N, domain.shape[1])
    )

    # evaluate the function at the random points
    f_values = np.apply_along_axis(func1d=f, axis=1, arr=X)

    # compute the domain volume and final integral
    domain_volume = np.prod(domain[:, 1] - domain[:, 0])
    integral_value = np.mean(f_values) * domain_volume

    return integral_value


# Example usage with unit circle
domain = np.array([[-1, 1], [-1, 1]])
N = 1e6
f = lambda arr: np.where(np.linalg.norm(arr) <= 1, 1, 0)
integral_value = monte_carlo_integration(domain, N, f)

בנוסף, קיימות דוגמאותתבנית:הערה וספריות פייתון שממשות גרסאות שונות של האלגוריתם.

קישורים חיצוניים

• תבנית:MathWorld