Skryté Markovove modely: Viterbiho algoritmus a nepoctivé kasíno

Nech N, M ≥ 1 sú celé čísla, S = {1, …, N} množina skrytých stavov a V = {1, …, M} množina pozorovateľných symbolov. Ďalej nech

• π = (π₁, …, π_N) ∈ ℝ^1×N je rozdelenie pravdepodobnosti na S (π_i ≥ 0, ∑_i π_i = 1) — počiatočné rozdelenie;
• A = (a_ij) ∈ ℝ^N×N je stochastická matica — prechodová matica;
• B = (b_i(o))_{i ∈ S, o ∈ V} ∈ ℝ^N×M je stochastická matica — emisné pravdepodobnosti.

Skrytý Markovov model (HMM) s parametrami λ = (π, A, B) je dvojica postupností náhodných premenných X₁, X₂, … : Ω → S (skryté stavy) a O₁, O₂, … : Ω → V (pozorovania) na pravdepodobnostnom priestore (Ω, 𝓕, P), pre ktorú platí:

1. (X_t)_{t ≥ 1} je homogénny Markovov reťazec s počiatočným rozdelením P(X₁ = i) = π_i a prechodovou maticou A, t. j. a_ij = P(X_t+1 = j | X_t = i);
2. pre každé T ≥ 1, t ∈ {1, …, T}, o ∈ V a všetky hodnoty x₁, …, x_T ∈ S a o_s ∈ V (s ≠ t), pre ktoré má podmienka kladnú pravdepodobnosť, P(O_t = o | X₁ = x₁, …, X_T = x_T, O_s = o_s pre s ∈ {1, …, T} ∖ {t}) = b_xt(o).

Slovami: skryté stavy sa menia ako obyčajný Markovov reťazec a pozorovanie v čase t závisí len od stavu X_t — nie od ostatných stavov ani od ostatných pozorovaní. Pozorujeme len O₁, O₂, …, stavy X_t nie.

Náš rad začína hodmi 1, 6, 3. Označme P poctivú a F falošnú kocku. História „trikrát falošná“:

P(FFF, 1, 6, 3) = 0,5 · 0,1 · 0,9 · 0,5 · 0,9 · 0,1 = 0,002025,

história „trikrát poctivá“:

P(PPP, 1, 6, 3) = 0,5 · 1/6 · 0,95 · 1/6 · 0,95 · 1/6 = 0,002089

(farby podľa vzorca vyššie: počiatočné rozdelenie, emisie, prechody). Poctivá história vyjde o chlp pravdepodobnejšia: šestka na falošnej kocke síce dostane 0,5 namiesto 1/6, ale jednotka a trojka len 0,1 namiesto 1/6. (V simulácii sa pritom všetky tri hody hádzali falošnou kockou — z troch hodov sa to spoznať nedá.)

Pravdepodobnosť pozorovania je súčet združenej pravdepodobnosti cez všetky možné histórie: P(O | λ) = ∑_x P(x, O | λ). Pre hody 1, 6, 3 ich je 2³ = 8:

Každá história stojí 2T − 1 = 5 násobení, dohromady 8 · 5 = 40. Pre T hodov je to (2T − 1) · 2^T — pri 300 hodoch nepredstaviteľne veľa.

Vstup: skrytý Markovov model λ = (π, A, B) so stavmi S = {1, …, N} a symbolmi V; pozorovaná postupnosť O = (O₁, …, O_T) ∈ V^T, T ≥ 1.
Výstup: P(O | λ) = P(O₁, …, O_T) a dopredné premenné

α_t(j) = P(O₁, …, O_t, X_t = j | λ),   t = 1, …, T,   j ∈ S.

1. Inicializácia: pre každé j ∈ S polož α₁(j) = π_j b_j(O₁).
2. Indukcia: pre t = 1, …, T − 1 a každé j ∈ S polož α_t+1(j) = [∑_i=1^N α_t(i) a_ij] · b_j(O_t+1).
3. Ukončenie: po spracovaní času T vráť P(O | λ) = ∑_j=1^N α_T(j).

Algoritmus skončí po T − 1 krokoch indukcie; každý stojí O(N²) operácií, spolu teda O(N² T) namiesto O(T · N^T) pri sčítaní cez všetky histórie. Aby čísla nepodtiekli, v každom čase t sa vypočíta c_t = ∑_j α_t(j) (s už normalizovaným α_t−1 v indukcii), α_t(j) sa nahradí α_t(j)/c_t a vracia sa log P(O | λ) = ∑_t=1^T log c_t.

Štart: α₁(P) = 0,5 · 1/6 = 0,0833, α₁(F) = 0,5 · 0,1 = 0,0500. Druhý hod je šestka (farby podľa vzorca vyššie):

α₂(P) = (0,0833 · 0,95 + 0,0500 · 0,1) · 1/6 = 0,0140,

α₂(F) = (0,0833 · 0,05 + 0,0500 · 0,9) · 0,5 = 0,0246.

