<jupyter_start><jupyter_text>Markov ModelsSome resources I found:* [A Tutorial on HMM and Selected Applications in Speech Recognition - Rabiner](http://www.ece.ucsb.edu/Faculty/Rabiner/ece259/Reprints/tutorial%20on%20hmm%20and%20applications.pdf)* [Hidden Markov Models - Brown](http://cs.brown.edu/research/ai/dynamics/tutorial/Documents/HiddenMarkovModels.html)* https://www.autonlab.org/_media/tutorials/hmm14.pdf* http://www.kanungo.com/software/hmmtut.pdf* Wiki articles about Markov Models are actually very detailedAlgorithms connected to HMM:* Forward algorithm - calculate the probability that a certain sequence of observations occurs* Viterbi algorithm - find the most probable sequence of hidden states given a sequence of observations* Baum-Welch algorithm - find the parameters of a model that maximize the likelihood of a certain sequence of observations occuringThe model must have a Markov property, which means that the distribution of the next hidden state and observation can only depend on the current state or a fixed number of past states. The model is also uniform, which means that the conditional distributions don't change. It should also have a property that the probability of getting from one state to any other one is always non zero.The model has weights determining the distribution of hidden states in each hidden state, weights for distributions of observations in each given state and also an initial state (or a vector of probabilities of initial state).The hidden states seem to be usually discrete, while the observations can be discrete continuous.Number of hidden states grow exponentially with a number of features that are modelled. That is, if the model of a person speech has N states and you need to additionaly include information whether the person is male or female, you now need 2N states.<jupyter_code>import functools as ft
import itertools as it
import json
import math
import operator as op
import os
import random

from IPython.display import display
from ipywidgets import interact, interact_manual, widgets
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import misc, stats
from sklearn import metrics<jupyter_output><empty_output><jupyter_text>Discrete ObservationsFirst I will consider a simpler case where hidden states and observations are discrete. Markov model is defined by choosing a number of hidden features $h$ and visible observations $v$ (we don't care about assigning symbols to them and will use numbers $0 \dots h$ and $0 \dots v$ as states. It's then necessary to define two matrices of probabilities.First is a matrix of probabilities of transitions between hidden states $H$ of size $h \times h$ where $h$ is a number of hidden states. $H_{i, j}$ is probability of going from hidden state $i$ to hidden state $j$, for all $i$ $\sum_{j=0}^{h-1} H_{i, j} = 1$ and for all $i, j$ $H_{i, j} > 0$.Second matrix is a matrix of probabilities of emitting observations while being in a certain state $V$ of size $h \times v$ where $v$ is a number of possible observations (visible states). $V_{i, k}$ is a probability of emitting observation $k$ while in hidden state $i$. For all $i$ $\sum_{k=0}^{v-1} V_{i, k} = 1$<jupyter_code>hidden_n = 4
visible_n = 3

# square matrix h x h, each row sums to 1
hidden_weights = np.array([[0.1, 0.4, 0.4, 0.1], 
                           [0.2, 0.2, 0.2, 0.4], 
                           [0.4, 0.2, 0.3, 0.1], 
                           [0.3, 0.3, 0.1, 0.3]]) 
# matrix h x v, each row sums to 1
visible_weights = np.array([[0.3, 0.4, 0.3], 
                            [0.1, 0.9, 0.0], 
                            [0.5, 0.4, 0.1],
                            [0.25, 0.1, 0.65]]) 

def hidden_selected(hidden_n, index):
    hidden_states = np.zeros(hidden_n, dtype=np.float64)
    hidden_states[index] = 1.0
    return hidden_states

print(hidden_selected(hidden_n, 2))

def hidden_uniform(hidden_n):
    return np.ones(hidden_n, dtype=np.float64) / hidden_n

print(hidden_uniform(hidden_n))

initial_hidden = np.array([0.1, 0.1, 0.7, 0.1])
print(initial_hidden)
def generate_visible(visible_weights, hidden_state):
    return hidden_weights @ hidden_state

print(generate_visible(visible_weights, hidden_selected(hidden_n, 2)))
print(generate_visible(visible_weights, hidden_uniform(hidden_n)))
print(generate_visible(visible_weights, initial_hidden))
def walk_max_likely(hidden_weights, visible_weights, initial_hidden):
    max_hidden_per_hidden = np.argmax(hidden_weights, axis=1)
    max_visible_per_hidden = np.argmax(visible_weights, axis=1)
    
    hidden_max = np.argmax(initial_hidden)
    while True:
        visible_max = max_visible_per_hidden[hidden_max]
        yield visible_max, hidden_max
        hidden_max = max_hidden_per_hidden[hidden_max]
        
walk = walk_max_likely(hidden_weights, visible_weights, initial_hidden)
print(list(it.islice(walk, 25)))


def walk_random(hidden_weights, visible_weights, initial_hidden):
    hidden_range = np.arange(hidden_weights.shape[1])
    visible_range = np.arange(visible_weights.shape[1])
    
    hidden_chosen = np.random.choice(hidden_range, p=initial_hidden)
    while True:
        visible_chosen = np.random.choice(visible_range, p=visible_weights[hidden_chosen, :])
        yield visible_chosen, hidden_chosen
        hidden_chosen = np.random.choice(hidden_range, p=hidden_weights[hidden_chosen, :])
        
walk = walk_random(hidden_weights, visible_weights, initial_hidden)
print(list(it.islice(walk, 25)))<jupyter_output>[(0, 2), (1, 0), (1, 1), (2, 3), (1, 0), (1, 1), (2, 3), (1, 0), (1, 1), (2, 3), (1, 0), (1, 1), (2, 3), (1, 0), (1, 1), (2, 3), (1, 0), (1, 1), (2, 3), (1, 0), (1, 1), (2, 3), (1, 0), (1, 1), (2, 3)]
[(2, 3), (2, 3), (1, 1), (0, 2), (1, 2), (0, 2), (0, 0), (0, 2), (1, 1), (0, 3), (2, 0), (0, 2), (0, 2), (1, 1), (2, 2), (2, 2), (1, 1), (0, 0), (1, 1), (0, 2), (2, 2), (1, 1), (2, 3), (1, 1), (2, 3)]<jupyter_text>Forward AlgorithmWiki articles on this are actually quite long and good: [wiki/Baum-Welch](https://en.wikipedia.org/wiki/Baum%E2%80%93Welch_algorithm), [wiki/Forward algorithm](https://en.wikipedia.org/wiki/Forward_algorithm).Calculate probability of a given state given that a certain sequence of observations occured in the past.To be exact it calculates $P(x_t, y_t \dots y_1)$. It does so with recursion $$P(x_t, y_t \dots y_1) = \sum_{x_{t-1}} P(x_t, x_{t-1}, y_t \dots y_1) = \sum_{x_{t-1}} P(y_t | x_t, x_{t-1}, y_{t-1} \dots y_1) P(x_t | x_{t-1}, y_{t-1} \dots y_1) P(x_{t-1}, y_{t-1} \dots y_1)$$Now, thanks to Markov property we can simplify$$P(y_t | x_t, x_{t-1}, y_{t-1} \dots y_1) = P(y_t | x_t) = V_{x_t, y_t}$$  $$P(x_t | x_{t-1}, y_{t-1} \dots y_1) = P(x_t | x_{t-1}) = H_{x_{t-1}, x_t}$$And the final equation is$$P(x_t, y_t \dots y_1) = V_{x_t, y_t} \sum_{x_{t-1}} H_{x_{t-1}, x_t} P(x_{t-1}, y_{t-1} \dots y_1) $$Where $P(x_{t-1}, y_{t-1} \dots y_1)$ can be recursively calculated in the same way until we get to $P(x_0, \emptyset) = \pi_1$ and $\pi$ is a given vector of size $h$ defining the initial state.Using Markov property reduces complexity from $\theta(nh^n)$ to $\theta(nh^2)$, where $n$ is a number of observations.I keep all intermediate states in my implementation to visualize the progress better. I don't really plan to make a production grade system here.<jupyter_code>def forward_probability(hidden_weights, visible_weights, initial_hidden, observations):
    joint_probs = np.zeros((observations.size + 1, initial_hidden.size), dtype=np.float64)
    joint_probs[0, :] = initial_hidden
    
    for i, observation in enumerate(observations):
        joint_probs[i + 1] = (joint_probs[i] @ hidden_weights) * visible_weights[:, observation]
            
    return joint_probs[-1].sum(), joint_probs
    
walk = walk_random(hidden_weights, visible_weights, initial_hidden)
observations = np.array([v for v, h in it.islice(walk, 16)])
print(*forward_probability(hidden_weights, visible_weights, initial_hidden, observations), sep='\n')

observations = np.array([2, 1, 0])
print(*forward_probability(hidden_weights, visible_weights, hidden_selected(4, 1), observations), sep='\n')

observations = np.array([2, 1, 0])
print(*forward_probability(hidden_weights, visible_weights, hidden_selected(4, 0), observations), sep='\n')<jupyter_output>3.56039699457e-07
[[  1.00000000e-01   1.00000000e-01   7.00000000e-01   1.00000000e-01]
 [  1.36000000e-01   2.07000000e-01   1.12000000e-01   1.50000000e-02]
 [  4.17200000e-02   1.10430000e-01   5.23600000e-02   1.12100000e-02]
 [  2.02260000e-02   4.73481000e-02   2.22412000e-02   5.69430000e-03]
 [  8.83879600e-03   2.13448950e-02   9.92072400e-03   2.48942500e-03]
 [  3.94719028e-03   9.48192273e-03   4.41186284e-03   1.11607375e-03]
 [  1.75626833e-03   4.22320982e-03   1.96417075e-03   4.96349653e-04]
 [  7.81936798e-04   1.87999951e-03   8.74414195e-04   2.21023273e-04]
 [  2.61079872e-04   0.00000000e+00   9.73201207e-05   6.39562225e-04]
 [  1.02761881e-04   2.84188177e-04   7.90336831e-05   2.27708667e-05]
 [  4.22234227e-05   1.08522346e-04   4.95718318e-05   1.38686087e-05]
 [  1.99664507e-05   4.74019085e-05   2.19408995e-05   5.67490465e-06]
 [  6.58675739e-06   0.00000000e+00   2.46167223e-06   1.61550804e-05]
 [  1.94696062e-06   0.00000000e+00   4.98871266e-07   3.73[...]<jupyter_text>Backward AlgorithmIt's similar to forward algorithm. It calculates $P(y_{t+1} \dots y_T | x_t)$ that is probability of future sequence of observations $y_{t+1} \dots y_T$ given initial state $x_t$ at time $t$. It starts with the final state $x_T$.$$P(y_{t+1} \dots y_{T} | x_t)= \sum_{x_{t+1}} P(y_{t+1} \dots y_{T}, x_{t+1} | x_t)= \sum_{x_{t+1}} P(y_{t+1} | y_{t+2} \dots y_{T}, x_{t+1}, x_t) P(y_{t+2} \dots y_{T} | x_{t+1}, x_t) P(x_{t+1} | x_t)$$Now we can hopefully simplify:$$P(y_{t+1} | y_{t+2} \dots y_{T}, x_{t+1}, x_t) = P(y_{t+1} | x_{t+1}) = V_{x_{t+1},y_{t+1}}$$$$P(y_{t+2} \dots y_{T} | x_{t+1}, x_t) = P(y_{t+2} \dots y_{T} | x_{t+1})$$$$P(x_{t+1} | x_t) = H_{x_t,x_{t+1}}$$And the end result is:$$P(y_{t+1} \dots y_{T} | x_t) =  \sum_{x_{t+1}} V_{x_{t+1},y_{t+1}} H_{x_t,x_{t+1}} P(y_{t+2} \dots y_{T} | x_{t+1})$$The probabilities are computed recursively. Base case is $P(\emptyset | x_T) = 1$<jupyter_code>def backward_probability(hidden_weights, visible_weights, observations):
    conditional_probs = np.zeros((observations.size + 1, hidden_weights.shape[0]), dtype=np.float64)
    conditional_probs[-1, :] = 1
    
    for i, observation in reversed(list(enumerate(observations))):
        conditional_probs[i] = (conditional_probs[i + 1] * visible_weights[:, observation]) @ hidden_weights.T
            
    return conditional_probs[0].sum(), conditional_probs
    
walk = walk_random(hidden_weights, visible_weights, initial_hidden)
observations = np.array([v for v, h in it.islice(walk, 16)])
print(*backward_probability(hidden_weights, visible_weights, observations), sep='\n')

observations = np.array([2, 2, 2, 2])
print(*backward_probability(hidden_weights, visible_weights, observations), sep='\n')

observations = np.array([1, 1, 1, 1])
print(*backward_probability(np.eye(3), np.array([[0.2, 0.8], [1.0, 0.0], [0.5, 0.5]]), observations), sep='\n')<jupyter_output>5.32360852066e-07
[[  1.44303052e-07   1.26907056e-07   1.51407262e-07   1.09743482e-07]
 [  4.63201642e-07   4.40443395e-07   5.15501399e-07   3.87559496e-07]
 [  1.75650861e-06   1.24308285e-06   1.61848723e-06   1.48342499e-06]
 [  4.16099094e-06   2.75610462e-06   3.53103658e-06   3.29054555e-06]
 [  7.77487839e-06   6.83753007e-06   8.30460889e-06   5.97475557e-06]
 [  2.64666440e-05   2.32697110e-05   2.77374648e-05   2.01039074e-05]
 [  8.46054189e-05   8.06485493e-05   9.46139183e-05   7.11902301e-05]
 [  3.25613410e-04   2.27982736e-04   2.94741055e-04   2.70779844e-04]
 [  7.31556909e-04   5.12128753e-04   6.61900222e-04   6.08074760e-04]
 [  1.64103516e-03   1.15025128e-03   1.48839188e-03   1.36823397e-03]
 [  3.71326108e-03   2.58607620e-03   3.31852224e-03   3.05537255e-03]
 [  8.00821555e-03   5.83637700e-03   7.62942115e-03   7.11293505e-03]
 [  2.11782050e-02   1.30002200e-02   1.45034850e-02   1.60450550e-02]
 [  1.82995000e-02   4.28100000e-02   2.91215000e-02   3.75[...]<jupyter_text>Forward-Backward AlgorithmFor any $k = 1 \dots t$* Forward algorithm calculates $P(x_t, y_t \dots y_1)$. It calculates probability of a hidden state and a sequence of past observations at time $t$.* Backward algorithm calculates $P(y_T \dots y_{t+1} | x_t)$. It calculates probability of a sequence of future observations $y_{t+1} \dots y_T$ given a current state $x_t$ at time $t$.* Forward-backward algorithm calculates $P(x_t | y_T \dots y_1)$. So it gives probability of any hidden state $x_t$ at any time taking into accout full sequence of observations. It can then at any time select the most likely state. Note that by selecting the most likely state we will not necessarily get the most likely sequence. Viterbi algorithms is needed for that.Forward-backward algorithm combines the results of the previous two algorithms.$$P(x_t, y_t \dots y_1) P(y_T \dots y_{t+1} | x_t) = P(y_t \dots y_1 | x_t) P(y_T \dots y_{t+1} | x_t) P(x_t) = P(y_T \dots y_1 | x_t) P(x_t) = P(y_T \dots y_1, x_t) = P(x_t | y_T \dots y_1) P(y_T \dots y_1)$$$$P(x_t | y_T \dots y_1) = \frac{P(x_t, y_t \dots y_1) P(y_T \dots y_{t+1} | x_t)}{P(y_T \dots y_1)}$$$P(y_T \dots y_1)$ can be calculated as a marginal distribution of $\sum_{x_T} P(x_T, y_T \dots y_1)$<jupyter_code>def forward_backward_probability(hidden_weights, visible_weights, initial_hidden, observations):
    _, forward_result = forward_probability(hidden_weights, visible_weights, initial_hidden, observations)
    _, backward_result = backward_probability(hidden_weights, visible_weights, observations)
    scaling_coefficient = forward_result[-1, :].sum()
    combined_result = forward_result * backward_result
    return combined_result / scaling_coefficient, combined_result

walk = walk_random(hidden_weights, visible_weights, initial_hidden)
observations = np.array([v for v, h in it.islice(walk, 16)])
print(*forward_backward_probability(hidden_weights, visible_weights, initial_hidden, observations), sep='\n')
print(forward_backward_probability(hidden_weights, visible_weights, initial_hidden, observations)[0].sum(axis=1))
print()

observations = np.array([2, 1, 0])
print(*forward_backward_probability(hidden_weights, visible_weights, hidden_selected(4, 1), observations), sep='\n')
print()

observations = np.array([2, 1, 0])
print(*forward_backward_probability(hidden_weights, visible_weights, hidden_selected(4, 0), observations), sep='\n')
print()<jupyter_output>[[ 0.10186777  0.0896863   0.73204128  0.07640464]
 [ 0.3085618   0.07162325  0.51207468  0.10774027]
 [ 0.35582238  0.07491621  0.44794118  0.12132024]
 [ 0.12190915  0.62833388  0.21392863  0.03582834]
 [ 0.31029617  0.          0.10113975  0.58856408]
 [ 0.21377889  0.56187313  0.18477054  0.03957744]
 [ 0.25325211  0.06750127  0.44891593  0.23033068]
 [ 0.29361925  0.42869215  0.24106035  0.03662825]
 [ 0.2410574   0.43381504  0.26596914  0.05915842]
 [ 0.27110046  0.42104765  0.2481342   0.0597177 ]
 [ 0.2062698   0.48005598  0.2635385   0.05013572]
 [ 0.28070503  0.06618658  0.44342586  0.20968253]
 [ 0.25175911  0.47043376  0.24918719  0.02861994]
 [ 0.23282503  0.08033678  0.5022377   0.18460049]
 [ 0.32057354  0.09561011  0.45672954  0.12708682]
 [ 0.17694599  0.14702262  0.49377463  0.18225676]
 [ 0.3964654   0.          0.13615381  0.46738079]]
[[  5.26025456e-09   4.63122702e-09   3.78011937e-08   3.94538757e-09]
 [  1.59335335e-08   3.69848583e-09   2.64425446e-08   5.5634[...]<jupyter_text>Baum-Welch AlgorithmEstimate the most likely weights $H$, $V$ and $\pi$ (initial state) of a model given sequences of observations.It fits the weights using Expectation Maximization. The forward-backward algorithm is used to find the likelihoods of hidden states given a sequence of observations.We proved previously that forward-backward algorithm can calculate $P(x_t | y_T \dots y_1)$ as:$$P(x_t | y_T \dots y_1) = \frac{P(x_t, y_t \dots y_1) P(y_T \dots y_{t+1} | x_t)}{P(y_T \dots y_1)}$$where $P(x_t, y_t \dots y_1)$ $P(y_T \dots y_{t+1} | x_t)$ are given to us by forward and backward algorithms and $P(y_T \dots y_1)$ can be calculated as a marginal distribution of $\sum_{x_T} P(x_T, y_T \dots y_1)$To perform Baum-Welch learning we additionally need $P(x_t, x_{t+1} | y_T \dots y_1)$, which is the probability of transitioning from state $x_t$ at time $t$ to state $x_{t+1}$ at time $t+1$ given the sequence of observations. (and implicitly given model parameters)$$P(x_t, x_{t+1} | y_T \dots y_1)= \frac{P(x_t, x_{t+1}, y_T \dots y_1)}{P(y_T \dots y_1)}$$$$P(x_t, x_{t+1}, y_T \dots y_1) \\= P(y_T \dots y_1 | x_t, x_{t+1}) P(x_{t+1} | x_t) P(x_t) \\= P(y_T \dots y_{t+2} | x_t, x_{t+1}) P(y_{t+1} \dots y_1 | x_t, x_{t+1}) P(x_{t+1} | x_t) P(x_t)$$$$P(y_{t+1} \dots y_1 | x_t, x_{t+1}) P(x_t) = P(y_{t+1} \dots y_1, x_t | x_{t+1})= P(y_{t+1} | x_{t+1}) P(y_t \dots y_1, x_t | x_{t+1})$$$$P(x_t, x_{t+1}, y_T \dots y_1) \\= P(y_T \dots y_{t+2} | x_t, x_{t+1}) P(y_{t+1} | x_{t+1}) P(y_t \dots y_1, x_t | x_{t+1}) P(x_{t+1} | x_t) \\= P(y_T \dots y_{t+2} | x_{t+1}) P(y_{t+1} | x_{t+1}) P(y_t \dots y_1, x_t) P(x_{t+1} | x_t)$$$$P(x_t, x_{t+1} | y_T \dots y_1)= \frac{P(y_T \dots y_{t+2} | x_{t+1}) P(y_t \dots y_1, x_t) V_{x_{t+1}, y_{t+1}} H_{x_t, x_{t+1}}}{P(y_T \dots y_1)}$$Then the updates are performed like this:$$\pi_i = P(X_0 = x_i | y_T \dots y_1)$$$$H_{x_t, x_{t+1}} = \frac{\sum_{t=0}^{T-1} P(x_t, x_{t+1} | y_T \dots y_1)}{\sum_{t=0}^{T-1} P(x_t | y_T \dots y_1)}$$$$V_{x_t, y_t} = \frac{\sum_{t=1}^T \mathbb{1}_{Y_t = y_t} P(x_t | y_T \dots y_1)}{\sum_{t=1}^T P(x_t | y_T \dots y_1)}$$<jupyter_code>def forward_backward_baum_welch(hidden_weights, visible_weights, initial_hidden, observations):
    _, forward_result = forward_probability(hidden_weights, visible_weights, initial_hidden, observations)
    _, backward_result = backward_probability(hidden_weights, visible_weights, observations)
    scaling_coefficient = forward_result[-1, :].sum()
    
    per_state_result = forward_result * backward_result # observation+1 x state
    per_state_result /= scaling_coefficient
    
    # state x state (we can sum over time/observations)
#     per_pair_result = (forward_result[:-1].T @ (backward_result[1:] * visible_weights[:, observations].T)) * hidden_weights
    per_pair_result = ((backward_result[1:] * visible_weights[:, observations].T).T @ forward_result[:-1]) * hidden_weights.T
    per_pair_result /= scaling_coefficient
    
    return per_state_result, per_pair_result

def estimate_params_baum_welch(per_state_result, per_pair_result, observations, visible_n):
    hidden_n = per_pair_result.shape[0]
    
    new_initial_hidden = per_state_result[0]
    
    per_state_rest = per_state_result[1:]
    per_state_rest_total = per_state_rest.sum(axis=0)
    
    new_hidden_weights = per_pair_result / per_state_rest_total[:, np.newaxis]
    
    new_visible_weights = np.zeros((hidden_n, visible_n), dtype=np.float64)
    for i in range(visible_n):
        new_visible_weights[:, i] = per_state_rest[observations == i].sum(axis=0) / per_state_rest_total
        
    return new_hidden_weights, new_visible_weights, new_initial_hidden

walk = walk_random(hidden_weights, visible_weights, initial_hidden)
observations = np.array([v for v, h in it.islice(walk, 16)])
result = forward_backward_baum_welch(hidden_weights, visible_weights, initial_hidden, observations)
print(*result, result[0].sum(axis=1), result[1].sum(axis=1), sep='\n')
print()

new_hidden_weights, new_visible_weights, new_initial_hidden = \
    estimate_params_baum_welch(*result, observations, visible_weights.shape[1])

print('new initial', new_initial_hidden, new_initial_hidden.sum(), sep='\n')
print('new hidden', new_hidden_weights, new_hidden_weights.sum(axis=0), new_hidden_weights.sum(axis=1), sep='\n')
print('new visible', new_visible_weights, new_visible_weights.sum(axis=0), new_visible_weights.sum(axis=1), sep='\n')
def grouper(iterable, n):
    args = [iter(iterable)] * n
    return zip(*args)

def estimate_weights_baum_welch(hidden_n, visible_n, observations_batch, iterations, batch_n=None, 
                                hidden_weights=None, visible_weights=None, initial_hidden=None):
    if hidden_weights is None:
        hidden_weights = np.random.uniform(size=(hidden_n, hidden_n))
        hidden_weights /= hidden_weights.sum(axis=1)[:, np.newaxis]
    if visible_weights is None:
        visible_weights = np.random.uniform(size=(hidden_n, visible_n))
        visible_weights /= visible_weights.sum(axis=1)[:, np.newaxis]
    if initial_hidden is None:
        initial_hidden = np.random.uniform(size=hidden_n)
        initial_hidden /= initial_hidden.sum()
    batch_n = batch_n or len(observations_batch)
        
    batch_hidden_weights = np.zeros((batch_n, hidden_n, hidden_n), dtype=np.float64)
    batch_visible_weights = np.zeros((batch_n, hidden_n, visible_n), dtype=np.float64)
    batch_initial_hidden = np.zeros((batch_n, hidden_n), dtype=np.float64)
    shuffled_observations = observations_batch[:]
    for i in range(iterations):
        random.shuffle(shuffled_observations)
        observation_groups = grouper(shuffled_observations, batch_n)
        for observation_group in observation_groups:
            for j, observations in enumerate(observation_group):
                per_state_result, per_pair_result = \
                    forward_backward_baum_welch(hidden_weights, visible_weights, initial_hidden, observations)
                batch_hidden_weights[j], batch_visible_weights[j], batch_initial_hidden[j] = \
                    estimate_params_baum_welch(per_state_result, per_pair_result, observations, visible_n)
            hidden_weights = batch_hidden_weights.mean(axis=0) 
            visible_weights = batch_visible_weights.mean(axis=0) 
            initial_hidden = batch_initial_hidden.mean(axis=0)
            
    return hidden_weights, visible_weights, initial_hidden
def test_baum_welch(hidden_weights, visible_weights, initial_hidden, observations_n, 
                    observations_min, observations_max, iterations, bw_kwargs):
    observations_batch = []
    for i in range(observations_n):
        walk = walk_random(hidden_weights, visible_weights, initial_hidden)
        length = int(np.random.uniform(observations_min, observations_max))
        observations = np.array([v for v, h in it.islice(walk, length)])
        observations_batch.append(observations)

    print('Generator', hidden_weights, visible_weights, initial_hidden, sep='\n', end='\n\n')

    hidden_n, visible_n = visible_weights.shape
    result = estimate_weights_baum_welch(hidden_n, visible_n, observations_batch, 
                                         iterations=iterations, **bw_kwargs)

    print('Result', *result, sep='\n', end='\n\n')
test_baum_welch(hidden_weights, visible_weights, initial_hidden, 100, 30, 40, 100)
test_baum_welch(hidden_weights, visible_weights, initial_hidden, 100, 30, 40, 100, bw_kwargs={'batch_n': 1})
test_baum_welch(hidden_weights, visible_weights, initial_hidden, 10, 300, 400, 10, 
                bw_kwargs={'batch_n': 1,
                           'hidden_weights': hidden_weights, 
                           'visible_weights': visible_weights, 
                           'initial_hidden': initial_hidden})
test_baum_welch(np.array([[0.5, 0.5], [0.5, 0.5]]), 
                np.array([[0.5, 0.5], [0.5, 0.5]]), 
                np.array([0.5, 0.5]),
                50, 50, 100, 100)
test_baum_welch(np.array([[0.9, 0.1], [0.1, 0.9]]), 
                np.eye(2), 
                np.array([0.5, 0.5]),
                6, 400, 500, 100, bw_kwargs={'batch_n': 2})<jupyter_output>Generator
[[ 0.9  0.1]
 [ 0.1  0.9]]
[[ 1.  0.]
 [ 0.  1.]]
[ 0.5  0.5]

Result
[[ 0.91131163  0.08868837]
 [ 0.09903204  0.90096796]]
[[ 0.99875008  0.00124992]
 [ 0.00479634  0.99520366]]
[ 0.49701092  0.50298908]<jupyter_text>It seems not to work very well. Well then.  Viterbi Algorithmhttps://www.vocal.com/echo-cancellation/viterbi-algorithm-in-speech-enhancement-and-hmm/Calculate the most likely sequence of hidden states given a sequence of observations.Using forward-backward algorithm we calculated probabilities for all states at all times given a sequence of observations P(x_t | y_T \dots y_1). We could just take the most probable state at each time but it would not necessarily give us the most probable sequence. We could select two states that were most likely separately, but that are very unlikely to occur together because the probability of transitioning from the first one to the second is very low or impossible (the corresponding probability in matrix $H$ is very low). Viterbi algorithm solves problem of finding a sequence $x_T \dots x_1$ that maximizes $P(x_T \dots x_1 | y_T \dots y_1)$. It uses a dynamic programming technique. It builds two matrices, let's call them $A$ and $B$. $A_{t,x_i} = max_{x_1 \dots x_{t-1}} P(x_1 \dots x_{t-1}, X_t = x_i, y_1 \dots y_t)$ so it contains the highest probability of single sequence of states $x_1 \dots x_{t-1}$ occuring and ending in state $x_i$ joint with a sequence of observations $y_1 \dots y_t$ occuring.$B_{t,x_i} = argmax_{x_1 \dots x_{t-1}} P(x_1 \dots x_{t-1}, X_t = x_i, y_1 \dots y_t)$ and it is used to recordwhich states were chosen as the states of maximum probability while calculating $A$. It is then used to recover the most likely path at the end. Note that $t > 1$, because at $t = 1$ we do not choose anything while calculating values in $A$.Given A_{t,x_i} we can recursively calculate probabilities for next time steps A_{t+1,x_i} as follows:$$P(x_1 \dots x_t, y_1 \dots y_t) \\= P(y_t | x_1 \dots x_t, y_1 \dots y_{t-1}) P(x_1 \dots x_t, y_1 \dots y_{t-1}) \\= P(y_t | x_1 \dots x_t, y_1 \dots y_{t-1}) P(x_t | x_1 \dots x_{t-1}, y_1 \dots y_{t-1}) P(x_1 \dots x_{t-1}, y_1 \dots y_{t-1}) \\= P(y_t | x_t) P(x_t | x_{t-1}) P(x_1 \dots x_{t-1}, y_1 \dots y_{t-1}) \\= V_{x_t, y_t} H_{x_{t-1}, x_t} P(x_1 \dots x_{t-1}, y_1 \dots y_{t-1})$$$$A_{t,x_i} \\= max_{x_1 \dots x_{t-1}} P(x_1 \dots x_{t-1}, X_t = x_i, y_1 \dots y_t) \\= max_{x_1 \dots x_{t-1}} V_{x_i, y_t} H_{x_{t-1}, x_i} P(x_1 \dots x_{t-1}, y_1 \dots y_{t-1}) \\= V_{x_i, y_t} max_{x_{t-1}}(H_{x_{t-1}, x_i} A_{t-1, x_{t-1}})$$Note that the base case of this recursion is $A_{1,x_i} = P(X_1 = x_i, y_1) = \pi_i V_{x_i, y_1}$Also regarding B $$B_{t, x_i} = argmax_{x_1 \dots x_{t-1}} P(x_1 \dots x_{t-1}, X_t = x_i, y_1 \dots y_t) = argmax_{x_{t-1}} H_{x_{t-1}, x_i} A_{t-1, x_{t-1}}$$After computing whole $A$ we can find the most probable path by selecting the state x_i with highest probability at time T $max_{x_i} A_{T, x_i}$ and reconstruct the whole sequence with $B$.<jupyter_code>def most_likely_hidden_viterbi(hidden_weights, visible_weights, initial_hidden, observations):
    joint_probs = np.zeros((observations.size, initial_hidden.size), dtype=np.float64)
    max_probs_ix = np.zeros((observations.size - 1, initial_hidden.size), dtype=np.int32)
    
    for i, observation in enumerate(observations):
        if i == 0:
            joint_probs[0] = initial_hidden * visible_weights[:, observation]    
        else:
            next_state_probs = joint_probs[i - 1, np.newaxis] * hidden_weights
            next_state_argmax = next_state_probs.argmax(axis=1)
            next_state_max = next_state_probs.max(axis=1)
            
            joint_probs[i] = next_state_max * visible_weights[:, observation]
            max_probs_ix[i - 1] = next_state_argmax
    
    most_likely_path = np.zeros(observations.size, dtype=np.int32)
    most_likely_path[-1] = joint_probs[-1].argmax()
    for i in range(observations.size - 1, 0, -1):
        most_likely_path[i - 1] = max_probs_ix[i - 1, most_likely_path[i]]
    
    return most_likely_path, joint_probs, max_probs_ix
def test_viterbi(hidden_weights, visible_weights, initial_hidden, length):
    walk = walk_random(hidden_weights, visible_weights, initial_hidden)
    observations_states = np.array(list(it.islice(walk, length))).T
    print(observations_states)
    print(*most_likely_hidden_viterbi(hidden_weights, visible_weights, initial_hidden, observations_states[0]), sep='\n')

test_viterbi(np.eye(2), np.eye(2), np.array([0.5, 0.5]), 5)
print()

test_viterbi(np.array([[0.5, 0.5], [0.5, 0.5]]), np.array([[0.5, 0.5], [0.5, 0.5]]), np.array([0.5, 0.5]), 8)
print()

test_viterbi(hidden_weights, visible_weights, initial_hidden, 8)<jupyter_output>[[0 0 0 0 0]
 [0 0 0 0 0]]
[0 0 0 0 0]
[[ 0.5  0. ]
 [ 0.5  0. ]
 [ 0.5  0. ]
 [ 0.5  0. ]
 [ 0.5  0. ]]
[[0 0]
 [0 0]
 [0 0]
 [0 0]]

[[1 0 0 0 1 1 0 1]
 [1 0 0 1 0 1 1 1]]
[0 0 0 0 0 0 0 0]
[[  2.50000000e-01   2.50000000e-01]
 [  6.25000000e-02   6.25000000e-02]
 [  1.56250000e-02   1.56250000e-02]
 [  3.90625000e-03   3.90625000e-03]
 [  9.76562500e-04   9.76562500e-04]
 [  2.44140625e-04   2.44140625e-04]
 [  6.10351562e-05   6.10351562e-05]
 [  1.52587891e-05   1.52587891e-05]]
[[0 0]
 [0 0]
 [0 0]
 [0 0]
 [0 0]
 [0 0]
 [0 0]]

[[0 1 0 2 1 1 0 1]
 [2 1 0 3 1 1 3 0]]
[2 1 0 3 1 0 2 1]
[[  3.00000000e-02   1.00000000e-02   3.50000000e-01   2.50000000e-02]
 [  5.60000000e-02   6.30000000e-02   4.20000000e-02   3.50000000e-03]
 [  7.56000000e-03   1.26000000e-03   1.12000000e-02   4.72500000e-03]
 [  1.34400000e-03   0.00000000e+00   3.36000000e-04   1.47420000e-03]
 [  5.89680000e-05   5.30712000e-04   2.15040000e-04   4.42260000e-05]
 [  8.49139200e-05   9.55281600e-05   4.24569600[...]
{'fraction_non_alphanumeric': 0.1391987060686727, 'fraction_numerical': 0.1091417596118206, 'mean_word_length': 3.9323583180987205, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 25, 'https://': 4, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 19, 'word_list_counts': {'cursed_substrings.json': 30, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '314722', 'license': 'MIT', 'n_tokens_mistral': 13153, 'n_tokens_neox': 11384, 'n_words': 2832, 'path': 'notebooks/MarkovModels.ipynb', 'repo_name': 'michalsosn/paprotka', 'size': '29677'}