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

تا اینجا با اساس و پیشنیازهای لازم برای بیان جزئیات الگوریتم فیلتر ذره‌ای آشنا شدیم و نیز ذکر شد که این فیلتر براساس الگوریتم SIR استوار است. شاید ذکر این نکته نیز برای شروع معرفی فیلتر ذره‌ای مناسب باشد که این فیلتر با تقریب زدن چگالی پسین به‌صورت تناوبی و براساس تخمین بیزین عمل می‌کند [?]. در حقیقت چون نمونه‌برداری مستقیم از چگالی پسین عملی نیست باید به دنبال راهکاری بود که این چگالی را تخمین زد، فیلتر ذره‌ای گزیری برای رفع این مشکل است.
اینک توضیح این فیلتر را با تخمین چگالی پسین به کمک ذره‌های بیرون کشیده شده از آن آغاز می‌کنیم [?]. همانگونه که بیان شد بر طبق تخمین مونت‌کارلو چگالی پسین را برپایه  N ذره بیرون کشیده شده از آن به‌صورت زیر می‌توان تقریب زد.
{\hat {P}}({\boldsymbol {x_{k}|Y_{k}}})={\frac {1}{N}}\sum _{i=1}^{N}{\delta ({\boldsymbol {x_{k}-x_{k}^{(i)}}})}
که در آن   \{x^{(i)}\}_{i=1}^{N} ذره‌های مستقل از هم و به‌صورت تصادفی بیرون کشیده شده از چگالی پسین است. بر پایه این تخمین امید ریاضی هر تابع دلخواه  f(.) را به‌صورت زیر می‌توان تقریب زد.
\mathbb {E} (f({\boldsymbol {x_{k}}}))\approx {\frac {1}{N}}\sum _{i=1}^{N}{f({\boldsymbol {x_{k}^{(i)}}})}
با این حال نمونه‌برداری از چگالی مزبور تنها حالتی خاص و غیرپیچیده در تخمین مونت‌کارلو است، در حالت کلی چگالی پسین دیگری در استخراج فیلتر ذره‌ای بکارگرفته می‌شود که چگالی کامل‌تری در روش فیلترکردن بیزین محسوب شده، به دلیل استفاده از اطلاعات افزون‌تر عملکرد بهتری خواهد داشت. این چگالی برابر   P({\boldsymbol {X_{k}|Y_{k}}}) بوده که در آن  {\boldsymbol {X_{k}}} مجموعه بردارهای حالت از ابتدا تا لحظه  kاست.
{\boldsymbol {X_{k}=\{x_{0},...,x_{k}\}}}
اما به دلیل اینکه نمونه‌برداری مستقیم از این چگالی قابل عملی نیست بر طبق روش نمونه‌برداری بااهمیت از یک توزیع پیشنهادی  q(.) بهره گرفته می‌شود؛ بنابراین تخمین تابع دلخواه f(.) با این توزیع پیشنهادی به صورت زیر بدست می‌آید.
{\begin{aligned}\mathbb {E} (f({\boldsymbol {x_{k}}}))&=\int f({\boldsymbol {x_{k}}})P({\boldsymbol {X_{k}|Y_{k}}})d{\boldsymbol {x_{k}}}=\int f({\boldsymbol {x_{k}}}){\frac {P({\boldsymbol {Y_{k}|X_{k}}})P({\boldsymbol {X_{k}}})}{q({\boldsymbol {X_{k}|Y_{k}}})P({\boldsymbol {Y_{k}}})}}q({\boldsymbol {X_{k}|Y_{k}}})d{\boldsymbol {x_{k}}}\\&={\frac {\int f({\boldsymbol {x_{k}}})W_{k}({\boldsymbol {X_{k}}})q({\boldsymbol {X_{k}|Y_{k}}})d{\boldsymbol {x_{k}}}}{\int P({\boldsymbol {Y_{k}|X_{k}}})P(X_{k}){\frac {q({\boldsymbol {X_{k}|Y_{k}}})}{q({\boldsymbol {X_{k}|Y_{k}}})}}d{\boldsymbol {x_{k}}}}}={\frac {\int f({\boldsymbol {x_{k}}})W_{k}({\boldsymbol {X_{k}}})q({\boldsymbol {X_{k}|Y_{k}}})d{\boldsymbol {x_{k}}}}{\int W_{k}({\boldsymbol {X_{k}}})q({\boldsymbol {X_{k}|Y_{k}}})d{\boldsymbol {x_{k}}}}}\\&={\frac {\mathbb {E} _{q}(f({\boldsymbol {x_{k}}})W_{k}({\boldsymbol {x_{k}}})|{\boldsymbol {Y_{k}}})}{\mathbb {E} _{q}(W_{k}({\boldsymbol {x_{k}}})|{\boldsymbol {Y_{k}}}))}}\end{aligned}}
اینک چون مخرج رابطه اخیر مقداری ثابت است آن را می‌توان به داخل اپراتور متوسط‌گیری صورت برد، در نتیجه می‌توان نوشت:
\mathbb {E} _{q}(f({\boldsymbol {x_{k}}}))=\mathbb {E} _{q}(f({\boldsymbol {x_{k}}}){\frac {W_{k}({\boldsymbol {X_{k}}})}{\mathbb {E} _{q}(W_{k}({\boldsymbol {X_{k}}}|{\boldsymbol {Y_{k}}}))}}|{\boldsymbol {Y_{k}}})
همان‌طور که پیشتر گفته شد، کمیت  W_{k}({\boldsymbol {X_{k}}})، وزن بااهمیت نامیده شده بر طبق رابطه زیر بدست می‌آید.
W_{k}({\boldsymbol {X_{k}}})={\frac {P({\boldsymbol {Y_{k}|X_{k}}})P({\boldsymbol {X_{k}}})}{q({\boldsymbol {X_{k}|Y_{k}}})}}
همچنین کمیت نرمالیزه شده وزن‌های بااهمیت با  {\tilde {W}}_{k}({\boldsymbol {X_{k}}}) نشان داده شد که از رابطه زیر بدست می‌آید.
{\tilde {W}}_{k}({\boldsymbol {X_{k}}})={\frac {W_{k}({\boldsymbol {X_{k}}})}{\mathbb {E} _{q}(W_{k}({\boldsymbol {X_{k}}})|{\boldsymbol {Y_{k}}})}}
پس به عبارتی چگالی پسین  P({\boldsymbol {X_{k}|Y_{k}}}) بر پایه تخمین نمونه‌برداری بااهمیت مونت‌کارلو از رابطه زیر تخمین زده می‌شود.
{\hat {P}}({\boldsymbol {X_{k}|Y_{k}}})\approx \sum _{i=1}^{N}{{\tilde {w}}_{k}^{(i)}\delta ({\boldsymbol {X_{k}-X_{k}^{(i)}}})}
که در آن ذره‌های  \{X_{k}^{(i)}\}_{i}^{N} به‌جای آنکه از چگالی نامعلوم و غیرقابل عملی  P({\boldsymbol {X_{k}|Y_{k}}}) نمونه‌برداری شده باشند از تویع پیشنهادی  q({\boldsymbol {X_{k}|Y_{k}}}) بیرون کشیده شده‌اند. وزن‌های نرمالیزه شده آن یعنی  \{{\tilde {w}}_{k}^{(i)}\}_{i}^{N} نیز طبق رابطه زیر بدست می‌آیند.
{\tilde {w}}_{k}^{(i)}={\tilde {W}}_{k}({\boldsymbol {X_{k}}})|_{{\boldsymbol {X_{k}}}={\boldsymbol {X_{k}^{(i)}}}}
پیرامون شرایط حیاتی انتخاب توزیع پیشنهادی در بخش [?] بحث شد، یکی از این شروط این بود که آن را بتوان به صورت زیر تجزیه کرد.
q({\boldsymbol {X_{k}|Y_{k}}})=q({\boldsymbol {X_{k-1}|Y_{k-1}}})q({\boldsymbol {x_{k}|X_{k-1},Y_{k}}})
در نتیجه با کمک گرفتن از این تجزیه و با درنظر گرفتن دو شرط اساسی \eqref{ASUM1} و \eqref{ASUM2} وزن‌های با اهمیت را به صورت تناوبی می‌توان محاسبه کرد، برای نشان دادن این مطلب توجه کنید که:
{\begin{aligned}W_{k}({\boldsymbol {X_{k}}})&={\frac {P({\boldsymbol {Y_{k}|X_{k}}})P({\boldsymbol {X_{k}}})}{q({\boldsymbol {X_{k}|Y_{k}}})}}={\frac {P({\boldsymbol {y_{k},Y_{k-1}|X_{k}}})P({\boldsymbol {x_{k}|X_{k-1}}})P({\boldsymbol {X_{k-1}}})}{q({\boldsymbol {X_{k-1}|Y_{k-1}}})q({\boldsymbol {x_{k}|X_{k-1},Y_{k}}})}}\\&={\frac {P({\boldsymbol {y_{k}|Y_{k-1},X_{k}}})P({\boldsymbol {Y_{k-1}|X_{k}}})P({\boldsymbol {x_{k}|X_{k-1}}})P({\boldsymbol {X_{k-1}}})}{q({\boldsymbol {X_{k-1}|Y_{k-1}}})q({\boldsymbol {x_{k}|X_{k-1},Y_{k}}})}}\end{aligned}}
اینک با درنظر داشتن این نکته که با داشتن بردار حالت در لحظه  k-1 و پیش از آن اطلاعات افزون‌تری در احتمال رخداد   {\boldsymbol {y_{k}}} بدست نمی‌دهند، داریم:
P({\boldsymbol {y_{k}|Y_{k-1},X_{k}}})=P({\boldsymbol {y_{k}|X_{k}}})=P({\boldsymbol {y_{k}|x_{k}}})
از طرف دیگر به دلیل علیت می‌دانیم بردار حالت در زمان  k تأثیری در احتمال رخداد بردار مشاهده در لحظه  k-1 ندارد. پس می‌توان نوشت:
P({\boldsymbol {Y_{k-1}|X_{k}}})=P({\boldsymbol {Y_{k-1}|X_{k-1}}})
بنابراین در مجموع رابطه زیر را به‌صورتی که درپی می‌آید می‌توان ادامه داد.
{\begin{aligned}W_{k}({\boldsymbol {X_{k}}})&={\frac {P({\boldsymbol {y_{k}|Y_{k-1},X_{k}}})P({\boldsymbol {Y_{k-1}|X_{k}}})P({\boldsymbol {x_{k}|X_{k-1}}})P({\boldsymbol {X_{k-1}}})}{q({\boldsymbol {X_{k-1}|Y_{k-1}}})q({\boldsymbol {x_{k}|X_{k-1},Y_{k}}})}}\\&=W_{k-1}({\boldsymbol {X_{k-1}}}){\frac {P({\boldsymbol {y_{k}|x_{k}}})P({\boldsymbol {x_{k}|x_{k-1}}})}{q({\boldsymbol {x_{k}|X_{k-1},Y_{k}}})}}\end{aligned}}
بنابراین وزن‌های متناظر با ذره‌های نمونه‌برداری شده از توزیع پیشنهادی به صورت تناوبی زیر قابل دستیابی هستند.
W_{k}({\boldsymbol {x_{k}^{(i)}}})=W_{k-1}({\boldsymbol {x_{k-1}^{(i)}}}){\frac {P({\boldsymbol {y_{k}|x_{k}^{(i)}}})P({\boldsymbol {x_{k}^{(i)}|x_{k-1}^{(i)}}})}{q({\boldsymbol {x_{k}^{(i)}|X_{k-1}^{(i)},Y_{k}}})}}
که در آن  {\boldsymbol {X_{k}^{(i)}}} مجموعه ذره‌های نمونه برداری شده از آغاز تا لحظه حاضر است، یعنی:
{\boldsymbol {X_{k}^{(i)}}}=\{{\boldsymbol {x_{k}^{(i)}}}\}_{i=1}^{N}
بنابراین فیلتر ذره‌ای برپایه نمونه‌برداری تناوبی مونت‌کارلو عمل کرده و به کمک ذره‌های بیرون کشیده شده از توزیع پیشنهادی و وزن‌های متناظر با آن که به صورت تناوبی از وزن‌های ماقبل خود بدست می‌آیند تخمینی از بردار حالت ارائه می‌دهد. شاید تکرار این نکته مفید باشد که طبق رابطه SIS وزن‌های بااهمیت در هر لحظه به کمک ذره‌های نمونه برداشته شده از توزیع پیشنهادی و ارزیابی احتمال درستنمایی و چگالی گذار که بر طبق فضای حالت دانسته هستند بدست می‌آیند؛ بنابراین تخمین MMSE به کمک این ذره‌ها و وزن‌های متناظر با آن به طریق زیر محاسبه می‌شوند.
{\hat {\boldsymbol {x}}}_{k}=\mathbb {E} _{P({\boldsymbol {X_{k}|Y_{k}}})}({\boldsymbol {x_{k}}})\approx \int {\boldsymbol {x_{k}}}{\hat {P}}({\boldsymbol {X_{k}|Y_{k}}})d{\boldsymbol {x_{k}}}=\sum _{i=1}^{N}{{\tilde {w}}_{k}^{(i)}{\boldsymbol {x_{k}}}^{(i)}}
در نهایت یادآوری این نکته حائز اهمیت است که واریانس وزن‌ها در طول الگوریتم رو به فزونی می‌گذارد، بگونه‌ای که مقدار بیشتر وزن‌ها به سمت صفر گراییده، مقدار تعداد اندکی از آنها بزرگ می‌شود، و درنهایت موجب تباه شدن وزن‌ها می‌شود که لازم می‌دارد از متدهای بازنمونه‌برداری که در بخش [?] پیرامون آن بحث شد استفاده شود. بنابر مطالب گفته شده الگوریتم فیلتر ذرهای در جدول زیر خلاصه شد.
نکته دیگری که باید گفته شود مربوط به گام آغازین الگوریتم است. در این لحظه ذره‌ها از روی چگالی   P({\boldsymbol {x_{0}}}) نمونه‌برداری می‌شوند که این چگالی احتمال با در نظر گرفتن تحرک جسم در لحظه اولیه تقریب زده می‌شود.
تا اینجا درمورد نحوه انتخاب توزیع پیشنهادی نکاتی بیان شد، اما گزینه‌هایی مطرح نشد. بر طبق مرجع [?] انتخابی بهینه از دید کمینه شدن واریانس وزن‌های بااهمیت طبق گزینه زیر است.
q({\boldsymbol {x_{k}|X_{k-1},Y_{k}}})=P({\boldsymbol {x_{k}|X_{k-1},Y_{k}}})
اما همان‌طور هم که در مرجع [?] مذکور بیان شده‌است نمونه‌برداری از این توزیع غیرعملی است. گزینه دیگر \cite{transition} استفاده از چگالی گذار است، یعنی:
q({\boldsymbol {x_{k}|X_{k-1},Y_{k}}})=P({\boldsymbol {x_{k}|x_{k-1}}})
با این انتخاب در هر iteration وزن‌های بااهمیت بر طبق رابطه [?] به صورت زیر بروز می‌شوند.
W_{k}^{(i)}=W_{k-1}^{(i)}P({\boldsymbol {y_{k}|x_{k}^{(i)}}})
بدین معنی که تنها دانستن احتمال درستنمایی در بروزکردن وزن‌ها کفایت می‌کند. روش دیگر تقریب زدن توزیع پیشنهادی با یک توزیع گاوسی است. انتخاب برای توزیع  q({\boldsymbol {x_{k}|X_{k-1},Y_{k}}}) می‌بایست حول احتمال رخداد بردار حالت در لحظه  k به شرط در اختیار داشتن کل مجموعه بردارهای مشاهدات تا لحظه  k و کل مجموعه بردارهای حالت تا لحظه  k-1 باشد. به عبارت دیگر باید دید چنانچه بردار مشاهدات و بردار حالت بترتیب تا لحظه  k و  k-1 در دسترس باشند، بردار حالت در لحظه  k چه توزیعی خواهد داشت آن توزیع را برای توزیع پیشنهادی برگزید. بنابر آنچه گفته شد می‌توان این توزیع را جداگانه به کمک فیلتر کالمن و با فرض گاوسی بودن این توزیع بدست آورد [?]. یعنی با انتشار Propagate ذره‌ها از لحظه پیشین به لحظه کنونی و به یاری معادلات کالمن این توزیع را تخمین زد. در زیر نحوه بدست آوردن میانگین و کواریانس این توزیع توضیح داده می‌شود.
نخست از روی ذره‌های کنونی و به کمک ماتریس گذار، ذره‌های پیش‌بینی طبق فرایند پیش‌بینی کالمن (Kalman filter prediction procedure) زیر بدست آورده می‌شوند.
{\boldsymbol {{\bar {x}}_{k/k-1}^{(i)}=A_{k}x_{k}^{(i)}}}
پس از آن ماتریس MSE پیش‌بینی کمینه(minimum prediction MSE matrix) به کمک کواریانس نویز پروسس محاسبه می‌شود.
{\boldsymbol {P_{k/k-1}^{(i)}=A_{k}P_{k-1}^{(i)}A_{k}^{T}+R_{k}^{u}}}
که در آن همان‌طور که گفته شد  {\boldsymbol {R^{u}}} ماتریس کواریانس نویز پروسس بوده از رابطه زیر حساب می‌شود.
{\boldsymbol {R_{k}^{u}}}=\mathbb {E} ({\boldsymbol {u_{k}u_{k}^{T}}})
گین کالمن (Kalman gain) بر پایه معادله زیر و به یاری ماتریس کواریانس نویز اندازه‌گیری بدست می‌آید.
{\begin{array}{c}{\boldsymbol {K_{k}^{(i)}=P_{k/k-1}^{(i)}H_{k}^{T}(R_{k}^{w}+H_{k}P_{k/k-1}^{(i)}H_{k}^{T})^{-1}}}\\{\boldsymbol {R_{k}^{w}}}=\mathbb {E} ({\boldsymbol {w_{k}w_{k}^{T}}})\end{array}}
سرانجام میانگین و ماتریس کواریانس توزیع پیشنهادی بترتیب به کمک ماتریس MMSE و فرایند تصحیح فیلتر کالمن طبق روابط زیر تعیین می‌شوند.
{\begin{array}{c}{\boldsymbol {{\bar {x}}_{k}^{(i)}={\bar {x}}_{k/k-1}^{(i)}+K_{k}^{(i)}(y_{k}-H_{k}{\bar {x}}_{k/k-1}^{(i)})}}\\{\boldsymbol {P_{k}^{(i)}=P_{k/k-1}^{(i)}-K_{k}^{(i)}H_{k}P_{k/k-1}^{(i)}}}\end{array}}
به عبارت دیگر انتخاب توزیع پیشنهادی به صورت زیر خواهد بود، که در آن میانگین این توزیع برابر بردار حالت اصلاح شده بوده و ماتریس کواریانس آن همان ماتریس MMSE است.
q({\boldsymbol {x_{k}^{(i)}|X_{k-1}^{(i)},Y_{k}}})={\mathcal {N}}({\boldsymbol {{\bar {x}}_{k}^{(i)},P_{k}^{(i)}}}),~~i=1,...,N
بنابراین گام نمونه‌برداری بااهمیت از الگوریتم فیلترذره‌ای در جدول زیر خلاصه شد.
به الگوریتم فیلتر ذره‌ای که توزیع پیشنهادی آن توسط فیلتر کالمن تخمین زده شود، فیلتر ذره‌ای کالمن (Kalman Particle Filter) یا به اختصار KPF می‌گویند.
تقریب مساحت نامعلوم  A_{s} (ناحیه ربع دایره) محصور (پَروسته) در مساحت معلوم مستطیلی شکل به کمک برآورد مونت کارلو
فرایند بازنمونه‌برداری