منبع اصلی نوشتار زیر در این لینک قرار دارد

مقدمه‌ای بر شبیه‌سازی در R، روش مونت‌کارلو

یکی از کاربرد‌های کمتر شناخته‌شده‌ی R، استفاده از آن برای شبیه‌سازی است. شبیه‌سازی سیستم‌های گسسته‌ پیش‌آمد(DES) ، سیستم‌های صف و روش‌های مانند مونت‌کارلو رایج‌ترین شبیه‌سازی‌هایی هستند که در R انجام می‌شوند. البته این را هم در نظر داشته باشید که R کمی تا قسمتی کند است و چیزی که آن را برای شبیه‌سازی مطلوب کرده این است که بسیاری از محققینی که زمینه‌ی آماری دارند R را به زبان‌ها و ابزارهای دیگر ترجیح می‌دهند، از طرف دیگر، متن‌باز بودن R  مزیت آن در این زمینه در مقابل Matlab است. از آنجایی که این وب‌سایت درباره‌ی R و مخلفاتش است، پرداختن به مبحث شبیه‌سازی در آن خالی از لطف نیست. تقریبا هر جایی که مبحث شبیه‌سازی مطرح می‌شود، روش مونت‌کارلو اولین روشی است که مورد بحث قرار می‌گیرد.

مونت‌کارلو را دانشمندانی که روی پروژه منهتن کار می‌کردند در دهه ۴۰ میلادی ابداع نمودند و از روی شهری به همین نام در موناکو نام‌گزاری شده است. ایده‌ی اصلی روش مونت‌کارلو، استفاده از نمونه‌هایی تصادفی از پارامترها یا ورودی‌ها، برای تحلیل رفتار یک سیستم یا فرایند یا حل یک مسئله است. بهتر است با ساده‌ترین مثالی که برای بیان روش مونت‌کارلو وجود دارد به تشریح این روش بپردازیم، تخمین عدد پی.

pi-1338548_960_720

می‌دانیم که مساحت یک دایره، برابر \pi r^{2} است، پس اگر شعاع و مساحت یک دایره را بدانیم می‌توانیم عدد پی را بدست آوریم. فرض کنید مربعی با ضلع دو را به یک دایره محیط کرده‌ایم، بدیهی است که شعاع دایره برابر یک خواهد بود.

cr.

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

mcpi <- function(n) {
x = runif(n,min = 0,max = 2);
y = runif(n,min = 0,max = 2);
i <- ((x-1)^2 + (y-1)^2 <= 1);
pi <- (sum(i)/n)*4;
plot(x,y,pch='.',col=ifelse(i,"red","grey"),xlab='',ylab='',asp=1,main=paste(" Pi =",pi));
}

تابع runif داده‌های تصادفی بر اساس توزیع یکنواخت به وجود می‌آورد، تابع plot برای رسم نمودار است، درباره‌ی این توابع در آینده‌ی نزدیک بیشتر خواهم نوشت. این کد را در پنجره source می‌نویسیم و اجرا می‌کنیم تا تابع mcpi در حافظه به وجود بیاید. تابعی که نوشتیم عدد n را به عنوان ورودی می‌گیرد، به تعداد n، داخل مربع نقاط تصادفی به وجود می‌آورد، سپس تعداد آنها را می‌شمارد و در قالب یک نمودار مختصات نمایش می‌دهد. این تابع را با اعداد مختلف امتحان می‌کنیم:

۱۰۰۰ نقطه: عدد ۳.۱۲ بدست آمده که خیلی دقیق نیست.

Rplot3

۱۰۰۰۰۰نقطه: با صد برابر نقطه نسبت به مرحله قبلی، تخمین ما فقط کمی به عدد اصلی نزدیکتر شده است.

Rplot2

۱۰۰۰۰۰۰نقطه: با یک میلیون نقطه، این روش تقریبا تا ۳ رقم اعشار دقیق است.

Rplot1

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