$\textbf{Expectation-Maximization algoritmus}\\$
Gyakorlatban sokszor az eloszlás bonyolult formája miatt nehéz ML-becslést adni egy rejtett Markov-modellhez tartozó eloszlás paraméterére, ezért azt közelítenünk kell. Egy lehetséges közelítés megadására szolgál a következő algoritmus:
Legyen adott egy (diszkrét idejű) $(X_n,Y_n)_{(1 \leq n \leq m)}$ rejtett Markov-modell, a megfigyelések $y=(y_1,\dots,y_m)$, azokhoz tartozó rejtett állapotok $x=(x_1,\dots,x_m)$, az eloszlás
$$(X,Y) \sim p((x,y) \mid \theta) = p_{\theta}(x,y)=\exp(\eta (\theta) \cdot s(x,y)) \cdot h(x,y)/c(\theta),$$
ahol $\eta$ monoton függvény, $\Theta \subseteq \mathbb{R}^n$ nyílt halmaz a vizsgált eloszlástípus lehetséges paramétereinek halmaza,
$$s: \mathbb{R}^d \rightarrow \mathbb{R}^k$$
elégséges statisztika,
$$h: \mathbb{R}^d \rightarrow \mathbb{R}^{\geq 0},$$
$$c: \Theta \rightarrow \mathbb{R}^{\geq 0}$$
normalizáló függvény. Tehát
$$c(\theta)= \int_{(X,Y)} \left(\exp( \eta (\theta) \cdot s(x,y))\cdot h(x,y)\right)d(x,y)$$
Célunk egy $\theta_{MLE} \in \Theta$ megadása, melyre
$$\theta_{MLE} \in \arg\max_{\theta}(p_{\theta}(y)),$$
vagyis maximalizálja a
$$p(\theta) = p_{\theta}(y)= \int p_{\theta}(x,y) dx$$
peremeloszlás likelihood-függvényét (legalábbis elég jól közelíti egy globális maximumhelyét). Feltehető tehát, hogy $\eta (\theta) = \theta$.
$\\ \textbf{Állítás}:\\$
$$\nabla \log(c(\theta))=\mathbb{E}_{\theta}s(x,y)$$
$\\ \textbf{Bizonyítás}:\\$
$\forall j$-re
$$(\partial/\partial \theta_j) \log(c(\theta))=\frac{1}{c(\theta)} \cdot (\partial/\partial \theta_j) c(\theta)$$
$$= \frac{1}{c(\theta)} \int s_j(x,y) p_{\theta}(x,y)d(x,y)= \mathbb{E}_{\theta}s_j(x,y)$$
(feltettük hogy a deriválás és az integrálás felcserélhető)
$\\ \textbf{Következmény}:\\$
$\theta$ a likelihood függvény szélsőértéke
$$\Rightarrow \mathbb{E}_{\theta}(s(X,Y) \mid Y=y)=\mathbb{E}_{\theta}(s(X,Y))$$
$\\ \textbf{Bizonyítás}:\\$
$$0 = (\partial/\partial \theta_i) \log(p_{\theta}(y)) = \frac{1}{p(\theta)} \cdot \int (\partial/\partial \theta_i) p_{\theta}(x,y)dx$$
$$=\frac{1}{p(\theta)} \cdot \int (\partial/\partial \theta_i) \exp(\theta s(x,y) - \log(c(\theta))) \cdot h(x,y) dx$$
$$=\frac{1}{p(\theta)} \cdot \int (s_i (x,y) - (\partial/\partial \theta_i)\log(c(\theta))) \cdot p_{\theta}(x,y)dx$$
$$= \frac{1}{p(\theta)} \cdot \int (s_i (x,y) - \mathbb{E}_{\theta} s_i(X,Y)) \cdot p_{\theta}(x,y)dx$$
$$= \int s_i (x,y) \cdot p_{\theta}(x \mid y) dx - \mathbb{E}_{\theta} s_i(X,Y)$$
$$= \mathbb{E}_{\theta}(s_i(X,Y) \mid Y=y)- \mathbb{E}_{\theta} s_i(X,Y) \iff \mathbb{E}_{\theta}(s(X,Y) \mid Y=y)=\mathbb{E}_{\theta}(s(X,Y))$$
Legyen $\forall t$-re
$$Q( \theta \mid \theta_t)= \mathbb{E}_{\theta_t}(\log(p_{\theta}(X,Y)) \mid Y=y)$$
Az algoritmus ötlete, hogy kiindulunk tetszőleges $\theta_0 \in \Theta$ becslésből és minden iterációban javítunk rajta úgy, hogy megoldjuk a
$$\theta_{t+1} = \arg\max_{\theta} Q(\theta \mid \theta_{t})$$
rekurziót. Vegyük észre, hogy az előző állítás miatt
$$\nabla Q(\theta \mid \theta_{t}) = \nabla ( \theta \cdot \mathbb{E}_{\theta_t}(s(X,Y) \mid Y=y) - \log(c(\theta)) - \text{konstans})$$
$$= \mathbb{E}_{\theta_t}(s(X,Y) \mid Y=y)- \mathbb{E}_{\theta}(s(X,Y))$$
Az eredeti rekurzió megoldásához tehát megoldandó az
$$\mathbb{E}_{\theta_{t-1}}(s(X,Y) \mid Y=y)=\mathbb{E}_{\theta_t}(s(X,Y))$$
rekurzió.
Konyhaszabály: addig iteráljuk, mígnem rögzített $t$-re $\mathbb{E}_{\theta_{t}}(s(X,Y) \mid Y=y)$ és $\mathbb{E}_{\theta_t}(s(X,Y))$ majdnem egyenlő (mondjuk az eltérésük $< 10^{-5}$), ekkor ugyanis a becslés előbbi következmény miatt közelít egy szélsőértékhez. A konvergencia abból adódik, hogy a $\theta_t$-khez tartozó értékekből alkotott sorozat monoton és korlátos, márpedig az
$$\mathbb{E}_{\theta_{t-1}}(s(X,Y) \mid Y=y)=\mathbb{E}_{\theta_t}(s(X,Y))$$
rekurzió miatt és mert
$$Q( \theta \mid \theta_t)= \mathbb{E}_{\theta_t}(\log(p_{\theta}(X,Y)) \mid Y=y)$$
sorozat konvergens, adódik, hogy $\exists t$: $\mathbb{E}_{\theta_{t}}(s(X,Y) \mid Y=y)$ és $\mathbb{E}_{\theta_t}(s(X,Y))$ majdnem egyenlő (mert a tagokból alkotott sorozatoknak ugyanaz a határértéke).
A kapott szélsőérték nem feltétlenül globális maximum, ennek érdekében különböző $\theta_0$ kezdőértékekből kell indítanunk az algoritmust. A legjobb becslést globális maximumnak feltételezzük.