شبیه سازی مونت کارلو در پایتون — راهنمای کاربردی

یکی از شیوه‌های شبیه‌سازی که برای حل بسیاری از مسائل بهینه‌سازی و محاسبه عددی انتگرال به کار می‌رود، «شبیه سازی مونت کارلو» (Monte Carlo Simulation) است. الگوریتم‌هایی که برمبنای این شیوه شبیه‌سازی پایه‌ریزی شده‌اند براساس نمونه‌گیری از توزیع‌های احتمال و آماری عمل می‌کنند. حیطه به کارگیری شبیه‌سازی مونت کارلو از حوزه محاسبات مالی تا تئوری‌های فیزیک گسترده است. برای مثال به منظور شبیه‌سازی اتفاقاتی که در یک رآکتور هسته‌ای رخ می‌دهد، می‌توان از شبیه‌سازی مونت‌کارلو استفاده کرد.

در این مطلب قصد داریم شیوه شبیه‌سازی مونت کارلو را معرفی کنیم. برای اجرای الگوریتم مربوطه نیز از زبان برنامه نویسی پایتون استفاده خواهیم کرد. به این منظور برای آشنایی با روش مونت کارلو بهتر است مطلب متد مونت‌کارلو – به زبان ساده و شبیه سازی مونت کارلو (Monte Carlo Simulation) – محاسبه انتگرال‌ به روش عددی را مطالعه کنید. از طرفی اطلاع از شیوه‌های نمونه‌گیری که در نوشتار روش‌ های نمونه‌گیری (Sampling) در آمار — به زبان ساده آمده است نیز خالی از لطف نیست. همچنین مطالعه مطلب توزیع های آماری — مجموعه مقالات جامع وبلاگ فرادرس نیز به منظور آگاهی از توزیع‌ها آماری و کاربردشان توصیه می‌شود.

شبیه سازی مونت کارلو

همانطور که گفته شد یکی از کاربردهای شبیه‌سازی مونت‌کارلو، محاسبه انتگرال به صورت عددی است. برای مثال فرض کنید قرار است انتگرال زیر را محاسبه کنیم.

$$\large I= \int_0^{\infty} \dfrac{e^{-x}}{1+(x-1)^2}dx$$

رابطه ۱

همانطور که مشخص است این انتگرال را نمی‌توان به کمک روش‌های تحلیلی محاسبه کرد. بنابراین استفاده از روش‌های عددی در اینجا ضروری به نظر می‌رسد. از همین رو در ادامه به شیوه محاسبه تقریبی انتگرال به کمک شبیه‌سازی مونت‌کارلو می‌پردازیم.

فیلم آموزش ابزار مونت کارلوی Geant4 – تکمیلی در فرادرس

کلیک کنید

محاسبه انتگرال به شیوه مونت‌کارلو

تکنیک مونت‌کارلو ساده (Crude Monte Carlo Simulation) یکی از ساده‌ترین تکنیک‌ها در گروه روش‌های شبیه‌سازی مونت‌کارلو است. البته نباید از این شیوه انتظار داشت که در انجام محاسبات دقت زیادی به خرج دهد. براساس قضیه مقدار میانگین مشخص است که برای تابع پیوسته و مشتق‌پذیر $$f$$ می‌توان رابطه زیر را نوشت:

$$\large f_{ave}=\dfrac{1}{b-a}\int_a^bf(x)dx$$

به این ترتیب به نظر می‌رسد که مقدار انتگرال را می‌توان به صورت ضرب میانگین مقدار انتگرال‌ده در فاصله بازه انتگرال محاسبه کرد.

$$\large f_{ave}\times (b-a) =\int_a^bf(x)dx$$

