فیلتر ذره‌ای

از testwiki
پرش به ناوبری پرش به جستجو

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

برآورد برپایه نمونه‌برداری Monte-Carlo

فرض کنید یک آزمایش تصادفی را به صورت مستقل از هم N بار انجام دهیم و پیشامد S را در طول این آزمایش مدنظر قرار دهیم. پس در این آزمایش تعداد دفعاتی که پیشامد S رخ می‌دهد را شمرده و باNs نشان می‌دهیم. در نتیجه احتمال رخداد پیشامد S به صورت زیر بدست خواهد آمد.الگو:سخ الگو:چپ‌چین Pr(S)=limNNsN الگو:سخ الگو:پایان چپ‌چین این رابطه بیان می‌کند که چنانچه آزمایش را با تعداد دفعات نامحدود و به صورت کاملاً مستقل از هم انجام دهیم، آنگاه می‌توان احتمال رخداد پیشامدی مانند S را با شمردن تعداد رویداد آن پیشامد و تقسیم کردن بر تعداد دفعات کل بدست آورد. اما در جهان واقع به دلیل اینکه تعداد دفعات آزمایش محدود است و وابستگی‌هایی بین آزمایش‌ها وجود دارد می‌توان احتمال رخداد پیشامدS را به‌صورت زیر تقریب زد.الگو:سخ الگو:چپ‌چین Pr^(S)=NsN الگو:سخ الگو:پایان چپ‌چین به این تقریب، تقریب مونت کارلو (Monte-Carlo) گفته می‌شود[۱]. اینک با توجه به اینکه تعداد رویدادهای پیشامد مورد نظر برای N‌های محدود یک متغیر تصادفی است، تقریب احتمال رخداد پیشامد S (رابطه بالا) نیز یک متغیر تصادفی خواهد بود. اما برای آنکه برآورد به روش مونت‌کارلو کیفیت مناسبی داشته باشد باید دارای ویژگی‌های زیر باشد:

ویژگی اول: برآورد برپایه مونت‌کارلو باید نااریب باشد، یعنی اگر با انجام آزمایش‌های متعدد مقدارهای برآورد متفاوتی بدست‌آوردیم باید میانگین این مقادیر دقیقاً برابر با مقدار مورد برآورد باشد؛ یعنی اگر به دفعات با به‌کارگیری از روش مونت کارلو، احتمال رخداد پیشامد S را برآورد کردیم و مجموعه‌ای از برآوردها (تخمین‌ها) بدست آوردیم، میانگین این برآوردها باید برابر با مقدار واقعی احتمال رخداد باشد.

ویژگی دوم: برآورد برپایه مونت‌کارلو باید سازگار (Consistent) باشد، یعنی واریانس مقادیر مختلف از مجموعه برآوردها بر پایه روش مونت‌کارلو با افزایش تعداد دفعات آزمایش به سمت صفر بگراید.

در نتیجه خطای برآورد نااریب و سازگار دارای میانگین آماری صفر بوده و واریانس آن با افزایش تعداد تکرر آزمایش به سمت صفر می‌گراید. به عنوان یک مثال کاربردی از روش برآورد مونت‌کارلو می‌توان به برآورد مساحت اشاره‌کرد. فرض کنید برآورد مساحت یک ناحیه برای ما مطلوب است و این ناحیه مورد برآورد در یک ناحیه مشخص دیگری پَروسته (محصور) شده باشد که مساحت آن دانسته‌است؛ بنابراین بر طبق روش مونت‌کارلو یک آزمایش تصادفی که جزئیات آن بیان می‌شود را انجام می‌دهیم؛ ریختن ذره‌هایی به صورت کاملاً مستقل و تصادفی همگی در درون ناحیه‌ای که مساحت آن را می‌دانیم؛ بنابراین آزمایش انجام شد. اینک پیشامد برآورد شده S را به‌صورت تعداد ذره‌های که در درون ناحیه‌ای که می‌خواستیم تخمین بزنیم می‌افتند تعریف می‌کنیم. اکنون تعداد ذره‌هایی را که در درون ناحیه با مساحت نامعلوم می‌افتد با Ns و تعداد ذره‌هایی را که در درون ناحیه با مساحت معلوم می‌افتد با N نشان می‌دهیم. شکل زیر یک مثال از ناحیه‌ای که مساحتش را می‌خواهیم برآورد کنیم و ناحیه با مساحت معلوم را نشان می‌دهد.

تقریب مساحت نامعلوم As (ناحیه ربع دایره) محصور (پَروسته) در مساحت معلوم مستطیلی شکل به کمک برآورد مونت کارلو

اگر مساحت ناحیه معلوم را با A و مساحت ناحیه برآورد شده را با As نشان دهیم، آنگاه نسبت مساحت ناحیه نامعلوم به کل مساحت ناحیه معلوم بر طبق این آزمایش تقریباً برابر است با:الگو:سخ الگو:چپ‌چین AsANsN الگو:پایان چپ‌چین الگو:سخ اینک با توجه به اینکه مقدار A معلوم است، مساحت ناحیه برآورد شده به صورت تقریبی از رابطه زیر بدست می‌آید. الگو:سخ الگو:چپ‌چین AsANsN الگو:پایان چپ‌چین الگو:سخ برای اینکه برآورد به صورت نااریب باشد لازم است ذره‌هایی که در درون ناحیه با مساحت معلوم قرار می‌گیرند دارای توزیع یکنواخت روی آن ناحیه باشند.

