Skryté Markovovy modely: Viterbiho algoritmus a nepoctivé kasino

Nechť N, M ≥ 1 jsou celá čísla, S = {1, …, N} množina skrytých stavů a V = {1, …, M} množina pozorovatelných symbolů. Dále nechť

• π = (π₁, …, π_N) ∈ ℝ^1×N je rozdělení pravděpodobnosti na S (π_i ≥ 0, ∑_i π_i = 1) — počáteční rozdělení;
• A = (a_ij) ∈ ℝ^N×N je stochastická matice — přechodová matice;
• B = (b_i(o))_{i ∈ S, o ∈ V} ∈ ℝ^N×M je stochastická matice — emisní pravděpodobnosti.

Skrytý Markovův model (HMM) s parametry λ = (π, A, B) je dvojice posloupností náhodných veličin X₁, X₂, … : Ω → S (skryté stavy) a O₁, O₂, … : Ω → V (pozorování) na pravděpodobnostním prostoru (Ω, 𝓕, P), pro kterou platí:

1. (X_t)_{t ≥ 1} je homogenní Markovův řetězec s počátečním rozdělením P(X₁ = i) = π_i a přechodovou maticí A, tj. a_ij = P(X_t+1 = j | X_t = i);
2. pro každé T ≥ 1, t ∈ {1, …, T}, o ∈ V a všechny hodnoty x₁, …, x_T ∈ S a o_s ∈ V (s ≠ t), pro které má podmínka kladnou pravděpodobnost, P(O_t = o | X₁ = x₁, …, X_T = x_T, O_s = o_s pro s ∈ {1, …, T} ∖ {t}) = b_xt(o).

Slovy: skryté stavy se mění jako obyčejný Markovův řetězec a pozorování v čase t závisí jen na stavu X_t — ne na ostatních stavech ani na ostatních pozorováních. Pozorujeme jen O₁, O₂, …, stavy X_t ne.

Naše řada začíná hody 1, 6, 3. Označme P poctivou a F falešnou kostku. Historie „třikrát falešná“:

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

historie „třikrát poctivá“:

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

(barvy podle vzorce výše: počáteční rozdělení, emise, přechody). Poctivá historie vyjde o chlup pravděpodobnější: šestka na falešné kostce sice dostane 0,5 místo 1/6, ale jednička a trojka jen 0,1 místo 1/6. (V simulaci se přitom všechny tři hody házely falešnou kostkou — ze tří hodů se to poznat nedá.)

Pravděpodobnost pozorování je součet sdružené pravděpodobnosti přes všechny možné historie: P(O | λ) = ∑_x P(x, O | λ). Pro hody 1, 6, 3 jich je 2³ = 8:

Každá historie stojí 2T − 1 = 5 násobení, dohromady 8 · 5 = 40. Pro T hodů je to (2T − 1) · 2^T — u 300 hodů nepředstavitelně mnoho.

Vstup: skrytý Markovův model λ = (π, A, B) se stavy S = {1, …, N} a symboly V; pozorovaná posloupnost O = (O₁, …, O_T) ∈ V^T, T ≥ 1.
Výstup: P(O | λ) = P(O₁, …, O_T) a dopředné proměnné

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

1. Inicializace: pro každé j ∈ S polož α₁(j) = π_j b_j(O₁).
2. Indukce: pro 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čení: po zpracování času T vrať P(O | λ) = ∑_j=1^N α_T(j).

Algoritmus skončí po T − 1 krocích indukce; každý stojí O(N²) operací, celkem tedy O(N² T) místo O(T · N^T) při sčítání přes všechny historie. Aby čísla nepodtekla, v každém čase t se spočítá c_t = ∑_j α_t(j) (s už normalizovaným α_t−1 v indukci), α_t(j) se nahradí α_t(j)/c_t a vrací se log P(O | λ) = ∑_t=1^T log c_t.

