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

קובץ:Monte-carlo-integration-animation-of-polynomial.gif

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

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

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

רקע

נניח כי ברצוננו לחשב את האינטגרל המסוים <math>I</math> של פונקציה חד או רב ממדית <math>g(\mathbf{x})</math> על התחום <math>\Omega \sub \mathbb{R}^n</math> כאשר <math>\Omega</math> היא קבוצה סופית.

<math display="block">I = \int_{\Omega} g( \mathbf{x}) d \mathbf{x}</math>נניח כי <math>\Omega</math> בעלת נפח סופי <math>V = \int_\Omega d \mathbf{x}</math>.

נגדיר משתנה מקרי <math>\vec X \sim \text{Unif}(\Omega)</math> שמתפלג בצורה אחידה על <math>\Omega</math> ונשים לב כי פונקציית הצפיפות <math>f</math> של <math>\vec X</math> תוגדר להיות<math display="block">f(x) = \begin{cases} p = 1/V & x \in \Omega \\ 0 &\text{otherwise} \end{cases}</math>

נשים לב כי, התוחלת של <math>\vec X</math> כאשר עליו מופעלת הפונקציה <math>g</math> היא

<math display="block">\mathbb{E}[g(\vec X)] = \int_{\bar \Omega} g(\mathbf{x}) f(\mathbf{x}) d\mathbf{x} = \int_{\bar \Omega} g(\mathbf{x}) p \ d\mathbf{x} = p I = I/V</math>

מאידך, עבור <math>N</math> דגימות מהמ״מ <math>\vec X</math>, <math>\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_N \in \bar \Omega</math> נשערך את <math>\mathbb{E}[g(\vec X)]</math> ע״י

<math display="block">\mathbb{E}[g(\vec X)] \approx {1 \over N} \sum_{i=1}^N g(\mathbf{x}_i)</math>ולבסוף נקבל

<math display="block">I = V \cdot \mathbb{E}[g(\vec X)] \approx Q_N \equiv {V \over N} \sum_{i=1}^N g(\mathbf{x}_i) </math>

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

<math display="block">\lim_{N \to \infty} Q_N = V \cdot \mathbb{E}[g(\vec X)] = I</math>

הערכת השגיאה

ראשית, נשים לב כי <math>Q_N</math> משערך חסר הטיה

<math display="block">\mathbb{E}[Q_N] = \mathbb{E} \Big[{V \over N } \sum_{i=1}^N g(\mathbf{x}_i) \Big] = {V \over N} \sum_{i=1}^N \mathbb{E} [g(\mathbf{x}_i)] = V \cdot \mathbb{E}[g(\vec X)] = I</math>

נחשב את המשערך חסר ההטייה <math>\hat \sigma^2_N</math> של השונות <math>\mathbb{V}[g(\vec X)]</math>

<math display="block">\hat \sigma^2_N \equiv {1 \over N-1} \sum_{i=1}^N (g(\mathbf{x}_i) - I)^2 </math>

קובץ:Error of a Monte Carlo integration to calculate pi.svg

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

והשונות של המשערך <math>Q_N</math>

<math display="block">\mathbb{V} [Q_N] = \mathbb{V} \Big[{V \over N} \sum_{i=1}^N g(\mathbf{x}_i) \Big] = {V^2 \over N^2} \sum_{i=1}^N \mathbb{V} [g(\mathbf{x}_i)] = {V^2 \over N} \cdot \mathbb{V}[g(\vec X)] \approx V^2 {\hat \sigma^2_N \over N}</math>

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

<math display="block">\text{RMSE}(Q_N) = \sqrt{\mathbb{V} [Q_N] + \text{bias}(Q_N)} = \sqrt{\mathbb{V} [Q_N]} \approx V {\sigma^2_N \over \sqrt{N}}</math>

דוגמה: שערוך הערך <math>\pi</math> ע״י חישוב שטח מעגל היחידה

קובץ:Estimacion de Pi por Montercarlo.gif

אנימציה של שערוך הערך <chem>\pi</chem> ע״י דגימת נקודות

הדגמה נפוצה של השיטה היא שערוך הערך <math>\pi</math> ע״י חישוב שטח מעגל היחידה <math>\mathbb{S}^1</math>. על אף שזו אינה דרך יעילה לחשיוב הערך של <math>\pi</math>, היא מדגימה את השיטה.

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

<math display="block">\mathbf{1}_{S^1} (x, y) = \begin{cases} 1 & \text{if }x^{2}+y^{2}\leq1\\ 0 & \text{else} \end{cases}</math>ומכך מקבלים:

<math display="block">\pi = \text{area}( S^1) = \int_{S^1} d \mathbf{x} = \int_{[-1, 1]^2} \mathbf{1}_{S^1} d \mathbf{x}</math>

והמשערך יהיה

<math display="block">Q_N(\pi) = {V \over N} \sum_{i=1}^N g(\mathbf{x}_i) = {\text{vol} ([-1, 1]^2) \over N} \sum_{i=1}^N \mathbf{1}_{S^1} = 4 \cdot { \# \{ \text{samples in } S^1 \} \over \# \{ \text{total samples} \}} </math>

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

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

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

ניתן כמובן לשערך את הערך <math>\pi</math> בעזרת חישוב שטחם של מעגלים בעלי רדיוס שונה מ-1

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

קובץ:Pseudorandom sequence 2D.svg

איור א׳ - 256 דגימות אחידות פסודו-אקראיות

קובץ:Halton sequence 2D.svg

256 הנקודות הראשונות ברצף הלטון 2,3

יתרונות

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

חסרונות

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

וריאציות

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

קובץ:Postscript-viewer-blue.svg ערך מורחב – דגימת שכבות

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

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

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

דגימת חשיבות

קובץ:Vegas-vs-naive-standard-normal-curve-on-50-50-relative-error-animation.gif

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

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

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

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

מימוש בקוד

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

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 (באנגלית)