الگوریتم متروپلیس-هیستینگز

این الگوریتم به افتخار نیکلاس متروپلیس (به انگلیسی: Nicholas Metropolis) که در سال ۱۹۵۳ در مقاله‌ای با عنوان «معادلهٔ محاسبهٔ حالت به وسیلهٔ ماشین‌های رایانش سریع» به همراه همکارانش (آریانا روزنبلوث (به انگلیسی: Arianna W. Rosenbluth)، مارشال روزنبلوث (به انگلیسی: Marshall N. Rosenbluth)، اگوستا تلر (به انگلیسی: Augusta H. Teller) و ادوارد تلر (به انگلیسی: Edward Teller) منتشر کرده‌است، نام‌گذاری شده. این مقاله الگوریتمی برای توزیع‌های پیشنهادی متقارن ارائه می‌دهد و سپس هیستینگز در سال ۱۹۷۰ الگوریتم را برای حالت‌های عمومی‌تر توسعه می‌دهد.

متروپولیس و اولام در سال ۱۹۴۹ برای تولید نمونه ای از یک تابع چگالی استفاده از الگوریتم متروپولیس را توسط یک تابع چگالی پیشنهادی q متقارن معرفی نمودند. متقارن بودن به این معنی است که به ازای هر x و y داشته باشیم: (q(x,y) = q(y,x. در سال ۱۹۵۳ متروپولیس و همکارانش به مطالعه بیشتر این الگوریتم پرداختند. بعدها در سال ۱۹۷۰ هیستینگز شرط متقارن بودن را حذف نمود و الگوریتم متروپلیس-هیستینگز نامیده شد. انواع دیگر این الگوریتم بر اساس ویژگی‌های تابع تولید کنندهٔ نمونه در سال‌های بعد معرفی شدند.

در مورد این که توسعهٔ این الگوریتم دقیقاً به چه کسی نسبت داده شود بحث‌هایی وجود دارد. متروپلیس عبارت «مونت کارلو» در این زمینه را قبل‌تر در مقاله‌ای به همراه استانیسلاو اولام استفاده کرده‌بوده، همچنین با بحث‌های محاسباتی این روش آشنا بود و گروهی را در بحث‌های نظری مربوط به طراحی و ساخت کامپیوتر مانیاک۱ که برای آزمایش‌های این روش در سال ۱۹۵۲ استفاده می‌شد، رهبری می‌کرد. هرچند تا قبل از سال ۲۰۰۳ هیچ شرح دقیقی از فرایند توسعهٔ الگوریتم وجود نداشت. در آن زمان مارشال روزنبلات، کمی قبل‌تر از مرگش، در کنفرانسی در سال ۲۰۰۳ در آزمایشگاه ملی لوس‌آلموس که برای پاس‌داشت پنجاهمین سالگرد انتشار مقاله برگزار شده‌بود، شرکت کرد و الگوریتم و فرایند توسعهٔ آن را در ارائه‌ای با عنوان «پیدایش الگوریتم مونت کارلو برای مکانیک آماری» شرح داد. شفاف‌سازی‌های تاریخی بعدی به وسیلهٔ گوبرناتیس در مقاله‌ای در سال ۲۰۰۵ انجام شد. روزنبلات به صورت شفافی بیان می‌کند که او و همسرش آریانا توسعهٔ الگوریتم را انجام داده‌اند و متروپلیس نقشی در توسعهٔ آن به جز فراهم کردن زمان کامپیوتر نداشته‌است.

این مطلب با حرف‌های ادوارد تلر که در خاطراتش اذعان می‌کند پنج نویسندهٔ مقالهٔ سال ۱۹۵۳ با هم برای روزها (و شب‌ها) کار می‌کرده‌اند تناقض دارد. در مقابل، شرح با جزئیات روزنبلات بیان می‌کند که تلر نقشی اساسی ولی اولیه در پیشنهاد ایدهٔ «بهره بردن از مکانیک آماری و گرفتن مجموعه‌ای از میانگین‌ها به جای استفادهٔ مفصل از سینماتیک» داشته‌است. روزنبلات می‌گوید این موضوع باعث شد او دربارهٔ عمومی کردن راهکار مونت کارلو بیاندیشد؛ موضوعی که وی معمولاً با فون نیومن دربارهٔ آن بحث می‌کرده‌است. آریانا بیان می‌کند (به گورنتیس در سال ۲۰۰۳) که آگوستا تلر کارهای کامپیوتری را آغاز کرد ولی آریانا آن را ادامه داد و کد مربوطه را از ابتدا نوشت. روزنبلات در یک تاریخ شفاهی که چندی قبل از مرگش ضبط شده یک بار دیگر نیز تلر را کسی معرفی می‌کند که مسئلهٔ اولیه را ارائه کرده و خودش را کسی می‌داند که مسئله را حل کرده و آریانا را برنامه‌نویس راه حل این مسئله معرفی می‌کند. از دید اعتبار سخن، دلایل چندانی برای این که اعتبار حرف‌های روزنبلات را زیر سؤال ببریم وجود ندارد. در یک شرح‌حال از زندگی روزنبلات فریمن دایسون دربارهٔ او می‌نویسد:

بارهای زیادی پیش روزنبلات آمدم و از او سؤالی پرسیدم [...] و در طول دو دقیقه پاسخی از او دریافت کردم. سپس برای من معمولاً یک هفته کار می‌برد تا جزئیات این که چرا پاسخ روزنبلات درست بوده‌است را متوجه شوم. او توانایی خارق‌العاده‌ای در دیدن عمق وضعیت یک مسئلهٔ فیزیکی پیچیده و رسیدن به پاسخ درستی به کمک استدلال‌های فیزیکی داشت. انریکو فرمی تنها فیزیکدان دیگری بود که من می‌شناسم که در درک شهودی فیزیک توانایی مشابهی با روزنبلات داشت.

هنگام انجام استنتاج بیزی، هدف ما محاسبه و استفاده از توزیع مشترک پسین کامل مجموعه ای از متغیرهای تصادفی است. متأسفانه، این روش اغلب نیاز به محاسبهٔ انتگرالی غیرقابل حل دارد. در چنین مواردی ممکن است ما در برابر حل معادلات تحلیلی تسلیم شویم و روش نمونه‌گیری مبتنی بر روش زنجیره مارکوف مونت کارلو را ادامه دهیم. وقتی از زنجیره مارکوف مونت کارلو استفاده می‌کنیم با استفاده از نمونه‌ای شبیه‌سازی شده از توزیع پیشین، توزیع پسین و انتگرال‌های غیرقابل حل را تخمین می‌زنیم.

معرفی روش‌های نمونه‌برداری برای تخمین توزیع احتمال

به صورت کلی روش‌های نمونه‌برداری (که روش‌های مونت کارلو نیز نامیده می‌شوند) به چهار دسته تقسیم می‌شوند:

• نمونه‌برداری پیش‌رو (به انگلیسی: forward sampling)
• نمونه‌برداری بازپس‌زننده (به انگلیسی: rejection sampling)
• نمونه‌برداری بر پایهٔ اهمیت (به انگلیسی: importance sampling)
• زنجیرهٔ مارکف مونت کارلو (به انگلیسی: Markov chain Monte Carlo)

مشکلات روش نمونه‌برداری پیش‌رو

یکی از ساده‌ترین روشهای نمونه‌برداری در شبکه‌های بیزی روش نمونه‌برداری پیش‌رو است که با داشتن یک مرتب‌سازی موضعی از گراف، نمونه‌برداری از متغیرها را انجام می‌دهد. در این روش اگر تعداد متغیرهای شواهد زیاد باشد، بسیاری از نمونه‌ها حذف می‌شوند و در نتیجه برای تخمین دقیق‌تر باید از نمونه‌های بیشتری استفاده شود.

در بسیاری از موارد به دلیل پیچیده بودن توزیع اصلی نمی‌توان به‌صورت مستقیم از توزیع موردنظر ${\displaystyle P(x)}$

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

در روش نمونه‌برداری بازپس‌زننده اگر توزیع پیشنهادی ${\displaystyle Q(x)}$

مشکلات نمونه‌برداری بر پایهٔ اهمیت

در روش نمونه‌برداری بر پایهٔ اهمیت اگر توزیع پیشنهادی ${\displaystyle Q(x)}$

معرفی روش‌های نمونه‌برداری بر اساس زنجیرهٔ مارکف

یکی از روشهای دیگر نمونه‌برداری روش زنجیرهٔ مارکف مونت کارلو است. در این روش از توزیع پیشنهادی ${\displaystyle Q(x^{\prime })}$

روش‌های زنجیرهٔ مارکف مونت کارلو به دو روش نمونه‌برداری گیبز و نمونه‌برداری متروپلیس-هیستینگز تقسیم می‌شوند.

محدودیت‌های روش نمونه‌برداری گیبز

برای استفاده از روش گیبز باید ${\displaystyle P(x_i|x_{-i})}$

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

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

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

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

شرط این که یک زنجیرهٔ مارکف دارای توزیع حالت پایدار باشد، داشتن خاصیت برگشت‌پذیری (به انگلیسی: Reversibility) است.

خاصیت برگشت‌پذیری

اگر شرط «تعادل جزء به جزء» (به انگلیسی: Detailed Balance) برای یک زنجیرهٔ مارکف برقرار باشد، زنجیره مارکوف برگشت‌پذیر نامیده می‌شود. این شرط به صورت زیر تعریف می‌شود:

${\displaystyle \pi (x^{\prime })T(x^{\prime }\to x)=\pi (x)T(x\to x^{\prime })}$

اگر این شرط برقرار باشد، می‌توانیم بنویسیم:

${\displaystyle \sum _x\pi (x^{\prime })T(x^{\prime }\to x)=\sum _x\pi (x)T(x\to x^{\prime })⟹\pi (x^{\prime })\sum _{x}T(x^{\prime }\to x)=\sum _x\pi (x)T(x\to x^{\prime })⟹\pi (x^{\prime })=\sum _x\pi (x)T(x\to x^{\prime })}$

خاصیت منظم بودن

شرط کافی و نه لازم برای اینکه توزیع حالت پایدار زنجیرهٔ مارکف یکتا باشد؛ خاصیت منظم بودن (به انگلیسی: Regularity) آن است. زنجیرهٔ مارکف درصورت داشتن دو ویژگی زیر منظم خواهد بود:

1. کاهش‌ناپذیری (به انگلیسی: Irreducibility): احتمال رفتن از هر وضعیت به وضعیت‌های دیگر بزرگتر از ۰ باشد.
2. نامتناوب بودن (به انگلیسی: Aperiodicity): در هر لحظه از زمان امکان بازگشت به هر وضعیت وجود داشته باشد.

به عبارت دیگر در صورتی یک زنجیرهٔ مارکف دارای خاصیت منظم بودن است که حداقل یک مقدار ${\displaystyle k>0}$

در الگوریتم متروپلیس-هیستینگز به جای نمونه‌برداری از ${\displaystyle P(x)}$

در این الگوریتم هر نمونهٔ ${\displaystyle x^{\prime }}$

${\displaystyle A(x\to x^{\prime })=min(1,\frac{P(x^{\prime })Q(x|x^{\prime })}{P(x)Q(x^{\prime }|x)})}$

فرایند اجرا

در ابتدا به صورت تصادفی از یک وضعیت اولیه ${\displaystyle x_0}$

به‌همین ترتیب در یک حلقه شروع به نمونه‌گیری از توزیع ${\displaystyle Q(x^{\prime }|x)}$

ادامه دادن این روند منجر به همگرایی توزیع زنجیرهٔ استفاده شده به توزیع اصلی (که در حال تخمین آن هستیم) خواهد شد.

نمونه‌هایی که در این فرایند دور ریخته می‌شوند به عنوان نمونهٔ سوخته شناخته می‌شوند و در ادامهٔ فرایند الگوریتم به کار نمی‌آیند. مجموعهٔ باقی مانده از مقادیر مورد قبول ${\displaystyle x}$

در زیر شبه‌کدی از این الگوریتم به نمایش درآمده است.

initialize starting sate x(0)
set t=0
while samples have not converged:
    x=x(t)
    sample x* according to Q(x*|x)
    A(x to x*)=min(1,[P(x*)Q(x|x*)]/[P(x)Q(x*|x)])
    generate a uniform random number u from [0,1]
    if u < A(x to x*):
        x(t+1)=x*
    else:
        x(t+1)=x(t)
    t=t+1

اثبات همگرایی الگوریتم متروپلیس-هیستینگز به توزیع اصلی

اگر با استفاده از توزیع ${\displaystyle Q(x^{\prime }|x)}$

${\displaystyle T(x\to x^{\prime })=Q(x^{\prime }|x)A(x\to x^{\prime })}$

${\displaystyle T(x\to x)=Q(x|x)+{\mathrm{\Sigma }}_{x\ne x^{\prime }}(Q(x^{\prime }|x)(1-A(x\to x^{\prime })))}$

اگر ${\displaystyle A(x\to x^{\prime })\le 1}$

${\displaystyle A(x\to x^{\prime })=\frac{P(x^{\prime })Q(x|x^{\prime })}{P(x)Q(x^{\prime }|x)}⟹P(x)Q(x^{\prime }|x)A(x\to x^{\prime })=P(x^{\prime })Q(x|x^{\prime })A(x^{\prime }\to x)⟹P(x)T(x\to x^{\prime })=P(x^{\prime })T(x^{\prime }\to x)}$

معادلهٔ آخر همان شرط تعادل دقیق است؛ یعنی الگوریتم متروپلیس-هیستینگز باعث می‌شود به یک توزیع حالت پایدار ${\displaystyle P(x)}$

از طرف دیگر اگر ${\displaystyle T(x\to x^{\prime })>0}$

بنابراین این الگوریتم در نهایت به توزیع اصلی همگرا می‌شود.

مقدار مناسب برای واریانس توزیع گاوسی اولیه

الگوریتم متروپلیس-هیستینگز در حالتی به بهترین نحو عمل می‌کند که احتمال شرطی پیشنهاد شده، با احتمالی که در حال تخمین توزیع آن هستیم (و نمونه‌برداری مستقیم از آن دشوار است) مطابقت داشته‌باشد؛ به عبارتی باید داشته‌باشیم${\displaystyle Q(x^{\prime }|x)\approx P(x^{\prime })}$

اگر مقدار${\displaystyle {\displaystyle {\sigma }^2}}$

روش نمونه‌گیری گیبز حالت خاصی از الگوریتم متروپلیس-هیستینگز است که توزیع پیشنهادی که در هر گام برای نمونه‌برداری استفاده می‌کند تنها برای متغیر ${\displaystyle i}$

${\displaystyle Q_i(x^{\prime }|x)=P(x_i^{\prime }|x_{-i})}$

با استفاده از توزیع پیشنهادی معرفی شده در الگوریتم متروپلیس-هیستینگز می‌توان ثابت کرد که احتمال پذیرش هر نمونه در روش گیبز برابر ۱ است:

${\displaystyle A(x\to x^{\prime })=min(1,\frac{P(x^{\prime })Q(x|x^{\prime })}{P(x)Q(x^{\prime }|x)})=min(1,\frac{P(x_i^{\prime },x_{-i}^{\prime })P(x_i|x_{-i}^{\prime })}{P(x_i,x_{-i})P(x_i^{\prime }|x_{-i})})x_{-i}^{\prime }=x_{-i}⟹min(1,\frac{P(x_i^{\prime }|x_{-i})P(x_{-i})P(x_i|x_{-i})}{P(x_i|x_{-i})P(x_{-i})P(x_i^{\prime }|x_{-i})})=1}$

یکی از کاربردهای مرسوم الگوریتم متروپلیس-هیستینگز در محاسبه‌ي یک انتگرال است. فضای ${\displaystyle \mathrm{\Omega }\subset \mathbb{R}}$

${\displaystyle P(E)={\int }_{\mathrm{\Omega }}A(x)P(x)dx}$

که در آن ${\displaystyle A(x)}$

${\displaystyle P(E)={\int }_{\mathrm{\Omega }}P(E|x)P(x)dx={\int }_{\mathrm{\Omega }}\delta (E-E(x))P(x)dx=E_X(P(E|X))}$

و بنابراین تخمین ${\displaystyle P(E)}$