Start: α₁(P) = 0,5 · 1/6 = 0,0833, α₁(F) = 0,5 · 0,1 = 0,0500. Druhý hod je šestka (barvy podle vzorce výše):

α₂(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.

Součet v posledním řádku je přesně součet osmi historií z předchozího příkladu, jen místo 40 násobení stačilo 2 + 2 · 6 = 14. Poslední sloupec c_t = P(O_t | O₁, …, O_t−1) je normalizační konstanta ze škálování: jak pravděpodobný byl hod t vzhledem k těm předchozím. Jejich součin je P(O | λ), takže log P(O | λ) = ∑_t log c_t.

Proč logaritmy? U 3 000 hodů vyjde P(O | λ) kolem 10^−2 278, a nejmenší kladné číslo, které počítač v běžné dvojité přesnosti (IEEE 754) uloží, je asi 4,9 · 10⁻³²⁴. Bez škálování by α podteklo na nulu.

Nechť (Ω, 𝓕, P) je pravděpodobnostní prostor, M₁, M₂ ∈ 𝓕 jsou jevy „platí model 1“ a „platí model 2“ s P(M₁) > 0, P(M₂) > 0 a O ∈ 𝓕 je jev „pozorovali jsme data“ s P(O | M₂) > 0. Bayesův faktor modelu M₁ proti M₂ je číslo

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

tj. podíl věrohodností dat v obou modelech. K > 1 znamená, že data jsou pravděpodobnější v modelu M₁. Je-li navíc P(O) > 0, plyne z Bayesovy věty napsané pro oba modely a vydělené (P(O) se zkrátí):

Výpočet. log₁₀ K = −230,03 − (−233,45) = 3,42 (z nezaokrouhlených hodnot 3,413), tedy K ≈ 10^3,413 ≈ 2 590, zaokrouhleně 2 600. Kdo předem považoval podvod a poctivou hru za stejně pravděpodobné (poměr šancí 1), má po 300 hodech šance 2 600 · 1 = 2 600 : 1 pro podvod. Podle běžně používané stupnice Kasse a Rafteryho (1995) je Bayesův faktor nad 150 „velmi silný“ důkaz.

Vstup: skrytý Markovův model λ = (π, A, B) se stavy S = {1, …, N}; pozorování O = (O₁, …, O_T) ∈ V^T, T ≥ 1.
Výstup: nejpravděpodobnější posloupnost skrytých stavů a její pravděpodobnost

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

Pro t = 1, …, T a j ∈ S označme δ_t(j) pravděpodobnost nejlepší cesty, která v čase t končí ve stavu j:

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

1. Inicializace: pro každé j ∈ S polož δ₁(j) = π_j b_j(O₁).
2. Indukce: pro 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 nejlepší předchůdce; při shodě vyber kterýkoli).
3. Ukončení: P* = max_{j ∈ S} δ_T(j) a x*_T = argmax_{j ∈ S} δ_T(j).
4. Zpětný průchod: pro t = T − 1, T − 2, …, 1 polož x*_t = ψ_t+1(x*_t+1); po t = 1 vrať x* = (x*₁, …, x*_T) a P*.

Algoritmus skončí po T − 1 krocích indukce a T − 1 krocích zpětného průchodu; čas je opět O(N² T). V praxi se místo součinů sčítají logaritmy, aby čísla nepodtekla.

Start je stejný jako u α: δ₁(P) = 0,0833, δ₁(F) = 0,0500. Pro šestku ve druhém hodu:

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

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

(barvy podle vzorce výše).

Na konci vyhraje P (0,002089 > 0,002025) a zpětně po odkazech ψ vyjde cesta PPP — přesně historie z prvního příkladu, která byla ze všech osmi nejpravděpodobnější. Dopředný algoritmus ale pro třetí hod dává P(F | O₁, O₂, O₃) = α₃(F)/α₃(P) + α₃(F) = 0,002283/0,004913 = 0,465. Nejlepší cesta a nejpravděpodobnější stav se tu shodují (P), ale „vítězství“ je těsné a tři hody na jistotu nestačí. Obecně se shodovat nemusejí — k tomu se vrátíme u aposteriorního dekódování.

