الگوریتم متروپلیس-هیستینگز
در علم آمار و فیزیک آماری، الگوریتم متروپلیس-هیستینگز یک روش زنجیره مارکوف مونت کارلو برای بدست آوردن ترتیبی از نمونه‌های تصادفی از یک توزیع احتمالی است که نمونه‌برداری مستقیم از آن دشوار می‌باشد. این ترتیب را می‌توان برای برآورد یک توزیع (به عنوان مثال تولید یک هیستوگرام) یا برای محاسبهٔ برخی انتگرال‌ها (به‌طور مثال یک امید ریاضی) استفاده کرد. متروپلیس-هیستینگز و دیگر الگوریتم‌های MCMC عموماً برای نمونه‌برداری از توزیع‌های چند بعدی استفاده می‌شوند، به خصوص زمانی که تعداد ابعاد بالا باشد. برای توزیع‌های تک بعدی معمولاً روش‌های دیگری وجود دارند (به عنوان مثال نمونه‌گیری عدم پذیرش تطبیقی) که می‌توانند به‌طور مستقیم نمونه‌های مستقل را از توزیع بازگردانند؛ معمولاً این روش‌ها با مشکل خودهمبستگی که در روش‌های MCMC به صورت ذاتی وجود دارند، رو به رو نیستند.
این الگوریتم به افتخار نیکلاس متروپلیس که در سال ۱۹۵۳ در مقاله‌ای با عنوان «معادلهٔ محاسبهٔ حالت به وسیلهٔ ماشین‌های رایانش سریع» به همراه همکارانش (آریانا روزنبلوث ، مارشال روزنبلوث ، اگوستا تلر و ادوارد تلر منتشر کرده‌است، نام‌گذاری شده. این مقاله الگوریتمی برای توزیع‌های پیشنهادی متقارن ارائه می‌دهد و سپس هیستینگز در سال ۱۹۷۰ الگوریتم را برای حالت‌های عمومی‌تر توسعه می‌دهد.
متروپولیس و اولام در سال ۱۹۴۹ برای تولید نمونه ای از یک تابع چگالی استفاده از الگوریتم متروپولیس را توسط یک تابع چگالی پیشنهادی "q" متقارن معرفی نمودند. متقارن بودن به این معنی است که به ازای هر "x" و "y" داشته باشیم: "(q(x,y) = q(y,x". در سال ۱۹۵۳ متروپولیس و همکارانش به مطالعه بیشتر این الگوریتم پرداختند. بعدها در سال ۱۹۷۰ هیستینگز شرط متقارن بودن را حذف نمود و الگوریتم متروپلیس-هیستینگز نامیده شد. انواع دیگر این الگوریتم بر اساس ویژگی‌های تابع تولید کنندهٔ نمونه در سال‌های بعد معرفی شدند.
در مورد این که توسعهٔ این الگوریتم دقیقاً به چه کسی نسبت داده شود بحث‌هایی وجود دارد. متروپلیس عبارت «مونت کارلو» در این زمینه را قبل‌تر در مقاله‌ای به همراه استانیسلاو اولام استفاده کرده‌بوده، همچنین با بحث‌های محاسباتی این روش آشنا بود و گروهی را در بحث‌های نظری مربوط به طراحی و ساخت کامپیوتر مانیاک۱ که برای آزمایش‌های این روش در سال ۱۹۵۲ استفاده می‌شد، رهبری می‌کرد. هرچند تا قبل از سال ۲۰۰۳ هیچ شرح دقیقی از فرایند توسعهٔ الگوریتم وجود نداشت. در آن زمان مارشال روزنبلات، کمی قبل‌تر از مرگش، در کنفرانسی در سال ۲۰۰۳ در آزمایشگاه ملی لوس‌آلموس که برای پاس‌داشت پنجاهمین سالگرد انتشار مقاله برگزار شده‌بود، شرکت کرد و الگوریتم و فرایند توسعهٔ آن را در ارائه‌ای با عنوان «پیدایش الگوریتم مونت کارلو برای مکانیک آماری» شرح داد. شفاف‌سازی‌های تاریخی بعدی به وسیلهٔ گوبرناتیس در مقاله‌ای در سال ۲۰۰۵ انجام شد. روزنبلات به صورت شفافی بیان می‌کند که او و همسرش آریانا توسعهٔ الگوریتم را انجام داده‌اند و متروپلیس نقشی در توسعهٔ آن به جز فراهم کردن زمان کامپیوتر نداشته‌است.
این مطلب با حرف‌های ادوارد تلر که در خاطراتش اذعان می‌کند پنج نویسندهٔ مقالهٔ سال ۱۹۵۳ با هم برای روزها (و شب‌ها) کار می‌کرده‌اند تناقض دارد. در مقابل، شرح با جزئیات روزنبلات بیان می‌کند که تلر نقشی اساسی ولی اولیه در پیشنهاد ایدهٔ «بهره بردن از مکانیک آماری و گرفتن مجموعه‌ای از میانگین‌ها به جای استفادهٔ مفصل از سینماتیک» داشته‌است. روزنبلات می‌گوید این موضوع باعث شد او دربارهٔ عمومی کردن راهکار مونت کارلو بیاندیشد؛ موضوعی که وی معمولاً با فون نیومن دربارهٔ آن بحث می‌کرده‌است. آریانا بیان می‌کند (به گورنتیس در سال ۲۰۰۳) که آگوستا تلر کارهای کامپیوتری را آغاز کرد ولی آریانا آن را ادامه داد و کد مربوطه را از ابتدا نوشت. روزنبلات در یک تاریخ شفاهی که چندی قبل از مرگش ضبط شده یک بار دیگر نیز تلر را کسی معرفی می‌کند که مسئلهٔ اولیه را ارائه کرده و خودش را کسی می‌داند که مسئله را حل کرده و آریانا را برنامه‌نویس راه حل این مسئله معرفی می‌کند. از دید اعتبار سخن، دلایل چندانی برای این که اعتبار حرف‌های روزنبلات را زیر سؤال ببریم وجود ندارد. در یک شرح‌حال از زندگی روزنبلات فریمن دایسون دربارهٔ او می‌نویسد:بارهای زیادی پیش روزنبلات آمدم و از او سؤالی پرسیدم [...] و در طول دو دقیقه پاسخی از او دریافت کردم. سپس برای من معمولاً یک هفته کار می‌برد تا جزئیات این که چرا پاسخ روزنبلات درست بوده‌است را متوجه شوم. او توانایی خارق‌العاده‌ای در دیدن عمق وضعیت یک مسئلهٔ فیزیکی پیچیده و رسیدن به پاسخ درستی به کمک استدلال‌های فیزیکی داشت. انریکو فرمی تنها فیزیکدان دیگری بود که من می‌شناسم که در درک شهودی فیزیک توانایی مشابهی با روزنبلات داشت.
هنگام انجام استنتاج بیزی، هدف ما محاسبه و استفاده از توزیع مشترک پسین کامل مجموعه ای از متغیرهای تصادفی است. متأسفانه، این روش اغلب نیاز به محاسبهٔ انتگرالی غیرقابل حل دارد. در چنین مواردی ممکن است ما در برابر حل معادلات تحلیلی تسلیم شویم و روش نمونه‌گیری مبتنی بر روش زنجیره مارکوف مونت کارلو را ادامه دهیم. وقتی از زنجیره مارکوف مونت کارلو استفاده می‌کنیم با استفاده از نمونه‌ای شبیه‌سازی شده از توزیع پیشین، توزیع پسین و انتگرال‌های غیرقابل حل را تخمین می‌زنیم.
به صورت کلی روش‌های نمونه‌برداری (که روش‌های مونت کارلو نیز نامیده می‌شوند) به چهار دسته تقسیم می‌شوند:
یکی از ساده‌ترین روشهای نمونه‌برداری در شبکه‌های بیزی روش نمونه‌برداری پیش‌رو است که با داشتن یک مرتب‌سازی موضعی از گراف، نمونه‌برداری از متغیرها را انجام می‌دهد. در این روش اگر تعداد متغیرهای شواهد زیاد باشد، بسیاری از نمونه‌ها حذف می‌شوند و در نتیجه برای تخمین دقیق‌تر باید از نمونه‌های بیشتری استفاده شود.
در بسیاری از موارد به دلیل پیچیده بودن توزیع اصلی نمی‌توان به‌صورت مستقیم از توزیع موردنظر formula_1نمونه‌برداری کرد، بلکه از توزیع پیشنهادی formula_2نمونه‌برداری می‌کنیم. روشهای نمونه‌برداری بازپس‌زننده و نمونه‌برداری بر پایهٔ اهمیت از این ایده استفاده می‌کنند.
در روش نمونه‌برداری بازپس‌زننده اگر توزیع پیشنهادی formula_2 تفاوت زیادی با formula_1 داشته باشد؛ اکثر نمونه‌ها پذیرفته نمی‌شوند. همچنین در این روش‌ها مشکل نفرین ابعاد وجود دارد؛ یعنی با افزایش ابعاد توزیع موردنظر، احتمال پذیرفته نشدن نمونه‌ها به صورت نمایی افزایش می‌یابد.
در روش نمونه‌برداری بر پایهٔ اهمیت اگر توزیع پیشنهادی formula_2تفاوت زیادی با formula_1داشته باشد؛ اکثر نمونه‌ها وزن بسیار کمی خواهند داشت و با توجه به وزن آن‌ها رسیدن به نمونه‌های مناسبی از توزیع دشوار خواهد بود.
یکی از روشهای دیگر نمونه‌برداری روش زنجیرهٔ مارکف مونت کارلو است. در این روش از توزیع پیشنهادی formula_7استفاده نمی‌شود، بلکه از توزیع پیشنهادی formula_8استفاده می‌شود که در آن formula_9نمونه جدید به دست آمده و formula_10آخرین نمونهٔ پذیرفته شده‌است. همزمان با تغییر formula_10، formula_8نیز تغییر می‌کند.
روش‌های زنجیرهٔ مارکف مونت کارلو به دو روش نمونه‌برداری گیبز و نمونه‌برداری متروپلیس-هیستینگز تقسیم می‌شوند.
برای استفاده از روش گیبز باید formula_13را محاسبه کنیم. اگرچه این توزیع برای مدل‌های گرافی گسسته قابل بدست آوردن است؛ اما برای مدل‌های گرافی پیوسته به راحتی قابل محاسبه نیست. محدودیت مهمتری که در این روش وجود دارد این است که زنجیرهٔ گیبز تنها حرکات محدودی در فضای دامنهٔ متغییرهای تصادفی توزیع انجام می‌دهد. در این روش هر بار تنها یک متغیر تغییر می‌کند. در مدل‌هایی که متغییرهای آن به شدت به هم وابسته هستند، چنین تغییری منجر به این خواهد شد که از حالت‌هایی که احتمال بالایی برای وقوع دارند به حالت‌هایی منتقل شویم که احتمال بسیار پایینی برای وقوع دارند. در این صورت حالت‌های با احتمال بالا بیشتر نمونه‌های تولید شده را دربر خواهندگرفت و این موضوع باعث می‌شود زنجیره به سختی بتواند از آن نمونه‌ها فاصله بگیرد که به معنی سرعت پایین همگرایی آن به توزیع اصلی خواهد بود. در این صورت ما معمولاً روش‌هایی را ترجیح می‌دهیم که حرکات آن‌ها برد بیشتری داشته باشد و در نتیجه اجازهٔ حرکت بیشتری در فضای دامنهٔ توزیع احتمالی به ما بدهد. الگوریتم متروپلیس-هیستینگز یکی از این راه‌کارها می‌باشد.
به صورت کلی محدودیت‌های بسیاری برای روش‌های بیزی وجود دارد. حتی با داشتن توزیع مشترک پسین، ممکن است استفاده از آن برای به دست آوردن توزیع شرطی برای هر یک از متغییرهای تصادفی موجود در مدل امکان‌پذیر و عملی نباشد. با داشتن احتمال شرطی پسین همچنان ممکن است آن‌ها از یک فرم شناخته شده نباشند و راه ساده‌ای برای نمونه‌گیری از آن‌ها وجود نداشته باشد.
روش‌های زنجیرهٔ مارکف مونت کارلو (MCMC) با قدم زدن تصادفی در طول یک زنجیرهٔ مارکف نمونه‌برداری از یک توزیع را انجام می‌دهند. به‌عبارت دیگر توزیع اصلی موردنظر را برابر توزیع حالت پایدار یک زنجیرهٔ مارکف در نظر می‌گیرند و بنابراین اگر بتوان زنجیرهٔ مارکفی طراحی کرد که دارای توزیع حالت پایدار باشد، می‌توان با شروع از یک وضعیت اولیهٔ تصادفی در زنجیره مارکف و حرکت مطابق با ماتریس انتقال بین وضعیت‌ها، درنهایت پس از تعدادی حرکت در زنجیرهٔ مارکف به نمونه‌ای از توزیع حالت پایدار و در نتیجه نمونه‌هایی از توزیع اصلی رسید.
باید توجه داشت که یک زنجیرهٔ مارکف برای این که مناسب نمونه‌برداری باشد، باید هم توزیع حالت پایدار داشته‌باشد و هم این توزیع برای آن صرف‌نظر از وضعیت شروع، یکتا باشد. این موضوع به عنوان یک قضیه در نظریهٔ مدل‌های گرافی احتمالی مطرح و اثبات می‌شود.
شرط این که یک زنجیرهٔ مارکف دارای توزیع حالت پایدار باشد، داشتن خاصیت برگشت‌پذیری است.
اگر شرط «تعادل جزء به جزء» برای یک زنجیرهٔ مارکف برقرار باشد، زنجیره مارکوف برگشت‌پذیر نامیده می‌شود. این شرط به صورت زیر تعریف می‌شود:
formula_14
اگر این شرط برقرار باشد، می‌توانیم بنویسیم:
formula_15 رابطهٔ اخیر رابطهٔ توزیع حالت پایدار است.
شرط کافی و نه لازم برای اینکه توزیع حالت پایدار زنجیرهٔ مارکف یکتا باشد؛ خاصیت منظم بودن آن است. زنجیرهٔ مارکف درصورت داشتن دو ویژگی زیر منظم خواهد بود:به عبارت دیگر در صورتی یک زنجیرهٔ مارکف دارای خاصیت منظم بودن است که حداقل یک مقدار formula_16وجود داشته باشد که به ازای هر دو وضعیت formula_10و formula_9بتوان با دقیقاً formula_19گام از وضعیت formula_10به formula_9 رفت.
در الگوریتم متروپلیس-هیستینگز به جای نمونه‌برداری از formula_1از توزیع شرطی formula_8نمونه‌برداری می‌کنیم و معمولاً formula_8به صورت توزیع نرمال formula_25در نظر گرفته می‌شود.
در این الگوریتم هر نمونهٔ formula_9که بعد از نمونهٔ formula_10تولید شده، احتمال مشخصی برای پذیرفته شدن دارد. این احتمال با formula_28نمایش داده شده و به صورت زیر محاسبه می‌شود:
formula_29
در ابتدا به صورت تصادفی از یک وضعیت اولیه formula_30شروع می‌کنیم و از توزیع formula_31نمونه‌برداری می‌کنیم و نمونهٔ formula_32را تولید می‌کنیم. این نمونه با احتمال formula_33پذیرفته می‌شود.به‌همین ترتیب در یک حلقه شروع به نمونه‌گیری از توزیع formula_8و تولید نمونهٔ formula_9می‌کنیم. هر نمونهٔ formula_9با احتمال formula_28شانس پذیرفته شدن دارد. دقت داریم که تابع A شامل یک کمینه‌گیری است و اگر مقدار پارامتر دیگر کمینه‌گیری بزرگ‌تر از یک شود، تابع مقدار یک را برخواهد گرداند و نمونه حتماً پذیرفته می‌شود. در صورتی که مقدار پارامتر دوم کمتر از یک شود، عددی تصادفی در بازهٔ [۰٬۱] تولید می‌کنیم و اگر کوچک‌تر یا برابر با مقدار پارامتر دوم شد، نمونه را می‌پذیریم و در غیر این صورت نمونهٔ جدید را برابر با نمونهٔ تولید شده در مرحلهٔ قبل قرار می‌دهیم. ادامه دادن این روند منجر به همگرایی توزیع زنجیرهٔ استفاده شده به توزیع اصلی (که در حال تخمین آن هستیم) خواهد شد.
نمونه‌هایی که در این فرایند دور ریخته می‌شوند به عنوان نمونهٔ سوخته شناخته می‌شوند و در ادامهٔ فرایند الگوریتم به کار نمی‌آیند. مجموعهٔ باقی مانده از مقادیر مورد قبول formula_10نمونه‌هایی از توزیع formula_1را نمایش می‌دهند.
در زیر شبه‌کدی از این الگوریتم به نمایش درآمده است. initialize starting sate x(0) set t=0 while samples have not converged: اگر با استفاده از توزیع formula_8نمونه‌یformula_9را به دست آوریم، آن‌گاه احتمال رفتن از وضعیت formula_10به formula_9برابر است با:
formula_44 (برای حالتی که formula_10و formula_9برابر نباشند)
formula_47 (برای حالتی که formula_10و formula_9برابر باشند)
اگر formula_50آن‌گاه formula_51و در نتیجه formula_52خواهد بود. در این صورت می‌توانیم بنویسیم:
formula_53
معادلهٔ آخر همان شرط تعادل دقیق است؛ یعنی الگوریتم متروپلیس-هیستینگز باعث می‌شود به یک توزیع حالت پایدار formula_1در زنجیرهٔ مارکف برسیم.
از طرف دیگر اگر formula_55باشد؛ آن‌گاه شرط منظم بودن زنجیرهٔ مارکف نیز برقرار است.
بنابراین این الگوریتم در نهایت به توزیع اصلی همگرا می‌شود.
الگوریتم متروپلیس-هیستینگز در حالتی به بهترین نحو عمل می‌کند که احتمال شرطی پیشنهاد شده، با احتمالی که در حال تخمین توزیع آن هستیم (و نمونه‌برداری مستقیم از آن دشوار است) مطابقت داشته‌باشد؛ به عبارتی باید داشته‌باشیمformula_56. اگر از یک توزیع احتمالی گاوسی برای این کار استفاده شود، باید پارامتر واریانس آن formula_57در طول فرایند انتخاب نمونه‌های مناسب تنظیم شود. معمولاً این کار با محاسبهٔ یک نرخ پذیرش صورت می‌گیرد. این نرخ نیز براساس تعداد نمونه‌هایی که در طول تعداد دفعات مشخصی از نمونه‌گیری در گذشته مورد قبول واقع شده‌اند؛ حساب می‌شود. مقدار نرخ پذیرش مطلوب بستگی به توزیع هدف دارد، هر چند از دید نظری نشان داده شده که نرخ پذیرش ایده‌آل برای یک توزیع گاوسی تک بعدی حدود ۵۰٪ است. این مقدار هنگامی که از یک توزیع N بعدی گاوسی استفاده شود به ۲۳٪ کاهش پیدا می‌کند.
اگر مقدارformula_58خیلی کوچک انتخاب شود، زنجیره با سرعت پایینی به تعادل خواهد رسید (به عبارت دیگر نرخ پذیرش بالا خواهد بود ولی نمونه‌های پیاپی گرفته شده، با سرعت پایینی در فضای دامنهٔ متغیر توزیع پخش خواهند شد و بنابراین زنجیره با سرعت پایینی به formula_1همگرا می‌شود). از طرفی اگرformula_58خیلی بزرگ انتخاب شود، نرخ پذیرش بسیار پایین خواهد بود؛ چرا که نمونه‌های پیشنهادی احتمالاً در نقاطی از دامنهٔ چگالیِ احتمال خواهند افتاد که مقدار احتمال بسیار پایینی دارد و بنابراین formula_61بسیار کوچک خواهد شد و دوباره سرعت همگرایی زنجیره بسیار پایین خواهد بود. با توجه به توضیحاتی که در مورد تخمین‌های نظری در این زمینه ارائه شد، این که توزیع پیشنهادی را طوری تنظیم کنیم که حدود ۳۰٪ از نمونه‌ها مورد پذیرش واقع شوند، روش مرسومی شناخته می‌شود.
روش نمونه‌گیری گیبز حالت خاصی از الگوریتم متروپلیس-هیستینگز است که توزیع پیشنهادی که در هر گام برای نمونه‌برداری استفاده می‌کند تنها برای متغیر formula_62م (که جزئی از متغییرهای تصادفی توزیع اصلی است) تعریف می‌شود و برابر است با:
formula_63
با استفاده از توزیع پیشنهادی معرفی شده در الگوریتم متروپلیس-هیستینگز می‌توان ثابت کرد که احتمال پذیرش هر نمونه در روش گیبز برابر ۱ است:
formula_64
یکی از کاربردهای مرسوم الگوریتم متروپلیس-هیستینگز در محاسبه‌ي یک انتگرال است. فضای formula_65 و یک توزیع احتمالی formula_1روی formula_67 و formula_68 را در نظر بگیرید. الگوریتم متروپلیس-هیستینگز می‌تواند انتگرالی در فرم زیر را تخمین بزند:
formula_69
که در آن formula_70تابعی دلخواه است. برای مثال آماره‌ي formula_71 و توزیع احتمالی آن formula_72 را که یک توزیع احتمالی حاشیه‌ای است، در نظر بگیرید. فرض کنیم هدف تخمین formula_72 برای formula_74براساس دنباله‌ای از formula_72باشد. formula_72 را می‌توانیم به صورت زیر بازنویسی کنیم:
formula_77
و بنابراین تخمین formula_72 می‌تواند با تخمین امید ریاضی مقدار تابع مشخصه‌یformula_79که مقدار آن تنها هنگامی یک خواهد بود کهformula_80(و درغیراینصورت صفر خواهدبود)، صورت پذیرد. به این دلیل که formula_74 روی دنباله‌ای از formula_72 محاسبه می‌شود؛ احتمال این که حالت formula_10با formula_71روی دنباله‌ای از formula_72به دست بیاید، متناسب با formula_72 است، که طبق تعریف مقدار کوچکی خواهد بود. الگوریتم متروپلیس-هیستینگ می‌تواند در این‌جا برای نمونه‌گیری از حالت‌هایی که نادر هستند استفاده شود و آن‌ها را با احتمال بالاتری تولید کند که منجر به افزایش تعداد نمونه‌های مورد استفاده برای تخمین formula_72خواهدشد. این روش می‌تواند به عنوان مثال با استفاده از یک توزیع مولد نمونه مثل formula_88که بیشتر به آن حالت‌ها شانس تولید می‌دهد، پیاده‌سازی شود (برای مثال formula_88 می‌تواند به شکل formula_90 با formula_91 تعریف شود).