אלגוריתם מיקסום התוחלת

דוגמה לשימוש באלגוריתם EM כדי ליצור אשכולות (Clustering). בדוגמה זו מוצג מידע על התפרצויות גייזר אולד פיית'פול. ציר ה-X מבטא את משך ההתפרצות של הגייזר, וציר ה-Y מבטא את הזמן שחלף מההתפרצות הקודמת. האלגוריתם מגיע להתכנסות ולסיווג לשני סוגים של התפרצויות הגייזר אולד פיית'פול.

השוואה של k-מרכזים ומקסום תוחלת על נתוני הדמיה. אלגוריתם מקסום תוחלת מגיע לתיאור יותר נכון של ההתפלגויות הנורמליות, בעוד k-מרכזים מפצל את הנתונים לתאי וורונוי . מרכז האשכול מסומן על ידי הסמל הבהיר והגדול יותר.

אלגוריתם מיקסום התוחלת (באנגלית: Expectation–maximization; ובראשי תיבות: EM) הוא שיטה איטרטיבית לאמידה של פרמטרים, במודלים סטטיסטיים שבהם משולבים משתנים מקריים סמויים שאין עבורם תצפיות. מודל איטרציה זה עובר בין שני שלבים:

1. התוחלת (שלב ה-E), אשר יוצר פונקציה של הפרמטרים לאמידה - תוחלת של לוג הנראות (של כל הנתונים, כולל הסמויים) כאשר מניחים כי למשתנים הסמויים יש התפלגות שנסמכת על האומדן הנוכחי עבור הפרמטרים.
2. שלב המיקסום (שלב ה-M), אשר מוצא עבור הפרמטרים את האומדנים אשר מביאים למקסימום את הפונקציה שחושבה בשלב התוחלת.

לאחר שלב המיקסום חוזרים חוזרים לשלב התוחלת עם האומדנים שנמצאו בשלב המיקסום וכן הלאה.

היסטוריה

אלגוריתם EM קיבל את שמו לראשונה בשנת 1977 במאמר מאת ארתור דמפסטר, נאן ליירד, ודונלד רובין. במאמר זה ציינו החוקרים כי בעבר השיטה כבר הועלתה על הכתב בידי מחברים אחרים. דמפסטר-ליירד-רובין הגדירו את שיטת מיקסום התוחלת ככלי חשוב לפתרון מגוון רחב של בעיות. עם זאת, ניתוח ההתכנסות של דמפסטר-ליירד-רובין היה שגוי. ניתוח מתוקן פורסם על ידי ס.פ ג'ף וו ב-1983.

מבוא

אלגוריתם EM משמש תחליף לאומד נראות מקסימלית במודלים סטטיסטיים שבהם קיים קושי הנובע מקיומם של משתנים מקריים סמויים או נתונים חסרים. במודלים אלו קשה מבחינה טכנית למצוא מקסימום לפונקציית הנראות.

אלגוריתם EM מתחיל עם אומדן התחלתי לפרמטרים (אפשר לקבלו באמצעות שיטות אמידה פחות מדויקות). מתוך האומדן הראשוני מקבלים את ההתפלגות עבור הנתונים החסרים ואת $Q$ תוחלת ההתפלגות של לוג הנראות עבור נתונים מלאים. זו עדיין פונקציה של הפרמטרים, כי האומדן משמש אותנו רק כדי לקבל את ההתפלגות לנתונים החסרים. עתה ניתן להשתמש ב-$Q$ במקום בפונקציית לוג הנראות, $l$, שבה משתמשים בשיטה של אומד נראות מקסימלי. האומדן החדש לפרמטרים יהיה זה המביא למקסימום $Q$. בניגוד למקרה של אומד נראות מקסימלי, האומדן שקיבלנו עדיין לא אופטימלי ועלינו להציבו שוב כדי לקבל התפלגות משופרת לנתונים החסרים. כך ממשיכים עם האיטרציות עד אשר מגיעים למצב שבו איטרציות נוספות אינן משנות את האומדנים באופן שהוא משמעותי.

תיאור מתמטי