Nechť λ = (π, A, B) je skrytý Markovův model se stavy S = {1, …, N} a O = (O₁, …, O_T) pozorování s P(O | λ) > 0. Aposteriorní pravděpodobnost stavu i ∈ S v čase t ∈ {1, …, T} je

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

K jejímu výpočtu slouží zpětné proměnné (pro t < T a P(X_t = i) > 0)

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

Dopředně-zpětný algoritmus spočítá α_t(i) dopředným algoritmem, β_t(i) odzadu (pro t = T − 1, …, 1) a z obou γ_t(i) podle vzorců níže. Slovy: γ_t(i) je pravděpodobnost, že v čase t byl model ve stavu i, když známe celou pozorovanou řadu — minulost i budoucnost.

Zůstaňme u hodů 1, 6, 3 a ptejme se na druhý hod (šestku). Zpětný krok od β₃ = 1 s třetím hodem (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ždém součinu je první činitel přechod a_ij, druhý emise trojky b_j(3) a β₃(j) = 1 se nepíše.) Podle vzorce výše (barvy odpovídají jeho členům):

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

Kdybychom znali jen první dva hody, bylo by P(F | O₁, O₂) = α₂(F)/α₂(P) + α₂(F) = 0,637. Trojka, která padla potom, je na falešné kostce méně pravděpodobná, a tak podezření na druhý hod zpětně sníží na 0,534.

Nechť λ je skrytý Markovův model se stavy S, O = (O₁, …, O_T) pozorování a γ_t(i) aposteriorní pravděpodobnosti stavů z předchozí definice. Aposteriorní dekódování je posloupnost x̂ = (x̂₁, …, x̂_T) ∈ S^T, kde

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

(stav se volí v každém čase zvlášť). Nechť dále x ∈ S^T je skutečná posloupnost stavů a y ∈ S^T libovolné dekódování (aposteriorní nebo Viterbiho). Přesnost dekódování y je podíl časů, ve kterých se shoduje se skutečností:

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

Aposteriorní dekódování maximalizuje očekávaný počet správně určených hodů, Viterbi pravděpodobnost, že je správně celá cesta naráz. Proto může aposteriorní posloupnost obsahovat i přechod, který má v modelu pravděpodobnost nula — u kasina takový není, ale obecně to nastat může.

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

Vstup: počet skrytých stavů N ≥ 1, množina symbolů V, pozorování O = (O₁, …, O_T) ∈ V^T s T ≥ 2, počáteční odhad parametrů λ⁽⁰⁾ = (π, A, B) s P(O | λ⁽⁰⁾) > 0 a počet iterací K ≥ 1 (nebo tolerance ε > 0).
Výstup: parametry λ^(K) s věrohodností P(O | λ^(K)) ≥ P(O | λ⁽⁰⁾).

Pro k = 0, 1, …, K − 1 opakuj s parametry λ^(k):

1. Krok E: dopředně-zpětným algoritmem spočítej α_t(i), β_t(i), P(O | λ^(k)), aposteriorní pravděpodobnosti γ_t(i) pro t = 1, …, T a pravděpodobnosti přechodů ξ_t(i, j) = P(X_t = i, X_t+1 = j | O, λ^(k)),   t = 1, …, T − 1,   i, j ∈ S.
2. Krok M: nové parametry λ^(k+1) jsou pro všechna 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. Zastavení: skonči po K iteracích, případně dřív, když log P(O | λ^(k+1)) − log P(O | λ^(k)) < ε.

Platí P(O | λ^(k+1)) ≥ P(O | λ^(k)) — věrohodnost v žádné iteraci neklesne. Algoritmus konverguje ke stacionárnímu bodu věrohodnosti, v praxi obvykle k lokálnímu maximu, které nemusí být globální; výsledek proto závisí na počátečním odhadu.