Súčet v poslednom riadku je presne súčet ôsmich histórií z predchádzajúceho príkladu, len namiesto 40 násobení stačilo 2 + 2 · 6 = 14. Posledný stĺpec c_t = P(O_t | O₁, …, O_t−1) je normalizačná konštanta zo škálovania: aký pravdepodobný bol hod t vzhľadom na tie predchádzajúce. Ich súčin je P(O | λ), takže log P(O | λ) = ∑_t log c_t.

Prečo logaritmy? Pri 3 000 hodoch vyjde P(O | λ) okolo 10^−2 278, a najmenšie kladné číslo, ktoré počítač v bežnej dvojitej presnosti (IEEE 754) uloží, je asi 4,9 · 10⁻³²⁴. Bez škálovania by α podtieklo na nulu.

Nech (Ω, 𝓕, P) je pravdepodobnostný priestor, M₁, M₂ ∈ 𝓕 sú javy „platí model 1“ a „platí model 2“ s P(M₁) > 0, P(M₂) > 0 a O ∈ 𝓕 je jav „pozorovali sme dáta“ s P(O | M₂) > 0. Bayesov faktor modelu M₁ proti M₂ je číslo

K = P(O | M₁)/P(O | M₂),

t. j. podiel vierohodností dát v oboch modeloch. K > 1 znamená, že dáta sú pravdepodobnejšie v modeli M₁. Ak je navyše P(O) > 0, vyplýva z Bayesovej vety napísanej pre oba modely a vydelenej (P(O) sa skráti):

Výpočet. log₁₀ K = −230,03 − (−233,45) = 3,42 (z nezaokrúhlených hodnôt 3,413), teda K ≈ 10^3,413 ≈ 2 590, zaokrúhlene 2 600. Kto vopred považoval podvod a poctivú hru za rovnako pravdepodobné (pomer šancí 1), má po 300 hodoch šance 2 600 · 1 = 2 600 : 1 v prospech podvodu. Podľa bežne používanej stupnice Kassa a Rafteryho (1995) je Bayesov faktor nad 150 „veľmi silný“ dôkaz.

Vstup: skrytý Markovov model λ = (π, A, B) so stavmi S = {1, …, N}; pozorovanie O = (O₁, …, O_T) ∈ V^T, T ≥ 1.
Výstup: najpravdepodobnejšia postupnosť skrytých stavov a jej pravdepodobnosť

x* = argmax_{x ∈ S^T} P(X₁ = x₁, …, X_T = x_T, O | λ),   P* = P(x*, O | λ).

Pre t = 1, …, T a j ∈ S označme δ_t(j) pravdepodobnosť najlepšej cesty, ktorá v čase t končí v stave j:

δ_t(j) = max_{x1, …, xt−1 ∈ S} P(X₁ = x₁, …, X_t−1 = x_t−1, X_t = j, O₁, …, O_t | λ).

1. Inicializácia: pre každé j ∈ S polož δ₁(j) = π_j b_j(O₁).
2. Indukcia: pre t = 1, …, T − 1 a každé j ∈ S polož δ_t+1(j) = max_{i ∈ S} [δ_t(i) a_ij] · b_j(O_t+1),   ψ_t+1(j) = argmax_{i ∈ S} δ_t(i) a_ij (ψ_t+1(j) je najlepší predchodca; pri zhode vyber ktoréhokoľvek).
3. Ukončenie: P* = max_{j ∈ S} δ_T(j) a x*_T = argmax_{j ∈ S} δ_T(j).
4. Spätný prechod: pre t = T − 1, T − 2, …, 1 polož x*_t = ψ_t+1(x*_t+1); po t = 1 vráť x* = (x*₁, …, x*_T) a P*.

Algoritmus skončí po T − 1 krokoch indukcie a T − 1 krokoch spätného prechodu; čas je opäť O(N² T). V praxi sa namiesto súčinov sčítavajú logaritmy, aby čísla nepodtiekli.

Štart je rovnaký ako pri α: δ₁(P) = 0,0833, δ₁(F) = 0,0500. Pre šestku v druhom hode:

δ₂(P) = max(0,0833 · 0,95; 0,0500 · 0,1) · 1/6 = 0,0792 · 1/6 = 0,0132   (vyhral P),

δ₂(F) = max(0,0833 · 0,05; 0,0500 · 0,9) · 0,5 = 0,0450 · 0,5 = 0,0225   (vyhral F)

(farby podľa vzorca vyššie).

Na konci vyhrá P (0,002089 > 0,002025) a spätne po odkazoch ψ vyjde cesta PPP — presne história z prvého príkladu, ktorá bola zo všetkých ôsmich najpravdepodobnejšia. Dopredný algoritmus však pre tretí hod dáva P(F | O₁, O₂, O₃) = α₃(F)/α₃(P) + α₃(F) = 0,002283/0,004913 = 0,465. Najlepšia cesta a najpravdepodobnejší stav sa tu zhodujú (P), ale „víťazstvo“ je tesné a tri hody na istotu nestačia. Vo všeobecnosti sa zhodovať nemusia — k tomu sa vrátime pri aposteriórnom dekódovaní.