בהינתן מודל סטטיסטי המורכב ממדידות $x$, משתנים סמויים שאינם נראים או נתונים חסרים $Z$ (משתנה מקרי), פרמטרים θ אותם מעוניינים לאמוד ופונקציית הצפיפות $f(\mathbf {x} ,\mathbf {z} ;{\boldsymbol {\theta }})$ בהינתן הנתונים המלאים. פונקציית הצפיפות השולית של הנתונים הקיימים היא $f(\mathbf {x} ;{\boldsymbol {\theta }})=\int f(\mathbf {x} ,\mathbf {z} ;{\boldsymbol {\theta }})d\mathbf {z}$ ולכן לוג הנראות בהינתן הנתונים הנראים היא $l(\theta )=log(f(\mathbf {x} ;\theta ))$. כאשר אומדים בשיטת הנראות המקסימלית, בדרך כלל, מוצאים $\theta \in \Omega$ שיביא למקסימום את $l(\theta )$. אולם בדרך כלל במקרים של נתונים חסרים בדרך כלל מדובר בבעיה טכנית קשה יחסית.

בתנאים מסוימים, אלגוריתם EM מאפשר לעשות זאת על ידי יישום איטרטיבי של שני הצעדים הבאים:

שלב התוחלת (שלב ה-$E$): חישוב תוחלת הלוג של פונקציית הנראות לנתונים המלאים, ביחס להתפלגות של $Z$ של הנתונים החסרים, בהינתן סדרת התצפיות $\mathbf {x}$ והאומדן הנוכחי של הפרמטרים ${\boldsymbol {\theta }}^{(t)}$:

$Q\left({\boldsymbol {\theta }};{\boldsymbol {\theta }}^{(t)}\right)=\int \log(f(\mathbf {x} \mid \mathbf {z} ;{\boldsymbol {\theta }}))f\left(\mathbf {z} ;{\boldsymbol {\theta }}^{(t)}\right)\ d\mathbf {z} =E_{Z|\theta ^{(t)}}[\log f(\mathbf {x} \mid \mathbf {Z} ;{\boldsymbol {\theta }})]$

שלב המיקסום (שלב ה-$M$): מציאת ה-${\boldsymbol {\theta }}$ אשר מביא למקסימום את הפונקציה שחושבה בצעד הראשון:

${\boldsymbol {\theta }}^{(t+1)}={\underset {{\boldsymbol {\theta }}\in \Omega }{\operatorname {arg\,max} }}\ Q\left({\boldsymbol {\theta }};{\boldsymbol {\theta }}^{(t)}\right)\,$.

מתחילים עם אומדן ראשוני, $\mathbf {\theta ^{(0)}}$ (שיכול להתקבל באמצעות שיטה אחרת, פחות מדויקת), ואז חוזרים שוב ושוב על שלב התוחלת ושלב המיקסום ומקבלים סדרה $\mathbf {\theta ^{(0)}} ,\mathbf {\theta ^{(1)}} ,...$ של אומדנים שבתנאים מסוימים הולכים ומשתפרים.

משפט התכנסות של האלגוריתם

ג'ף וו הוכיח את משפט ההתכנסות הבא:

בהינתן שההנחות הבאות מתקיימות:

• $\Omega$ הוא תת-קבוצה של המרחב האוקלידי ה-$r$ ממדי $\Omega \subseteq \mathbb {R} ^{r}$,
• הקבוצה $\Omega _{\theta ^{(0)}}=\left\{\theta :l(\theta )\geq l\left(\theta ^{(0)}\right)\right\}$ קומפקטית לכל $l\left(\theta ^{(0)}\right)>-\infty$,
• $l$ רציפה על $\Omega$ וגזירה בפנים של $\Omega$,
• הפונקציה $Q({\boldsymbol {\theta }}';{\boldsymbol {\theta }})$ רציפה גם ב-${\boldsymbol {\theta }}$ וגם ב-${\boldsymbol {\theta }}'$;

אז כל נקודות הגבול של סדרה $\mathbf {\theta ^{(0)}} ,\mathbf {\theta ^{(1)}} ,...$ הן נקודות קריטיות של $l$ (נקודה בה הנגזרת מתאפסת או לא מוגדרת) והסדרה $l\left(\theta ^{(0)}\right),l\left(\theta ^{(1)}\right),...$ מתכנסת באופן מונוטוני ל-$l\left(\theta ^{*}\right)$ עבור נקודה $\theta ^{*}$, שהיא אחת מהנקודות הקריטיות של $l$.

אם נניח בנוסף לכך את ההנחות הבאות:

• לפונקציה $l$ יש נקודה קריטית יחידה $\theta ^{*}$ שבה מתקבל מקסימום גלובלי,
• וקטור הנגזרות החלקיות של $Q({\boldsymbol {\theta }}';{\boldsymbol {\theta }})$ רציף גם ב-${\boldsymbol {\theta }}$ וגם ב-${\boldsymbol {\theta }}'$;

אז הסדרה $\theta ^{(0)},\theta ^{(1)},...$ מתכנסת ל-$\theta ^{*}$.

דוגמה

התפלגות תערובת של שתי התפלגויות רב-נורמליות

נתון מדגם $\mathbf {x} =(\mathbf {x} _{1},\mathbf {x} _{2},\ldots ,\mathbf {x} _{n})$ של $n$ ווקטורי תצפיות בלתי תלויים מתערובת של שתי התפלגויות רב-נורמליות ממימד $d$, ונניח שקיימים משתנים מקריים סמויים $\mathbf {Z} =(Z_{1},Z_{2},\ldots ,Z_{n})$ אשר קובעים מאיזה מרכיב של התערובת מגיעה כל תצפית.

$X_{i}\mid Z_{i}=0\sim {\mathcal {N}}_{d}({\boldsymbol {\mu }}_{0},\Sigma _{0})$ ו-

$X_{i}\mid Z_{i}=1\sim {\mathcal {N}}_{d}({\boldsymbol {\mu }}_{1},\Sigma _{1})$,

כאשר

$\operatorname {P} (Z_{i}=1)=\tau \,$ ו-

$\operatorname {P} (Z_{i}=0)=1-\tau$.

המטרה היא לאמוד את המשקלים של הצירוף הקמור בין שתי ההתפלגויות המרכיבות ואת התוחלות והשונויות של כל אחד מהם:

$\theta ={\big (}\tau ,{\boldsymbol {\mu }}_{0},{\boldsymbol {\mu }}_{1},\Sigma _{0},\Sigma _{1}{\big )}$,

פונקציית הנראות עבור הנתונים המלאים הכוללים וקטור תצפיות $\mathbf {z}$:

$f(\mathbf {x} ,\mathbf {z} \mid \theta )=\prod _{i=1}^{n}\ [\left((1-\tau )f(\mathbf {x} _{i};{\boldsymbol {\mu }}_{0},\Sigma _{0})\right)^{1-z_{i}}\left(\tau f(\mathbf {x} _{i};{\boldsymbol {\mu }}_{1},\Sigma _{1})\right)^{z_{i}}]$,

כאשר $f(\mathbf {x} _{i};{\boldsymbol {\mu }}_{j},\Sigma _{j})$ הן פונקציות צפיפות של של התפלגויות רב-נורמליות.

הלוג של פונקציית הנראות הוא:

$\log f(\mathbf {x} ,\mathbf {z} ;\theta )=\sum _{i=1}^{n}\ [(1-z_{i})\left(\log(1-\tau )+\log f(\mathbf {x} _{i};{\boldsymbol {\mu }}_{0},\Sigma _{0})\right)+z_{i}\left(\log \tau +\log f(\mathbf {x} _{i};{\boldsymbol {\mu }}_{1},\Sigma _{1})\right)]$.

מאחר ש-$\mathbf {z}$ מבטא נתונים חסרים נציב במקומו משתנים מקריים:

$\log f(\mathbf {x} ,\mathbf {Z} ;\theta )=\sum _{i=1}^{n}\ [(1-Z_{i})\left(\log(1-\tau )+\log f(\mathbf {x} _{i};{\boldsymbol {\mu }}_{0},\Sigma _{0})\right)+Z_{i}\left(\log \tau +\log f(\mathbf {x} _{i};{\boldsymbol {\mu }}_{1},\Sigma _{1})\right)]$.

שלב התוחלת

בהינתן האומדן $\theta ^{(t)}$, ההתפלגות המותנה של $Z_{i}$ נקבעת על ידי חוק בייס:

$\tau _{i,1}^{(t)}:=\operatorname {P} (Z_{i}=1\mid \mathbf {x} _{i};\theta ^{(t)})={\frac {\tau ^{(t)}\ f(\mathbf {x} _{i};{\boldsymbol {\mu }}_{1}^{(t)},\Sigma _{1}^{(t)})}{(1-\tau ^{(t)})\ f(\mathbf {x} _{i};{\boldsymbol {\mu }}_{0}^{(t)},\Sigma _{0}^{(t)})+\tau ^{(t)}\ f(\mathbf {x} _{i};{\boldsymbol {\mu }}_{1}^{(t)},\Sigma _{1}^{(t)})}}$.

לכן,

$\operatorname {E} _{Z|\theta ^{(t)}}[Z_{i}=1\mid X_{i}=\mathbf {x} _{i};\theta ^{(t)}]=\tau _{i,1}^{(t)}$.

נסמן גם $\tau _{i,0}^{(t)}=1-\tau _{i,1}^{(t)}$. לכן,

$Q(\theta \mid \theta ^{(t)})=\operatorname {E} _{Z|\theta ^{(t)}}[\log f(\mathbf {x} ,\mathbf {Z} ;\theta )]=\sum _{i=1}^{n}\ [\tau _{i,0}^{(t)}\left(\log(1-\tau )+\log f(\mathbf {x} _{i};{\boldsymbol {\mu }}_{0},\Sigma _{0})\right)+\tau _{i,1}^{(t)}\left(\log \tau +\log f(\mathbf {x} _{i};{\boldsymbol {\mu }}_{1},\Sigma _{1})\right)]$

שלב המיקסום

בשלב זה מקבלים אומדנים חדשים לפרמטרים, כאלה אשר מביאים את $Q$ למקסימום. נעשה זאת, במקרה המסוים הזה, על ידי גזירה(חלקית) לפי הפרמטרים והשוואה ל-0,

$.{\boldsymbol {\tau }}^{(t+1)}={\underset {\boldsymbol {\tau }}{\operatorname {arg\,max} }}\ Q(\theta \mid \theta ^{(t)})={\underset {\boldsymbol {\tau }}{\operatorname {arg\,max} }}\ \left\{\left[\sum _{i=1}^{n}\tau _{i,0}^{(t)}\right]\log(1-\tau )+\left[\sum _{i=1}^{n}\tau _{i,1}^{(t)}\right]\log \tau \right\}$

מקבלים,

$\tau ^{(t+1)}={\frac {1}{n}}\sum _{i=1}^{n}\tau _{i,1}^{(t)}$.

עבור $j=0,1$, מציאת האומדנים החדשים של ${\boldsymbol {\mu }}_{j}$ ו- ${\boldsymbol {\Sigma }}_{j}$,

$.({\boldsymbol {\mu }}_{j}^{(t+1)},{\boldsymbol {\Sigma }}_{j}^{(t+1)})={\underset {{\boldsymbol {\mu }}_{j},\Sigma _{j}}{\operatorname {arg\,max} }}\ Q(\theta \mid \theta ^{(t)})={\underset {{\boldsymbol {\mu }}_{j},\Sigma _{j}}{\operatorname {arg\,max} }}\ \sum _{i=1}^{n}\tau _{i,j}^{(t)}\left\{-{\tfrac {1}{2}}\log |\Sigma _{j}|-{\tfrac {1}{2}}(\mathbf {x} _{i}-{\boldsymbol {\mu }}_{j})^{\top }\Sigma _{j}^{-1}(\mathbf {x} _{i}-{\boldsymbol {\mu }}_{j})\right\}$

מקבלים,

${\boldsymbol {\mu }}_{j}^{(t+1)}=\left(\sum _{i=1}^{n}\tau _{i,j}^{(t)}\right)^{-1}\sum _{i=1}^{n}\tau _{i,j}^{(t)}\mathbf {x} _{i}$     וגם     ${\boldsymbol {\Sigma }}_{j}^{(t+1)}=\left(\sum _{i=1}^{n}\tau _{i,j}^{(t)}\right)^{-1}\sum _{i=1}^{n}\tau _{i,j}^{(t)}(\mathbf {x} _{i}-{\boldsymbol {\mu }}_{j}^{(t+1)})(\mathbf {x} _{i}-{\boldsymbol {\mu }}_{j}^{(t+1)})^{\top }$.

הכללה

ניתן להכליל את האלגוריתם המודגם לעיל עבור תערובות של יותר משתי התפלגויות רב-נורמליות.