دقت کنید که در اینجا منظور از $$f_{ave}$$ همان میانگین مقدار انتگرال‌ده در بازه $$[a,b]$$‌است. تکنیک شبیه‌سازی مونت‌کارلو دقیقا بر همین اصل و قضیه مقدار میانگین عمل می‌کند. بنابراین به جای محاسبه انتگرال به صورت تحلیلی، از روش تقریبی به کمک میانگین‌گیری تابع انتگرال‌ده استفاده می‌کنیم. البته باید توجه داشت که در بعضی از موارد شیوه‌های انتگرال‌گیری عددی بوسیله روش‌های هندسی مانند روش نیوتن-رافسون نیز صورت می‌گیرد.

در شیوه محاسبه انتگرال با شبیه‌سازی مونت‌کارلو، مقادیری از مقدار متغیر ($$x$$) را براساس تولید اعداد تصادفی در بازه قابل تعریف برای انتگرال (در اینجا فاصله $$a$$ تا $$b$$) ایجاد می‌کنیم. برای مثال فرض کنید تابع $$f$$ به شکل ساده $$y=2x$$ باشد. بنابراین برای محاسبه انتگرال زیر در بازه $$[0,2]$$ چندین عدد تصادفی برای متغیر $$x$$ در فاصله $$[0,2]$$ ایجاد کرده و مقدار $$y=۲x$$ را محاسبه می‌کنیم. انتگرال موردنظر تقریباً برابر با میانگین مقدارهای $$y$$ خواهد بود.

با توجه به تصویر بالا، مقدار انتگرال که توسط میانگین مقدارهای $$y$$ حاصل می‌شود برابر با $$1.984$$ است. در حقیقت می‌دانیم مقدار واقعی این انتگرال توسط رابطه زیر محاسبه می‌شود.

$$\large \int_0^2 2xdx= \dfrac{2x^2}{2}|^{2}_{0}=4-0=4$$

حال به مسئله اصلی بر می‌گردیم که محاسبه انتگرال رابطه ۱ بود. برای انجام محاسبات در این مرحله از کدهای پایتون استفاده می‌کنیم.

1import numpy as np
2import math
3import random
4from matplotlib import pyplot as plt
5from IPython.display import clear_output
6
7e = 2.71828

همانطور که مشاهده می‌کنید ابتدا کتابخانه‌های مورد نیاز بارگذاری شده‌اند. از کتابخانه math برای تعریف تابع و از random‌ نیز به منظور تولید نمونه تصادفی استفاده خواهیم کرد. همچنین از کتابخانه matplotlib نیز برای ایجاد و نمایش نمودارها استفاده می‌شود.

دستورات بعدی به منظور تولید اعداد تصادفی ایجاد شده‌اند. توجه داشته باشید که دو پارامتر min_value و max_value پارامترهای توزیع یکنواخت را مشخص کرده تا محدوده اعداد تصادفی تعیین شوند.

1def get_rand_number(min_value, max_value):
2    """
3    This function gets a random number from a uniform distribution between
4    the two input values [min_value, max_value] inclusively
5    Args:
6    - min_value (float)
7    - max_value (float)
8    Return:
9    - Random number between this range (float)
10    """
11    range = max_value - min_value
12    choice = random.uniform(0,1)
13    return min_value + range*choice

در اینجا از توزیع احتمالی یکنواخت برای تولید اعداد تصادفی استفاده شد ولی ممکن است با توجه به شکل تغییرات تابع $$f(x)$$ از توزیع‌های دیگر آماری برای تولید اعداد تصادفی استفاده کرد که با تابع $$f(x)$$ شبیه‌تر باشند. این نکته را در بخش «اهمیت نمونه‌گیری» شرح خواهیم داد.