Nech λ = (π, A, B) je skrytý Markovov model so stavmi S = {1, …, N} a O = (O₁, …, O_T) pozorovanie s P(O | λ) > 0. Aposteriórna pravdepodobnosť stavu i ∈ S v čase t ∈ {1, …, T} je

γ_t(i) = P(X_t = i | O₁, …, O_T, λ).

Na jej výpočet slúžia spätné premenné (pre t < T a P(X_t = i) > 0)

β_t(i) = P(O_t+1, …, O_T | X_t = i, λ),   β_T(i) = 1.

Dopredno-spätný algoritmus vypočíta α_t(i) dopredným algoritmom, β_t(i) odzadu (pre t = T − 1, …, 1) a z oboch γ_t(i) podľa vzorcov nižšie. Slovami: γ_t(i) je pravdepodobnosť, že v čase t bol model v stave i, keď poznáme celý pozorovaný rad — minulosť aj budúcnosť.

Zostaňme pri hodoch 1, 6, 3 a pýtajme sa na druhý hod (šestku). Spätný krok od β₃ = 1 s tretím hodom (trojkou):

β₂(P) = 0,95 · 1/6 + 0,05 · 0,1 = 0,1633,   β₂(F) = 0,1 · 1/6 + 0,9 · 0,1 = 0,1067.

(V každom súčine je prvý činiteľ prechod a_ij, druhý emisia trojky b_j(3) a β₃(j) = 1 sa nepíše.) Podľa vzorca vyššie (farby zodpovedajú jeho členom):

γ₂(F) = 0,024583 · 0,106667/0,004913 = 0,534.

Keby sme poznali len prvé dva hody, bolo by P(F | O₁, O₂) = α₂(F)/α₂(P) + α₂(F) = 0,637. Trojka, ktorá padla potom, je na falošnej kocke menej pravdepodobná, a tak podozrenie na druhý hod spätne zníži na 0,534.

Nech λ je skrytý Markovov model so stavmi S, O = (O₁, …, O_T) pozorovanie a γ_t(i) aposteriórne pravdepodobnosti stavov z predchádzajúcej definície. Aposteriórne dekódovanie je postupnosť x̂ = (x̂₁, …, x̂_T) ∈ S^T, kde

x̂_t = argmax_{i ∈ S} γ_t(i)   pre každé t = 1, …, T

(stav sa volí v každom čase zvlášť). Nech ďalej x ∈ S^T je skutočná postupnosť stavov a y ∈ S^T ľubovoľné dekódovanie (aposteriórne alebo Viterbiho). Presnosť dekódovania y je podiel časov, v ktorých sa zhoduje so skutočnosťou:

presnosť(y, x) = |{t ∈ {1, …, T} : y_t = x_t}|/T ∈ [0, 1].

Aposteriórne dekódovanie maximalizuje očakávaný počet správne určených hodov, Viterbi pravdepodobnosť, že je správna celá cesta naraz. Preto môže aposteriórna postupnosť obsahovať aj prechod, ktorý má v modeli pravdepodobnosť nula — pri kasíne taký nie je, ale vo všeobecnosti to nastať môže.

Je to algoritmus EM (Baum a kol., okolo roku 1970).

Vstup: počet skrytých stavov N ≥ 1, množina symbolov V, pozorovanie O = (O₁, …, O_T) ∈ V^T s T ≥ 2, počiatočný odhad parametrov λ⁽⁰⁾ = (π, A, B) s P(O | λ⁽⁰⁾) > 0 a počet iterácií K ≥ 1 (alebo tolerancia ε > 0).
Výstup: parametre λ^(K) s vierohodnosťou P(O | λ^(K)) ≥ P(O | λ⁽⁰⁾).

Pre k = 0, 1, …, K − 1 opakuj s parametrami λ^(k):

1. Krok E: dopredno-spätným algoritmom vypočítaj α_t(i), β_t(i), P(O | λ^(k)), aposteriórne pravdepodobnosti γ_t(i) pre t = 1, …, T a pravdepodobnosti prechodov ξ_t(i, j) = P(X_t = i, X_t+1 = j | O, λ^(k)),   t = 1, …, T − 1,   i, j ∈ S.
2. Krok M: nové parametre λ^(k+1) sú pre všetky i, j ∈ S a o ∈ V π_i = γ₁(i),   a_ij = ∑_t=1^T−1 ξ_t(i, j)/∑_t=1^T−1 γ_t(i),   b_i(o) = ∑_{t : Ot = o} γ_t(i)/∑_t=1^T γ_t(i).
3. Zastavenie: skonči po K iteráciách, prípadne skôr, keď log P(O | λ^(k+1)) − log P(O | λ^(k)) < ε.

Platí P(O | λ^(k+1)) ≥ P(O | λ^(k)) — vierohodnosť v žiadnej iterácii neklesne. Algoritmus konverguje k stacionárnemu bodu vierohodnosti, v praxi obvykle k lokálnemu maximu, ktoré nemusí byť globálne; výsledok preto závisí od počiatočného odhadu.