اما روش مونت‌کارلو به اینجا محدود نمی‌شود، فرض کنید ذره‌های {x(i)}i=1N به یک طریقی و به‌طور مستقل از هم و تصادفی از توزیع احتمال P(x) بیرون کشیده شده باشند. اینک برپایه روش مونت‌کارلو تابع توزیع P(x) به کمک ذره‌های بیرون کشیده شده به صورت زیر تقریب زده می‌شود. الگو:سخ الگو:چپ‌چین P(x)1Ni=1Nδ(xx(i)) الگو:پایان چپ‌چین الگو:سخ که در آن δ(.) تابع دلتای دیراک و N تعداد ذره‌ها است. در نتیجه مقدار متوسط هر تابع دلخواه f(x) را می‌توان به صورت زیر تقریب زد. الگو:سخ الگو:چپ‌چین الگو:NumBlk الگو:پایان چپ‌چین الگو:سخ رابطه اخیر بیان می‌دارد که برای محاسبه امید ریاضی هر تابع دلخواه f(x) کافیست N ذره از تابع توزیع P(x) نمونه داشت و مقدار آنها تحت تابعیت f(x) را محاسبه کرد. در صورتی‌که 𝔼(f(x)) وجود داشته باشد بر طبق قانون ضعیف اعداد بزرگ برای هر کمیت کوچک اپسیلون رابطه زیر را می‌توان نوشت[۲]. الگو:سخ الگو:چپ‌چین limN(|f^n𝔼(f(x))|ϵ)=0 الگو:پایان چپ‌چین الگو:سخ بر طبق این رابطه هرچه N بزرگتر شود احتمال اینکه برآورد متوسط تابع f(x) از مقدر واقعی آن یعنی 𝔼(f(x)) دور شود کمتر است. از طرف دیگر این برآورد یک برآورد نااریب است چراکه داریم:

𝔼(f^N)=𝔼(1Ni=1Nf(x(i)))=1Ni=1N𝔼(f(x(i)))=𝔼(f(x))

همچنین این برآورد یک برآورد سازگار است، برای فهم آن نیز داریم:

الگو:چپ‌چین Var(f^N)=𝔼(f^N2)𝔼(f(x))2=𝔼(1N2i=1Nj=1Nf(x(i))f(x(j)))𝔼(f(x))2=1N2((N2N)𝔼(f(x))2+N𝔼(f(x)2))𝔼(f(x))2=1N(𝔼(f(x)2)𝔼(f(x))2)=1NVar(f(x)) الگو:پایان چپ‌چین الگو:سخ که در بدست آوردن رابطه اخیر، معادله زیر و با توجه به i.i.d. بودن ذره‌های بیرون کشیده شده بکار برده شد.