از طرفی برای تعریف تابع محاسباتی مورد نظر در انتگرال نیز از کد زیر استفاده کرده‌ایم. (f_of_x(x الگوی محاسباتی را برای انتگرال‌ده (Integrand)، مشخص می‌کند.

1def f_of_x(x):
2    """
3    This is the main function we want to integrate over.
4    Args:
5    - x (float) : input to function; must be in radians
6    Return:
7    - output of function f(x) (float)
8    """
9    return (e**(-1*x))/(1+(x-1)**2)

اجرای شبیه‌سازی مونت کارلو

به منظور انجام شبیه‌سازی مونت کارلو برای محاسبه انتگرال بالا، باید مراحل زیر طی یک الگوریتم به اجرا درآیند.

1. یک  نمونه تصادفی از یک توزیع آماری (مثلا توزیع یکنواخت) در بازه انتگرال‌گیری تولید شود.
2. انتگرال‌ده براساس مقدار حاصل از مرحله ۱ محاسبه شود.
3. مراحل ۱ و ۲ به تعداد دلخواه (برای رسیدن به دقت مناسب) تکرار شوند.
4. مقدارهای حاصل از مرحله ۳ میانگین‌‌گیری شده و مقدار این میانگین در فاصله بین کران‌های انتگرال ضرب شود.

1def crude_monte_carlo(num_samples=10000):
2    """
3    This function performs the Crude Monte Carlo for our
4    specific function f(x) on the range x=0 to x=5.
5    Notice that this bound is sufficient because f(x)
6    approaches 0 at around PI.
7    Args:
8    - num_samples (float) : number of samples
9    Return:
10    - Crude Monte Carlo estimation (float)
11    
12    """
13    lower_bound = 0
14    upper_bound = 5
15    
16    sum_of_samples = 0
17    for i in range(num_samples):
18        x = get_rand_number(lower_bound, upper_bound)
19        sum_of_samples += f_of_x(x)
20    
21    return (upper_bound - lower_bound) * float(sum_of_samples/num_samples)

کدی که در بالا قابل مشاهده است، تابعی است که بر اساس 10000 نمونه تصادفی (num_samples=10000) محاسبه انتگرال رابطه ۱ را در بازه ۰ تا ۵ انجام می‌دهد. مقدار حاصل از اجرای این کد برای انتگرال برابر است با $$0.6994$$.

محاسبه واریانس برآوردگر

سوالی که در این میان مطرح می‌شود، این است که چگونه می‌توان خطای محاسباتی را مشخص کرد. آیا افزایش تعداد نمونه‌های تصادفی (مثلا رساندن تعداد آن‌ها به ۱۰۰۰۰ مشاهده) دقت محاسبات را به شکل محسوسی تغییر می‌دهد. یا به طور دقیق‌تر می‌توان این سوال را مطرح کرد که چه رابطه‌ای بین تعداد نمونه‌ها و خطای حاصل از محاسبات وجود دارد.

فیلم آموزش نرم افزار شبیه سازی مونت کارلوی GATE در مهندسی پزشکی در فرادرس

کلیک کنید

از آنجایی که مقدار واقعی انتگرال مشخص نیست نمی‌توان به راحتی فاصله محاسباتی بین مقدار واقعی و مقدار تقریبی را مبنا قرار داده و خطای محاسبه انتگرال را تعیین کرد. در عوض می‌توانیم براساس نمونه‌های مختلف، محاسبات را انجام داده و خطای میانگین (واریانس میانگین نتایج) را مبنای اندازه‌ خطا قرار دهیم. به این ترتیب واریانس یا میانگین مرابعات انحراف از میانگین را به عنوان خطا در نظر خواهیم گرفت.

$$\large \sigma^2 = E(I^2)-E^2(I)$$

همانطور که مشخص است رابطه بالا مربوط به شیوه محاسبه واریانس برای متغیر تصادفی $$I$$‌ است. از طرفی می‌دانیم که با تغییر نمونه‌ها و تکرار الگوریتم مقدار متفاوتی برای انتگرال بدست خواهد آمد. واریانس مقدارهای محاسبه شده برای انتگرال را بر مبنای دسته‌های نمونه‌های متفاوت، مبنای خطا در نظر می‌گیریم. به این ترتیب رابطه زیر را خواهیم داشت.

$$\large \sigma^2=[\dfrac{b-a}{N}\sum f^2(x_i)]-[\sum_{i=1}^N\dfrac{b-a}{N}f(x_i)]^2$$

کدی که در ادامه مشاهده می‌کنید به منظور محاسبه واریانس (خطای) شبیه‌سازی مونت‌کارلو نوشته شده است.

1def get_crude_MC_variance(num_samples):
2    """
3    This function returns the variance fo the Crude Monte Carlo.
4    Note that the inputed number of samples does not neccissarily
5    need to correspond to number of samples used in the Monte
6    Carlo Simulation.
7    Args:
8    - num_samples (int)
9    Return:
10    - Variance for Crude Monte Carlo approximation of f(x) (float)
11    """
12    int_max = 5 # this is the max of our integration range
13    
14    # get the average of squares
15    running_total = 0
16    for i in range(num_samples):
17        x = get_rand_number(0, int_max)
18        running_total += f_of_x(x)**2
19    sum_of_sqs = running_total*int_max / num_samples
20    
21    # get square of average
22    running_total = 0
23    for i in range(num_samples):
24        x = get_rand_number(0, int_max)
25        running_total = f_of_x(x)
26    sq_ave = (int_max*running_total/num_samples)**2
27    
28    return sum_of_sqs - sq_ave

براساس 10000 نقطه، مقدار واریانس براساس کد بالا برابر با $$0.266$$ است که می‌توان خطا را در این حالت به صورت زیر محاسبه کرد.

$$\large \operatorname{Error} =\left(\dfrac{\sigma^2}{N}\right)^{\tfrac{1}{2}}$$

پس مطابق با رابطه بالا میزان خطا برای محاسبه این انتگرال برابر با $$0.005$$ است. البته مقدار خطا زیاد به نظر نمی‌رسد و می‌توان به عنوان یک روش ساده و کارا، محاسبه به کمک تکنیک مونت‌کارلو را مد نظر قرار داد. برای افزایش میزان دقت برآورد باید تعداد نمونه‌ها را بیشتر کرد. ولی از طرفی افزایش تعداد نقاط، بار محاسباتی و زمان بیشتری برای اجرای الگوریتم را به همراه دارد. سوالی که در این قسمت مطرح می‌شود این است که چگونه می‌توانیم با هزینه زمانی کمتر و زمان کوتاه‌تر به پاسخ مناسب‌تر با خطای کمتری در محاسبه انتگرال به روش مونت‌کارلو رسید.

اهمیت نمونه‌گیری

یکی از روش‌های کاهش خطای در شبیه‌سازی مونت‌کارلو، استفاده از نمونه‌هایی است که از محاسبه مقدار انتگرال از اهمیت بیشتری برخوردار هستند. به این ترتیب نقاط با اهمیت در نمونه متمایز شده و از آن‌ها برای محاسبات انتگرال استفاده می‌شود. در نتیجه بدون آنکه تعداد نمونه‌ها افزایش یابد، می‌توان با دقت بیشتری عملیات شبیه‌سازی مونت‌کارلو را انجام داد.

فیلم آموزش شبیه سازی مونت کارلو و کرنر در کیدنس Cadence در فرادرس

کلیک کنید

اهمیت نمونه‌گیری را می‌توان به شکلی تاثیر نمونه‌گیری در محاسبه انتگرال در نظر گرفت. بنابراین به آن گاهی «نمونه‌گیری با اهمیت» نیز می‌گویند. انتخاب نقاطی که دارای اهمیت بیشتری هستند به توزیع داده‌ها بستگی دارد. برای مثال فرض کنید تابع پله‌ای $$f(x)$$ مطابق با تصویر زیر مورد نظر است. مشخص است که در فاصله $$[0,2]$$ این تابع مقداری برابر با ۱ داشته و در بعد از آن صفر است.