𝔼(f(x(i))f(x(j)))={𝔼(f(x)2)i=j,𝔼(f(x))2ij

اینک با توجه به اینکه بسیاری از کمیت‌ها و پارامترهایی را که می‌خواهیم برآورد کنیم می‌توان به صورت امید ریاضی نوشت و امید ریاضی را نیز می‌توان به کمک روش بیان شده مونت‌کارلو برآورد کرد، گستره کاربرد روش مونت‌کارلو روشن می‌شود. برای مثال فرض کنید بخواهیم کمیت x را برپایه تخمین MMSE و با در نظر گرفتن کمیت مشاهده شدهy تخمین بزنیم. بر طبق مرجع[۳] می‌دانیم که این برآورد به رابطه زیر می‌انجامد.

x^MMSE=𝔼(x|y)

بنابراین باید تابع f(x)=x را بر روی توزیع P(x|y) متوسط بگیریم. بدین منظور کافیست ذره‌های مستقل و تصادفی {x(i)}i=1N را از توزیع P(x|y) بیرون کشیده و طبق برآورد مونت‌کارلو در رابطه زیر قرار دهیم.

x^MMSE1Ni=1Nx(i)

در ادامه پیرامون نمونه‌برداری بااهمیت به عنوان شاخه‌ای از برآورد مونت‌کارلو بحث می‌کنیم.

فریافت نمونه‌برداری بااهمیت

نمونه‌برداری بااهمیت (Importance Sampling) به عنوان شاخه‌ای از تخمین برپایه مونت‌کارلو اولین بار توسط مارشال[۴] معرفی شد. ایده اصلی آن بر پایه نمونه‌برداری از توزیع احتمال تنها در نواحی پراهمیت است. اینکه در این روش به‌جای نمونه‌برداری از همه نواحی توزیع احتمال تنها از نواحی پراهمیت یا مدنظر، نمونه‌برداری انجام می‌گیرد به دلیل کاهش حجم محاسبات یا سهولت در نمونه‌برداری می‌تواند باشد. این روش به‌ویژه برای فضاهایی با بعد بالا که معمولاً هم داده‌ها تنک بوده و ناحیه مطلوب کوچکتر از کل فضای داده‌است بسیار کاراست. ایده روش نمونه‌برداری بااهمیت، گزینش یک توزیع پیشنهادی q(x) به‌جای استفاده از توزیع احتمال اصلی P(x) است. همانگونه که گفته شد دلیل این جایگزینی به دلیل عدم لزوم یا سختی در نمونه‌برداری از توزیع احتمال P(x) است.

f(x)P(x)dx=f(x)P(x)q(x)q(x)dx

در نتیجه نمونه‌برداری بااهمیت مونت‌کارلو به استفاده از N ذره مستقل و تصادفی بیرون کشیده شده از q(x) تلقی می‌شود تا مقدار متوسط تابع f(x) را بگونه زیر تخمین بزند.

𝔼^(f(x))=1Ni=1Nw(x(i))f(x(i))

که در آن w(x) وزن بااهمیت (Importance Weight) گفته شده و به گونه زیر بدست می‌آید.

w(x(i))=P(x(i))q(x(i))

درحقیقت دلیل تقریب زدن متوسط تابع f(x) بر طبق روابط بالا در این است که توزیع q(x) (و نه P(x)) را برپایه تخمین مونت‌کارلو به گونه زیر می‌توان تقریب زد.

q(x)1Ni=1Nδ(xx(i))

که این‌بار ذره‌های {x(i)}i=1N به‌طور مستقل از هم و تصادفی از توزیع چگالی احتمال q(x) بیرون کشیده شده‌اند. سپس بر طبق رابطه [?] تخمین متوسط f(x) بگونه زیر بدست می‌آید. الگو:چپ‌چین الگو:NumBlk الگو:پایان چپ‌چین

با مقایسه تخمین مونت‌کارلو از رابطه الگو:EquationNote و تخمین با نمونه‌های با اهمیت از رابطه الگو:EquationNote فهمیده می‌شود که تخمین برپایه نمونه‌برداری بااهمیت مونت کارلو به کمک ذره‌های بیرون کشیده شده از توزیع دست یافتنی q(x) و نیز وزن‌های بااهمیت مرتبط با آن w(x(i)) انجام می‌گیرد؛ بنابراین به طریقی نمونه‌برداری از توزیع P(x) با نمونه‌برداری از توزیع قابل حصول q(x) جایگزین شده‌است. در استفاده‌های عملی از این روش مانند فیلتر ذره‌ای که پیرامون آن در همین فصل توضیح داده می‌شود از وزن‌های نرمالیزه (Normalized importance weights) شده استفاده می‌شود. تعریف وزن‌های بااهمیت نرمالیزه از رابطه زیر بدست می‌آید.

w~(x(i))=w(x(i))j=1Nw(x(j))

در این صورت تخمین متوسط تابع f(x) را به گونه زیر می‌توان نوشت [?].

f^N=1Nj=1Nw(x(j))f(x(j))j=1Nw(x(j))=j=1Nw~(x(j))f(x(j))

با توجه به رابطه [?] به راحتی می‌توان فهمید که این تخمین‌گر نیز نااریب است چراکه داریم:

𝔼q(x)(f^N)=𝔼q(x)(1Ni=1NP(x(i))q(x(i))f(x(i)))=1Ni=1N𝔼q(x)(P(x)q(x)f(x))=1Ni=1N𝔼P(x)(f(x))=𝔼(f(x))

واریانس این تخمین نیز به کمک روابط زیر محاسبه می‌شود [?].

Varq(x)(f^N)=𝔼q(x)((1Ni=1Nw(x(i))f(x(i))𝔼q(x)(P(x)q(x)f(x)))2)=1N2𝔼q(x)(i,jNP(x(i))P(x(j))q(x(i))q(x(j))f(x(i))f(x(j)))𝔼q(x)((P(x)q(x)f(x))2)=1N2((N2N)𝔼q(x)(P(x)q(x)f(x))2+N𝔼q(x)((P(x)q(x)f(x))2))𝔼q(x)(P(x)q(x)f(x))2=1N(f(x)P(x))2q(x)dx𝔼P(x)(f(x))2N

در بدست آوردن رابطه اخیر از i.i.d. بودن ذره‌های بیرون کشیده شده از توزیع q(x) استفاده شد، یعنی:

𝔼(P(x(i))P(x(j))q(x(i))q(x(j))f(x(i))f(x(j)))={𝔼((P(x)q(x)f(x))2)i=j,𝔼(P(x)q(x)f(x))2ij

با توجه به رابطه [?] مشخص می‌شود که واریانس این تخمین متفاوت از رابطه [?] مربوط به تخمین مونت‌کارلو بدست می‌آید، اما این واریانس در صورت برقراری یکی از دو شرط زیر می‌تواند کوچک شود [?].

  • اگر q(x) متناسب با P(x) برگزیده شود به نحوی که تقریبی از واریانس توزیع اصلی را بدست دهد.
  • اگر q(x) متناسب با |f(x)|P(x) برگزیده شود به نحوی که واریانس توزیع اصلی را نیز بر طبق رابطه بالا کمتر کند.

بنابراین نمونه‌برداری با اهمیت از دو جهت می‌تواند حائز اهمیت باشد، یکی اینکه نمونه‌برداری از توزیع q(x) می‌تواند به مراتب راحت‌تر و کم محاسبه‌تر از نمونه‌برداری از توزیع P(x) باشد و دیگر اینکه این نوع نمونه‌برداری می‌تواند حتی واریانس تخمین را کمتر از واریانس تخمین مونت‌کارلو کند.

اما گزینش توزیع q(x) خود یک مسئله بسیار حیاتی در این مقوله است، انتخاب نادرست این توزیع موجب افزایش واریانس تخمین‌گر شده باعث می‌شود هزینه زیادی در قبال راحت‌تر بودن نمونه‌برداری از آن در قبال توزیع P(x) بپردازیم. در حقیقت اگر q(x) مناسب برگزیده نشود توزیع وزن‌ها بسیار ناهمگون و پراکنده خواهد شد و در نتیجه بسیاری از وزن‌ها به دلیل توزیع خاصشان بی‌استفاده می‌شوند. همچنین در این حالت تخمین برای فضاهای بزرگ تنها توسط تعداد کمی از ذره‌ها با وزن‌های بزرگ تعیین‌کننده می‌شود. بنابر آنچه گفته شد، گزینش توزیع پیشنهادی q(x) باید بر طبق بندهای زیر باشد [?]:

  • توزیع q(x) باید به ازای xهایی که P(x)0 است بزرگتر از صفر باشد، q(x)0.
  • توزیع q(x) باید در نواحی مورد محاسبه متناسب با |f(x)|P(x) ویا متناسب با P(x) باشد.
  • نمونه‌برداری کردن از q(x) راحت باشد.
  • محاسبه q(x) برای هر ذره x(i) ساده باشد.

با این همه در برخی تجربه‌های عملی نشان داده شده‌است[۵] که q(x) یک شرط دیگر باید داشته باشد، و آن اینکه باسرعت کمتری نسبت به P(x) به سمت صفر بگراید تا حساسیتی کمتر نسبت به داده‌های پرت (Outlier) داشته باشد یا به عبارتی دیگر این توزیع باید دمبی (tail) بزرگ داشته باشد، یک معیار برای سنجش میزان دمب یک توزیع، تابع رتبه (score) است که رابطه آن در زیر آمده‌است.

ϕ(x)=dlog(P(x))dx

نشان داده می‌شود که اگر ϕ(x) شرط محدود بودن را داشته باشد آنگاه توزیع P(x) دارای دمبی بزرگ خواهد بود[۶]، در غیر این صورت خیر. اما محدود بودن برای تابع مزبور بدین معنی است که به ازای هر x حقیقی‌ای وجود داشته باشد کمیت‌های حقیقی L و U که داشته باشیم:

x:Lϕ(x)U

نمونه‌برداری بااهمیت تناوبی

همانگونه که در بخش پیشین بحث شد، گزینش توزیع پیشنهادی در بازدهی و کارایی تخمین برپایه نمونه‌برداری بااهمیت بسیار مهم و حیاتی است؛ بنابراین شیوه گزینش یک توزیع مناسب کلید بهره‌بری از یک تخمین‌گر موفق شمرده می‌شود. با این حال در اغلب موارد یافتن یک توزیع متناسب به‌ویژه در مواردی که با فضاهایی با بعد بالا سروکار داریم کار سختی است [?]. از این رو یک راه مقابله با این مشکل ساختن توزیع پیشنهادی به صورت تناوبی و پی‌درپی است که این ایده اصلی روش نمونه‌برداری با اهمیت تناوبی (Sequential Importance Sampling) یا به اختصار SIS محسوب می‌شود. برای این منظور بایسته‌است که توزیع پیشنهادی q(x) به نحوی برگزیده شود تا خاصیت زیر را داشته باشد [?].

q(x0:n|y0:n)=q(x0)t=1Nq(xt|x0:t1,y0:t)

برپایه این انتخاب، نمونه‌برداری بااهمیت به صورت تناوبی قابل حصول است. برای نشان دادن این مطلب طبق قانون تلسکوپی (Telescope law) در احتمالات برای توزیع P(x) می‌توان نوشت.

P(x0:n)=P(x0)P(x1|x0)...P(xn|x0,...,xn1)

نیز برای توزیع پیشنهادی q(x) طبق خاصیتی که برای آن مقرر شد داریم:

q(x0:n)=q(x0)q(x1|x0)...q(xn|x0,...,xn1)

در نتیجه وزن‌های بااهمیت را می‌توان به صورت زیر بازنویسی کرد.

w(x0:n)=P(x0)P(x1|x0)...P(xn|x0,...,xn1)q(x0)q(x1|x0)...q(xn|x0,...,xn1)

که در این صورت این وزن‌ها به‌صورت بازگشتی طبق رابطه زیر می‌توانند محاسبه شوند.

w(x0:n)=w(x0:n1)P(xn|x0:n1)q(xn|x0:n1)

از جمله مزیت‌های این روش در این است که این وزن‌ها در طول زمان به‌صورت تناوبی محاسبه می‌شوند و بازدهی آن در حین الگوریتم افزوده می‌شود، اما عیب مهم آن در این است که واریانس وزن‌ها ممکن است بزرگ باشد که موجب کاهش دقت در تخمین می‌شوند. به عبارت دیگر در مرجع [?] نشان داده شده‌است که واریانس وزن‌های بااهمیت در طول زمان افزوده شده و پس از چند تکرر (iteration) از الگوریتم تعداد بسیاری از وزن‌ها مقداری نزدیک به صفر و تنها اندکی از وزن‌ها ناصفر خواهند بود. این پدیده که به مشکل تباهیدگی وزن‌ها (Weight degeneracy problem) موسوم است از معایب مهم نمونه‌برداری تناوبی شمرده می‌شود چرا که موجب هدررفت زمان محاسبات در بروزکردن وزن‌ها از روی وزن‌های پیشین خود می‌شود. راه حل مقابله با این مشکل استفاده از متدهای بازنمونه‌برداری است؛ بنابراین بازنمونه‌برداری به عنوان یکی از گام‌های اصلی در الگوریتم SIS شمرده می‌شود که هدف از آن حذف ذره‌هایی با وزن‌های کوچک و تقویت ذره‌های با وزن بزرگ است. یا به عبارت دیگر این فرایند موجب تشدید ذره‌های با احتمال بالا و تضعیف ذره‌های با احتمال کوچک می‌شود. همراهی روش SIS با فرایند بازنمونه‌برداری به فرایند Sequential Importance Resampling موسوم است. شکل زیر فرایند بازنمونه‌برداری را به صورت نمادین نمایش می‌دهد.

فرایند بازنمونه‌برداری

در این شکل گوی‌ها معرف ذره‌ها و ابعاد آن وزن متناظر با آن را نشان می‌دهد، گوی‌های پررنگ در ابعاد مختلف هستند، گوی بزرگ بدین معنی است که وزن ذره متناظر با آن بزرگ است یا به عبارتی احتمال آن ذره بالاست، از طرف دیگر گوی کوچک به معنی اینست که وزن ذره متناظر با آن کوچک است یا به عبارتی احتمال رخداد آن ذره کم است. گوی‌های کم‌رنگ نیز پس از فرایند بازنمونه‌برداری حاصل شده‌اند. در این شکل تشدید ذره‌ها با وزن‌های بالا و تضعیف ذره‌ها با وزن‌های کم توسط فرایند بازنمونه‌برداری به خوبی نشان داده شده‌است؛ بنابراین بر طبق مطالب گفته شده الگوریتم SIR را می‌توان به صورت زیر خلاصه کرد[?].

بنابراین این گام نقشی اساسی در الگوریتم نمونه‌برداری بااهمیت ایفا می‌کند [?] چراکه از طرفی با پدیده پراکندگی ناموزون وزن‌ها که موجب اتلاف محاسباتی می‌شود مقابله کرده و از تباهیدگی وزن‌ها جلوگیری می‌کند و از طرف دیگر هنگامیکه وزن‌ها مورب (Skewed) باشند موجب گزینش وزن‌های بااهمیت‌تر می‌شود. برای همین چه بسا این گام خود موجب تغییراتی در واریانس نمونه‌ها شده و لزوماً تخمین را بهبود ندهد. برای همین این گام با ملاحظاتی انجام می‌شود که خود موجب پیدایش الگوریتم‌های مختلف برای بازنمونه‌برداری شده‌است. نکته قابل تأمل در میان این الگوریتم‌ها در این است که اگر معیار، واریانس ذره‌های تولیدی باشد روش بازنمونه‌برداری‌ای که به روش سیستماتیک (systematic) موسوم است، واریانس کمتری را ایجاد کرده و از این لحاظ ارجح است [?]. در زیر به‌صورت خلاصه این روش بیان شد.

الگو:چپ‌چین

SIR Algorithm
1. Draw N random sample {x(i)}iN from proposal distribution q(.).
2. Calculate the importance weights for each particle: w(x(i))=P(x(i))q(x(i)).
3. Normalize the importance weights to obtain w~(i) as: w~(i)=w(x(i))j=1Nw(x(j)).
4. Resampling step: Multiply/Suppress the particle x(i) with high/low importance weights w~(i) , respectively, according to P(x).

الگو:پایان چپ‌چین

  • فرض کنید ذره‌ها و وزن‌های مرتبط با آن به صورت روبرو باشد. {x(i),w(i)}i=1N
  • اینک به ازای هر ذره نخست متغیر ui به‌صورت یکنواخت در بازه [۰٬۱] تولید شده و به کمک آن متغیر u^i به صورت روبرو ایجاد می‌شود. u^i=ui+i1N
  • به ازای i=1,...,N، وند ji به‌صورت روبرو برگزیده می‌شود. ji:ji=iPr(x(i1))u^iPr(x(i))
  • ذره و وزن‌های جدید برپایه وندji گزیده می‌شوند. x^(ji)x(ji),w(ji)=1N
  • اینک ذره و وزن‌های جدید به‌صورت رندوم تولید می‌شوند: {x^(ji),w(ji)}i=1N.

برای درک بهتر نحوه عملکرد فرایند بازنمونه‌برداری سیستماتیک خواننده را به مرجع [?] ارجاع می‌دهیم.

الگوریتم فیلتر ذره‌ای

تا اینجا با اساس و پیشنیازهای لازم برای بیان جزئیات الگوریتم فیلتر ذره‌ای آشنا شدیم و نیز ذکر شد که این فیلتر براساس الگوریتم SIR استوار است. شاید ذکر این نکته نیز برای شروع معرفی فیلتر ذره‌ای مناسب باشد که این فیلتر با تقریب زدن چگالی پسین به‌صورت تناوبی و براساس تخمین بیزین عمل می‌کند [?]. در حقیقت چون نمونه‌برداری مستقیم از چگالی پسین عملی نیست باید به دنبال راهکاری بود که این چگالی را تخمین زد، فیلتر ذره‌ای گزیری برای رفع این مشکل است. اینک توضیح این فیلتر را با تخمین چگالی پسین به کمک ذره‌های بیرون کشیده شده از آن آغاز می‌کنیم [?]. همانگونه که بیان شد بر طبق تخمین مونت‌کارلو چگالی پسین را برپایه N ذره بیرون کشیده شده از آن به‌صورت زیر می‌توان تقریب زد.

P^(𝒙𝒌|𝒀𝒌)=1Ni=1Nδ(𝒙𝒌x𝒌(𝒊))

که در آن {x(i)}i=1N ذره‌های مستقل از هم و به‌صورت تصادفی بیرون کشیده شده از چگالی پسین است. بر پایه این تخمین امید ریاضی هر تابع دلخواه f(.) را به‌صورت زیر می‌توان تقریب زد.

𝔼(f(𝒙𝒌))1Ni=1Nf(x𝒌(𝒊))

با این حال نمونه‌برداری از چگالی مزبور تنها حالتی خاص و غیرپیچیده در تخمین مونت‌کارلو است، در حالت کلی چگالی پسین دیگری در استخراج فیلتر ذره‌ای بکارگرفته می‌شود که چگالی کامل‌تری در روش فیلترکردن بیزین محسوب شده، به دلیل استفاده از اطلاعات افزون‌تر عملکرد بهتری خواهد داشت. این چگالی برابر P(𝑿𝒌|𝒀𝒌) بوده که در آن 𝑿𝒌 مجموعه بردارهای حالت از ابتدا تا لحظه kاست.

𝑿𝒌={𝒙0,...,𝒙𝒌}

اما به دلیل اینکه نمونه‌برداری مستقیم از این چگالی قابل عملی نیست بر طبق روش نمونه‌برداری بااهمیت از یک توزیع پیشنهادی q(.) بهره گرفته می‌شود؛ بنابراین تخمین تابع دلخواهf(.) با این توزیع پیشنهادی به صورت زیر بدست می‌آید.

𝔼(f(𝒙𝒌))=f(𝒙𝒌)P(𝑿𝒌|𝒀𝒌)d𝒙𝒌=f(𝒙𝒌)P(𝒀𝒌|𝑿𝒌)P(𝑿𝒌)q(𝑿𝒌|𝒀𝒌)P(𝒀𝒌)q(𝑿𝒌|𝒀𝒌)d𝒙𝒌=f(𝒙𝒌)Wk(𝑿𝒌)q(𝑿𝒌|𝒀𝒌)d𝒙𝒌P(𝒀𝒌|𝑿𝒌)P(Xk)q(𝑿𝒌|𝒀𝒌)q(𝑿𝒌|𝒀𝒌)d𝒙𝒌=f(𝒙𝒌)Wk(𝑿𝒌)q(𝑿𝒌|𝒀𝒌)d𝒙𝒌Wk(𝑿𝒌)q(𝑿𝒌|𝒀𝒌)d𝒙𝒌=𝔼q(f(𝒙𝒌)Wk(𝒙𝒌)|𝒀𝒌)𝔼q(Wk(𝒙𝒌)|𝒀𝒌))

اینک چون مخرج رابطه اخیر مقداری ثابت است آن را می‌توان به داخل اپراتور متوسط‌گیری صورت برد، در نتیجه می‌توان نوشت:

𝔼q(f(𝒙𝒌))=𝔼q(f(𝒙𝒌)Wk(𝑿𝒌)𝔼q(Wk(𝑿𝒌|𝒀𝒌))|𝒀𝒌)

همان‌طور که پیشتر گفته شد، کمیت Wk(𝑿𝒌)، وزن بااهمیت نامیده شده بر طبق رابطه زیر بدست می‌آید.

Wk(𝑿𝒌)=P(𝒀𝒌|𝑿𝒌)P(𝑿𝒌)q(𝑿𝒌|𝒀𝒌)

همچنین کمیت نرمالیزه شده وزن‌های بااهمیت با W~k(𝑿𝒌) نشان داده شد که از رابطه زیر بدست می‌آید.

W~k(𝑿𝒌)=Wk(𝑿𝒌)𝔼q(Wk(𝑿𝒌)|𝒀𝒌)

پس به عبارتی چگالی پسین P(𝑿𝒌|𝒀𝒌) بر پایه تخمین نمونه‌برداری بااهمیت مونت‌کارلو از رابطه زیر تخمین زده می‌شود.

P^(𝑿𝒌|𝒀𝒌)i=1Nw~k(i)δ(𝑿𝒌X𝒌(𝒊))

که در آن ذره‌های {Xk(i)}iN به‌جای آنکه از چگالی نامعلوم و غیرقابل عملی P(𝑿𝒌|𝒀𝒌) نمونه‌برداری شده باشند از تویع پیشنهادی q(𝑿𝒌|𝒀𝒌) بیرون کشیده شده‌اند. وزن‌های نرمالیزه شده آن یعنی {w~k(i)}iN نیز طبق رابطه زیر بدست می‌آیند.

w~k(i)=W~k(𝑿𝒌)|𝑿𝒌=X𝒌(𝒊)

پیرامون شرایط حیاتی انتخاب توزیع پیشنهادی در بخش [?] بحث شد، یکی از این شروط این بود که آن را بتوان به صورت زیر تجزیه کرد.

q(𝑿𝒌|𝒀𝒌)=q(𝑿𝒌1|𝒀𝒌1)q(𝒙𝒌|𝑿𝒌1,𝒀𝒌)

در نتیجه با کمک گرفتن از این تجزیه و با درنظر گرفتن دو شرط اساسی \eqref{ASUM1} و \eqref{ASUM2} وزن‌های با اهمیت را به صورت تناوبی می‌توان محاسبه کرد، برای نشان دادن این مطلب توجه کنید که:

Wk(𝑿𝒌)=P(𝒀𝒌|𝑿𝒌)P(𝑿𝒌)q(𝑿𝒌|𝒀𝒌)=P(𝒚𝒌,𝒀𝒌1|𝑿𝒌)P(𝒙𝒌|𝑿𝒌1)P(𝑿𝒌1)q(𝑿𝒌1|𝒀𝒌1)q(𝒙𝒌|𝑿𝒌1,𝒀𝒌)=P(𝒚𝒌|𝒀𝒌1,𝑿𝒌)P(𝒀𝒌1|𝑿𝒌)P(𝒙𝒌|𝑿𝒌1)P(𝑿𝒌1)q(𝑿𝒌1|𝒀𝒌1)q(𝒙𝒌|𝑿𝒌1,𝒀𝒌)

اینک با درنظر داشتن این نکته که با داشتن بردار حالت در لحظه k1 و پیش از آن اطلاعات افزون‌تری در احتمال رخداد 𝒚𝒌 بدست نمی‌دهند، داریم:

P(𝒚𝒌|𝒀𝒌1,𝑿𝒌)=P(𝒚𝒌|𝑿𝒌)=P(𝒚𝒌|𝒙𝒌)

از طرف دیگر به دلیل علیت می‌دانیم بردار حالت در زمان k تأثیری در احتمال رخداد بردار مشاهده در لحظه k1 ندارد. پس می‌توان نوشت:

P(𝒀𝒌1|𝑿𝒌)=P(𝒀𝒌1|𝑿𝒌1)

بنابراین در مجموع رابطه زیر را به‌صورتی که درپی می‌آید می‌توان ادامه داد.

Wk(𝑿𝒌)=P(𝒚𝒌|𝒀𝒌1,𝑿𝒌)P(𝒀𝒌1|𝑿𝒌)P(𝒙𝒌|𝑿𝒌1)P(𝑿𝒌1)q(𝑿𝒌1|𝒀𝒌1)q(𝒙𝒌|𝑿𝒌1,𝒀𝒌)=Wk1(𝑿𝒌1)P(𝒚𝒌|𝒙𝒌)P(𝒙𝒌|𝒙𝒌1)q(𝒙𝒌|𝑿𝒌1,𝒀𝒌)

بنابراین وزن‌های متناظر با ذره‌های نمونه‌برداری شده از توزیع پیشنهادی به صورت تناوبی زیر قابل دستیابی هستند.

Wk(x𝒌(𝒊))=Wk1(x𝒌1(𝒊))P(𝒚𝒌|x𝒌(𝒊))P(x𝒌(𝒊)|x𝒌1(𝒊))q(x𝒌(𝒊)|X𝒌1(𝒊),𝒀𝒌)

که در آن X𝒌(𝒊) مجموعه ذره‌های نمونه برداری شده از آغاز تا لحظه حاضر است، یعنی:

X𝒌(𝒊)={x𝒌(𝒊)}i=1N

بنابراین فیلتر ذره‌ای برپایه نمونه‌برداری تناوبی مونت‌کارلو عمل کرده و به کمک ذره‌های بیرون کشیده شده از توزیع پیشنهادی و وزن‌های متناظر با آن که به صورت تناوبی از وزن‌های ماقبل خود بدست می‌آیند تخمینی از بردار حالت ارائه می‌دهد. شاید تکرار این نکته مفید باشد که طبق رابطه SIS وزن‌های بااهمیت در هر لحظه به کمک ذره‌های نمونه برداشته شده از توزیع پیشنهادی و ارزیابی احتمال درستنمایی و چگالی گذار که بر طبق فضای حالت دانسته هستند بدست می‌آیند؛ بنابراین تخمین MMSE به کمک این ذره‌ها و وزن‌های متناظر با آن به طریق زیر محاسبه می‌شوند.

𝒙^k=𝔼P(𝑿𝒌|𝒀𝒌)(𝒙𝒌)𝒙𝒌P^(𝑿𝒌|𝒀𝒌)d𝒙𝒌=i=1Nw~k(i)𝒙𝒌(i)

در نهایت یادآوری این نکته حائز اهمیت است که واریانس وزن‌ها در طول الگوریتم رو به فزونی می‌گذارد، بگونه‌ای که مقدار بیشتر وزن‌ها به سمت صفر گراییده، مقدار تعداد اندکی از آنها بزرگ می‌شود، و درنهایت موجب تباه شدن وزن‌ها می‌شود که لازم می‌دارد از متدهای بازنمونه‌برداری که در بخش [?] پیرامون آن بحث شد استفاده شود. بنابر مطالب گفته شده الگوریتم فیلتر ذرهای در جدول زیر خلاصه شد.

الگو:چپ‌چین

Particle Filter Algorithm
Initialization: k=0
For i=1,... ,N:الگو:سخ

Draw the particles x0(𝒊) from density P(𝒙0).الگو:سخ End

For k=1,2,...
Importance sampling step:الگو:سخ

For i=1,... ,Nالگو:سخ sample x^𝒌(𝒊) from q(x𝒌(𝒊)|X𝒌1(𝒊),𝒀𝒌)الگو:سخ End

For i=1,... ,Nالگو:سخ

Evaluate the importance weights: wk(i)=wk1(i)P(𝒚𝒌|x^𝒌(𝒊))P(x^𝒌(𝒊)|x𝒌1(𝒊))q(x^𝒌(𝒊)|X𝒌1(𝒊),𝒀𝒌)الگو:سخ End

For i=1,... ,Nالگو:سخ

Normalize the importance weights: w~k(i)=wk(i)i=1Nwk(i)الگو:سخ End

Resampling step:الگو:سخ

Eliminate/Multiply particles x^𝒌(𝒊) according to its normalized weights w~k(i) to obtain particles x𝒌(𝒊)

Output: the optimal MMSE estimation is given as: 𝒙^𝒌1N𝒊=1𝑵𝒙𝒌(𝒊)
End

الگو:پایان چپ‌چین

نکته دیگری که باید گفته شود مربوط به گام آغازین الگوریتم است. در این لحظه ذره‌ها از روی چگالی P(𝒙0) نمونه‌برداری می‌شوند که این چگالی احتمال با در نظر گرفتن تحرک جسم در لحظه اولیه تقریب زده می‌شود. تا اینجا درمورد نحوه انتخاب توزیع پیشنهادی نکاتی بیان شد، اما گزینه‌هایی مطرح نشد. بر طبق مرجع [?] انتخابی بهینه از دید کمینه شدن واریانس وزن‌های بااهمیت طبق گزینه زیر است.

q(𝒙𝒌|𝑿𝒌1,𝒀𝒌)=P(𝒙𝒌|𝑿𝒌1,𝒀𝒌)

اما همان‌طور هم که در مرجع [?] مذکور بیان شده‌است نمونه‌برداری از این توزیع غیرعملی است. گزینه دیگر \cite{transition} استفاده از چگالی گذار است، یعنی:

q(𝒙𝒌|𝑿𝒌1,𝒀𝒌)=P(𝒙𝒌|𝒙𝒌1)

با این انتخاب در هر iteration وزن‌های بااهمیت بر طبق رابطه [?] به صورت زیر بروز می‌شوند.

Wk(i)=Wk1(i)P(𝒚𝒌|x𝒌(𝒊))

بدین معنی که تنها دانستن احتمال درستنمایی در بروزکردن وزن‌ها کفایت می‌کند. روش دیگر تقریب زدن توزیع پیشنهادی با یک توزیع گاوسی است. انتخاب برای توزیع q(𝒙𝒌|𝑿𝒌1,𝒀𝒌) می‌بایست حول احتمال رخداد بردار حالت در لحظه k به شرط در اختیار داشتن کل مجموعه بردارهای مشاهدات تا لحظه k و کل مجموعه بردارهای حالت تا لحظه k1 باشد. به عبارت دیگر باید دید چنانچه بردار مشاهدات و بردار حالت بترتیب تا لحظه k و k1 در دسترس باشند، بردار حالت در لحظه k چه توزیعی خواهد داشت آن توزیع را برای توزیع پیشنهادی برگزید. بنابر آنچه گفته شد می‌توان این توزیع را جداگانه به کمک فیلتر کالمن و با فرض گاوسی بودن این توزیع بدست آورد [?]. یعنی با انتشار Propagate ذره‌ها از لحظه پیشین به لحظه کنونی و به یاری معادلات کالمن این توزیع را تخمین زد. در زیر نحوه بدست آوردن میانگین و کواریانس این توزیع توضیح داده می‌شود. نخست از روی ذره‌های کنونی و به کمک ماتریس گذار، ذره‌های پیش‌بینی طبق فرایند پیش‌بینی کالمن (Kalman filter prediction procedure) زیر بدست آورده می‌شوند.

x¯𝒌/𝒌1(𝒊)=𝑨𝒌x𝒌(𝒊)

پس از آن ماتریس MSE پیش‌بینی کمینه(minimum prediction MSE matrix) به کمک کواریانس نویز پروسس محاسبه می‌شود.

P𝒌/𝒌1(𝒊)=𝑨𝒌P𝒌1(𝒊)A𝒌𝑻+R𝒌𝒖

که در آن همان‌طور که گفته شد 𝑹𝒖 ماتریس کواریانس نویز پروسس بوده از رابطه زیر حساب می‌شود.

R𝒌𝒖=𝔼(𝒖𝒌u𝒌𝑻)

گین کالمن (Kalman gain) بر پایه معادله زیر و به یاری ماتریس کواریانس نویز اندازه‌گیری بدست می‌آید.

K𝒌(𝒊)=P𝒌/𝒌1(𝒊)H𝒌𝑻(R𝒌𝒘+𝑯𝒌P𝒌/𝒌1(𝒊)H𝒌𝑻)1R𝒌𝒘=𝔼(𝒘𝒌w𝒌𝑻)

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

x¯𝒌(𝒊)=x¯𝒌/𝒌1(𝒊)+K𝒌(𝒊)(𝒚𝒌𝑯𝒌x¯𝒌/𝒌1(𝒊))P𝒌(𝒊)=P𝒌/𝒌1(𝒊)K𝒌(𝒊)𝑯𝒌P𝒌/𝒌1(𝒊)

به عبارت دیگر انتخاب توزیع پیشنهادی به صورت زیر خواهد بود، که در آن میانگین این توزیع برابر بردار حالت اصلاح شده بوده و ماتریس کواریانس آن همان ماتریس MMSE است.

q(x𝒌(𝒊)|X𝒌1(𝒊),𝒀𝒌)=𝒩(x¯𝒌(𝒊),P𝒌(𝒊)),i=1,...,N

بنابراین گام نمونه‌برداری بااهمیت از الگوریتم فیلترذره‌ای در جدول زیر خلاصه شد.

الگو:چپ‌چین

Importance Sampling with Kalman Filter
The importance sampling step in the particle filter:
For i=1,... ,N
Update the particles by the Kalman filter:الگو:سخ

x¯𝒌/𝒌1(𝒊)=𝑨𝒌x𝒌(𝒊)الگو:سخ P𝒌/𝒌1(𝒊)=𝑨𝒌P𝒌1(𝒊)A𝒌𝑻+R𝒌𝒖الگو:سخ K𝒌(𝒊)=P𝒌/𝒌1(𝒊)H𝒌𝑻(R𝒌𝒘+𝑯𝒌P𝒌/𝒌1(𝒊)H𝒌𝑻)1الگو:سخ x¯𝒌(𝒊)=x¯𝒌/𝒌1(𝒊)+K𝒌(𝒊)(𝒚𝒌𝑯𝒌x¯𝒌/𝒌1(𝒊))الگو:سخ P𝒌(𝒊)=P𝒌/𝒌1(𝒊)K𝒌(𝒊)𝑯𝒌P𝒌/𝒌1(𝒊)الگو:سخ Sample x^𝒌(𝒊) from 𝒩(x¯𝒌(𝒊),P𝒌(𝒊))

End

الگو:پایان چپ‌چین

به الگوریتم فیلتر ذره‌ای که توزیع پیشنهادی آن توسط فیلتر کالمن تخمین زده شود، فیلتر ذره‌ای کالمن (Kalman Particle Filter) یا به اختصار KPF می‌گویند.

منابع

الگو:پانویسالگو:یادکرد-ویکی

  1. A. Doucet and X.Wang, “Monte Carlo methods for signal processing: A review in the statistical signal processing context”, IEEE Signal Process. Mag. , vol. 22, no. 6, pp. 152–170, Nov. 2005
  2. Z. Chen, “Bayesian filtering: From kalman filters to particle filters, and beyond”, Technical report, McMaster University, 2003.
  3. S. M. Kay, “Fundamentals of statistical signal processing: estimation theory”, Englewood Cliffs, NJ: Prentice Hall International, Inc. , 1993.
  4. ?
  5. D. J. C. Mackay, “Introduction to Monte Carlo methods”, In Proceedings of the NATO Advanced Study Institute on Learning in graphical models, pp.175-204, March 1998.
  6. Z. Chen, “Bayesian filtering: From kalman filters to particle filters, and beyond”, Technical report, McMaster University, 2003.