# CS 530 - Lecture 05

## Hidden Markov Models

Bernhard Firner

2026-02-05

---

## Review

* We've been discussing how to make decisions about discrete event sequences
* A sequence of observations, with probabilities assigned to each sequence, is called a Markov Chain
* Observe tokens, $X_1,...,X_T$
  * $T$ is the last step in the chain
  * $p(X_T) = p(X_1)\prod_{t=2}^{T}p(X_t|X_{t-1})$

---

## Problems

* We had some fun using Markov chains to predict future tokens in a sequence
* If we want to make better predictions, we need longer sequences
  * But the probability of any arbitrary sequence gets smaller with the size of the sequence
  * This means it is difficult to collect data for long sequences
  * This is related to the *curse of dimensionality*

---

## Solutions to Missing Data

* Add more data!
* Google has n-gram models fit to 5 by gathering one trillion words
  * 100GB text dataset
    * [https://research.google/blog/all-our-n-gram-are-belong-to-you/](https://research.google/blog/all-our-n-gram-are-belong-to-you/)
    * [Google Books Ngram Viewer](https://books.google.com/ngrams/)
* Or we could change our approach

---

## Better Probabilities

* Observations may be complicated
  * Could be continuous, or made up of related tokens
* So what if we considered the source of observed tokens rather than the tokens themselves?
* By considering the source of our observations, we can perform probabilistic reasoning
  * e.g. the *Viterbi algorithm* finds a most likely true state of the world given a sequence of observations

---

## Viterbi Algorithm

* A dynamic programming approach to find the most likely sequence of states
  * This is useful when our observations are only probabilistic
  * For example, we can use this with DNN predictions

---

## Other Techniques

* The Viterbi algorithm assumed that we had information about the *emission* and *transition* matrices
  * The emission matrix has the probabilities of observations, given a current state
  * The transition matrix has likelihoods of changing from one state to another
* Sometimes, we don't know those

---

## Hidden Markov Models

* We may not know much about a system
  * Perhaps a guess about the number of states
* Observations are `emitted` by the hidden states
* The hidden states may `transition` without us knowing

---

## Hidden States

* With K hidden states, observation likelihoods are a $K\times V$ matrix
  * K rows, one for each state
  * V columns, one for each token in our vocabulary
* Called the emission matrix, B

---

## Similarities

* If we are learning everything from scratch, this is a bit like clustering
  * Attempting to find similar patterns in the data
  * Will use an algorithm that maximizes the likelihood
* But similarities come from sequential observations, not sets of features

---

## Clustering and HMM

* If there were no sequential information
  * Hidden states would be clusters
  * Parameters would be optimized to match the observations through the iterative EM algorithm
* But there *are* temporal dependencies

---

## Quick Review

* Clustering
  * Is K-Means familiar?
  * How about mixture models?

---

## Probabilities with Hidden States

* Given a sequence of observations, we'll have to estimate three things:
  * $\pi$: the initial state distribution
  * A: the transition matrix
  * B: the emission matrix
* All together, we call those parameters $\Theta$

---

## The initial state distribution

* Probability of the hidden states
  * $\pi(i) = p(z_1=i)$
* Probabilities are just the relative frequency of the hidden states
  * So if we can estimate the likelihood of being in different states, we get this for free

---

## The Transition Matrix

* $A(i,j) = p(z_t=j|z_{t-1}=i)$
* Transitions are between hidden states
* One advantage over Markov chains is that observations won't be sparse
  * Remember that getting likelihoods for long word sequences was difficult

---

## The Emission Matrix

* $B(i,j) = p(x_i|z_j)$
  * The probability of seeing token $x_i$ given we are in hidden state $z_j$
* This generates the observations, but now it has no temporal dependence

---

## Observation Probability

* We want to maximize the likelihood of our observations
  * So what is the probability of an observation?
* Given an observation at time $t$, $x_t$:
  * What is the probability of seeing $x_t$ over all of the possible hidden states?
* Similar to the shared responsibility from clustering

---

## forwards-backwards algorithm

* The probability for each hidden state, $z$, given all observations is
  * $p(z_t=j|x_{1:T})$
* From $1:t$:
  * What are the hidden state probabilities given the past observations?
* From $t+1:T$
  * What are the probabilities we were in those hidden states given future observations?

---

## Sometimes unfair coin example

* Let's say we're flipping a coin
  * $p(tails)$ and $p(head)$ are both 0.5
* Unbeknownst to us, sometimes we grab an unfair coin from our pocket
  * $p(tails) = 0.1$ and $p(head) = 0.9$

<p><math display="block" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><mi>B</mi><mo>=</mo><mrow><mo stretchy="true" form="prefix">(</mo><mtable><mtr><mtd columnalign="center" style="text-align: center"><mn>0.5</mn></mtd><mtd columnalign="center" style="text-align: center"><mn>0.5</mn></mtd></mtr><mtr><mtd columnalign="center" style="text-align: center"><mn>0.1</mn></mtd><mtd columnalign="center" style="text-align: center"><mn>0.9</mn></mtd></mtr></mtable><mo stretchy="true" form="postfix">)</mo></mrow></mrow><annotation encoding="application/x-tex">B = 
  \begin{pmatrix}
    0.5 & 0.5\\
    0.1 & 0.9
  \end{pmatrix}</annotation></semantics></math></p>

---

## Unfair Coin Transitions

* The transition probabilities come from how many flips we do at a time and how many coins are in our pocket
* Let's say that they are

<p><math display="block" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><mi>A</mi><mo>=</mo><mrow><mo stretchy="true" form="prefix">(</mo><mtable><mtr><mtd columnalign="center" style="text-align: center"><mn>0.8</mn></mtd><mtd columnalign="center" style="text-align: center"><mn>0.2</mn></mtd></mtr><mtr><mtd columnalign="center" style="text-align: center"><mn>0.2</mn></mtd><mtd columnalign="center" style="text-align: center"><mn>0.8</mn></mtd></mtr></mtable><mo stretchy="true" form="postfix">)</mo></mrow></mrow><annotation encoding="application/x-tex">A = 
  \begin{pmatrix}
    0.8 & 0.2\\
    0.2 & 0.8
  \end{pmatrix}</annotation></semantics></math></p>

---

## Initial State Distribution

* We have the transition matrix, A, and the emission matrix B
* What is the initial state distribution?
  * Let's say we're equally likely to pull out the fair or unfair coins
  * $\pi = [0.5, 0.5]$

---

## Observation Probability

* We observe THH (tails and then two heads)
  * In vector form, $O = [[1, 0], [0, 1], [0, 1]]$
  * What is the probability of that heads at t?
* Begin at t=0, calculate for both hidden states
  * $\pi BO[0]$
    * $p(z_i=0|O(0)) = 0.5 \times 0.5 = 0.25$
    * $p(z_i=1|O(0)) = 0.5 \times 0.1 = 0.05$
  * $p(x_0) = [0.25, 0.05]$

---

## Forward Pass

* $BO[1] = [0.5, 0.9]$
* Need to multiply by the probability that we were in those hidden states
* Given likelihood of the hidden states at the previous time, the probability of being in state 0 or 1 now
  * $p(x_{t-1})ABO[t] = [0.25, 0.05]A[0.5, 0.9]$
  * $= [0.2+0.01,0.05+0.04][0.5, 0.9] = [0.105, 0.081]$

---

## Observation

* The first column of A captures the probability that we were in $z_0$ and stayed in $z_0$, or that we were in $z_1$ and switched
* This calculation will only ever depend upon the state at $t-1$
  * So the forward pass will take $T$ steps
    * $k$ operations at each step, where $k$ is the number of hidden states
  * $k^2T$ in total for the coin example

---

## Finished Forward Pass

* $\alpha = [[0.25, 0.05], [0.105, 0.081], [0.0501 , 0.07722]]$
* Tells us the hidden state probabilities given the past observations at each t
* Notice that long sequences will have tiny probabilities

---

## Backward Pass

* What are the probabilities we were in those hidden states given future observations?
* The initial probabilities begin at 1
* Iterate beginning at T-2, check transitions in pairs

---

## Backward Pass

* What's the probability of going from $O[-1]$ to $O[-2]$?
* $\beta = A B O[t] p(x_{t+1}|x_{t+1:T})$
  * $BO[t]$ is the probability of the observation for each hidden state
  * $A$ (the transition matrix) governs how likely each state is

---

## Backward Pass

* Initialize $p(x_{t+1}|x_{t+1:T})$ to 1 at T
* $\beta = A B O[t+1]1 = [0.58, 0.82]$
  * Makes sense; given $O[-1]$ was heads, hidden state 2 was more likely

---

## Forward and Backward

* $\alpha = [[0.25, 0.05], [0.105, 0.081], [0.0501 , 0.07722]]$
* $\beta = [[0.141288, 0.121032], [0.3364, 0.6724], [0.58, 0.82]]$
* These show the likelihood of $O$ for each hidden state
  * e.g. hidden states $[0, 1, 1]$ look the most probable
  * If the transition from fair to unfair was really low, that second state could stay fair

---

## Expectations

* We should combine $\alpha$ and $\beta$ to get expectations
* One expectation for the states, and one for the transitions
  * Remember we need to update A, B, and $\pi$
  * A is the state transitions
  * B is the emission probabilities
  * $\pi$ are the initial state probabilities

---

## Expectation Estimation

* There are other ways, but we'll use the following
* Going to normalize values by the chain probability
  * $p(O) = \prod(\alpha)$
  * We know that the observations occurred, so that is the maximum probability we can achieve
* Call state estimate $\gamma$, call transition estimate $\xi$

---

## Expectation continued

* $\gamma = \frac{\alpha\beta}{p(O)}$
  * Just the normalized probability of the observations given our current parameters
* $\xi$ will be size $[T,k,k]$
  * Transition probability at each time
  * Will be built from the forward pass probability, the current transition probability, the probability of the next observation, and the backward probabilities

---

## $\xi$ continued

* Recall $A B O[t]$ was how we estimated the probability of transitioning into a state and observing $O[t]$
* $\xi[t] \leftarrow \frac{\alpha[t] A B O[t] \beta[t+1]}{p(O)}$
  * Read this as:
    * the probability of the left chain
    * the probability of seeing the observation given our model
    * the probability of seeing this given the right chain
    * normalized by the entire chain probability

---

## Expectation-Maximization

* Getting $\gamma$ and $\xi$ finishes our expectation calculation
  * Analogous to estimating the normalized probability of all points during clustering
  * Here, the relative responsibilities of the hidden states for each point were calculated whenever we multiplied $A B O[t]$

---

## Maximization Updates

* $\pi \leftarrow \gamma[0]$
  * Use the initial state estimate as our initial state distribution
* $A \leftarrow \frac{\sum\xi}{\sum\gamma}$
  * For each $A_{jk}$
    * The expected number of transitions from j to k (from $\xi$, the $T\times N\times N$ matrix)
    * Divided by the expected number of transitions from j to anything

---

## Maximizing Emission (B)

* Compute the emission for each symbol in the vocabulary
* Sum the expected number of times we are in state $j$ and see symbol $l$
* Divide by expected number of times we are in state $j$
* Select indices with $\mathbb{I}(x_{i,t}=l)$, but we'll use a loop in reality
* Assume $j$ is the index of vocabulary $l$
* $B[j] \leftarrow \frac{\sum\gamma[j,(where O[t]=l)]}{\gamma[j]}$

---

## EM Algorithm

* Now repeat
* The EM algorithm for HMMs is called the Baum-Welch algorithm
* Concern:
  * Those chain probabilities will be small
  * So we'll change our math to use logarithms

---

## Forward Code

```python
def forward_pass_log(log_pi, log_A, log_B, observations):
    """Computes the forward probabilities (alpha) in log space.
    This is the first part of the equation below:
      p(z_t = j|x_{1:T}) = p(z_t = j|x_{1:t})p(x_{t+1:T}|z_t = j)
    """
    T = len(observations)
    # The number of hidden states
    N = log_A.shape[0]
    # Initialize with negative infinity (log(0))
    log_alpha = np.full((T, N), -np.inf)

# Initialization (t=0)
    # Begin with the probabilities of seeing the first observation for the emissions of each hidden state
    log_alpha[0, :] = log_pi + log_B[:, observations[0]]

# Iterate over the history (t=1 to T-1)
    # This fills in the forward probabilities at each value
    for t in range(1, T):
        # log_B_t is the log emission probability for the observation at time t
        # That's p(z_t = j | x_t)
        log_B_t = log_B[:, observations[t]]

# Loop over the hidden states
        # See 17.54 through 17.58 of Murphy's book
        for j in range(N):
            # Remember, non-log space these are multiplications
            # We are computing p(z_t = j|x_{1:t}), but with a history of 1
            # log_alpha[t-1, i] + log_A[i, j] + log_B_t[j]
            # Sum log probabilities (using np.logaddexp.reduce for numerical stability)
            log_sum = np.logaddexp.reduce(log_alpha[t-1, :] + log_A[:, j])
            log_alpha[t, j] = log_sum + log_B_t[j]

# Termination: Log-likelihood of the observation sequence
    log_likelihood = np.logaddexp.reduce(log_alpha[T-1, :])

return log_alpha, log_likelihood
```

---

## Backward Code

```python
def backward_pass_log(log_A, log_B, observations):
    """Computes the backward probabilities (beta) in log space.
    This is the second part of the equation below:
      p(z_t = j|x_{1:T}) = p(z_t = j|x_{1:t})p(x_{t+1:T}|z_t = j)
    """
    T = len(observations)
    # The number of hidden states
    N = log_A.shape[0]
    # Initialize with negative infinity (log(0))
    log_beta = np.full((T, N), -np.inf)

# Initialization (t=T-1): P(nothing | S_T-1) = 1, so log(1) = 0
    log_beta[T-1, :] = 0.0

# Iterate from (t=T-2 down to 0)
    # This fills in the backward probabilities at each value
    for t in range(T - 2, -1, -1):
        # log_B_t_plus_1 is the log emission probability for the observation at time t+1
        log_B_t_plus_1 = log_B[:, observations[t+1]]

for i in range(N):
            # Remember, non-log space these are multiplications
            # We are computing p(x_{t+1:T}|z_t = j), but only looking ahead 1
            # log_A[i, j] + log_B_t_plus_1[j] + log_beta[t+1, j]
            log_sum = log_A[i, :] + log_B_t_plus_1 + log_beta[t+1, :]
            log_beta[t, i] = np.logaddexp.reduce(log_sum)

return log_beta
```

---

## Expectation Code

```python
def compute_expectations_log(log_alpha, log_beta, log_A, log_B, observations, log_likelihood):
    """
    Computes gamma (single-state expectation) and xi (two-state transition expectation)
    in log space.
    These are from Murphy 17.101 and 17.102
      gamma is p(z_t = j|x_{i,1:T_i},\Theta)
      xi is p(z_{t-1} = j,z_{t}=k|x_{i,1:T_i},\Theta)
    Keep in mind that T is actually just as long as our history, although we'll
    loop over all values in this function.
    """
    # T is the number of observations, N is the number of hidden states
    T, N = log_alpha.shape

# Log P(O) is the log_likelihood from the forward pass
    log_PO = log_likelihood

# log_gamma[t, i] = log P(St=i | O)
    # Gamma is analogous to the soft clustering via assignment from mixture models
    # alpha*beta has the forward-backward probabilities, normalized by the likelihood of the observations
    log_gamma = log_alpha + log_beta - log_PO

# log_xi[t, i, j] = log P(St=i, St+1=j | O)
    log_xi = np.full((T - 1, N, N), -np.inf)

for t in range(T - 1):
        # Sum of logs for the numerator (P(St=i, St+1=j, O))
        # log(alpha_t(i)) + log(A_ij) + log(B_j(O_t+1)) + log(beta_t+1(j))
        log_B_t_plus_1 = log_B[:, observations[t+1]]

# For all N hidden states, this is the probability from the forward pass * transition probability * emission probability * backward probabilities
        # This is 17.66 from Murphy's book
        # These are the transition probabilities over every possible transition
        # Also called the (smoothed) two-slice marginal
        log_terms = (
            log_alpha[t, :, np.newaxis] +   # log alpha_t(n),    shape (N, 1)
            log_A +                         # log A_ij,          shape (N, N)
            log_B_t_plus_1[np.newaxis, :] + # log B_j(x_{t+1}),  shape (1, N)
            log_beta[t+1, np.newaxis, :]    # log beta_{t+1}(j), shape (1, N)
        )

# Divide (or subtract in log space) by the probability of this observation, normalizing this observation
        log_xi[t, :, :] = log_terms - log_PO

return log_gamma, log_xi
```

---

## Maximization Code

```python
def maximization_step(log_gamma, log_xi, observations, vocabulary):
    """
    Compute the new HMM parameters (pi_new, A_new, B_new) from the expectations.
    """
    T = len(observations)
    # Shape is T, N, N
    N = log_xi.shape[1]

# 1. Update Initial Probabilities (pi)
    # pi'_i = gamma_0(i)
    # From 17.103 in Murphy
    log_pi = log_gamma[0, :]

# 2. Update Transition Matrix (A)
    # Numerator: \sum_{t=0}^{T-2}(xi_t(i, j))
    # Denominator: \sum_{t=0}^{T-2}(gamma_t(i))
    # From 17.103 in Murphy
    log_xi_sum = np.logaddexp.reduce(log_xi)
    log_gamma_sum = np.logaddexp.reduce(log_gamma[:-1, :])

# For A, we take the number of expected transitions from j to k,
    # and divide by the expected transitions from j to any other state.
    # This guarantees that each row sums to 1.

# Use broadcasting for A_new[i, j] = xi_sum[i, j] / gamma_sum[i]
    # np.newaxis ensures the denominator broadcasts correctly across columns j
    # From 17.103 in Murphy
    log_A = log_xi_sum - log_gamma_sum[:, np.newaxis]

# 3. Update Emission Matrix (B)
    # Numerator: sum(gamma_t(i) * I(O_t = k))
    # Denominator: sum(gamma_t(i)) from t=0 to T-1 (total expected time in state i)
    # Numerator from 17.104 in Murphy
    # It is the expectation of the multinomial, E[M_{jl}]
    # This is conceptually simple. It's the expected number of times we observe symbol l while in state j
    # We normalize by dividing by the expected number of times we are in state j
    log_gamma_sum_all_t = np.logaddexp.reduce(log_gamma)

# If there was only one state, this would just be the expected number of
    # times we see symbol l divided by the total number of observations.
    # That would give us the emission probability into l for markov model
    # without a hidden state.

# Population the emission matrix. Initialize with 0 probabilities.
    log_B = np.zeros((N, len(vocabulary))) - np.inf

for k in vocabulary:
        # Find indices where observation O_t equals k
        indices_k = (observations == k)

if len(indices_k) > 0:
            # Sum log_gamma over time t where O_t=k
            log_gamma_sum_k = np.logaddexp.reduce(log_gamma[indices_k, :])
            # log(B'_jk) = log_gamma_sum_k[j] - log_gamma_sum_all_t[j]
            log_B[:, k] = log_gamma_sum_k - log_gamma_sum_all_t
        else:
            # If an observation never appears, its log-probability will be -inf (or very small)
            #log_B[:, k] = -np.inf
            # We already initialized these to 0
            pass

return log_pi, log_A, log_B
```

---

## Estimating Hidden States

* We want to understand a system from observations
* Simple example
  * Can we determine if a coin being flipped is fair or fake?

---

## Data Generation

```python
def getCoinObservations(k=50):
    # true emission probabilities for a fair coin
    fair_emission = [0.5, 0.5]
    # true emission probabilities for an unfair coin
    unfair_emission = [0.1, 0.9]

actual_emission_matrix = np.array([fair_emission, unfair_emission])

actual_transition_matrix = [ [0.8, 0.2],
                                 [0.2, 0.8]]

# Using the values above, we can generate some data

actual_states = []
    observations = []

cur_state = 0

# 0 for tail, 1 for heads
    vocabulary=list(range(2))

for i in range(k):
        # Generate the next state from the current state and the transition matrix
        new_state = random.choices([0, 1], weights=actual_transition_matrix[cur_state])[0]
        # Generate a die roll from the new state and the emission matrix
        new_roll = random.choices(vocabulary, weights=actual_emission_matrix[new_state])[0]
        # Remember the results
        actual_states.append(new_state)
        observations.append(new_roll)
        cur_state = new_state
    observations = np.array(observations)
    return actual_states, observations, vocabulary
```

-v-

## Full Code

```python
import argparse
import random
import numpy as np

## Variable names from Murphy's "Machine Learning: A Probabilistic Approach"

def getDieObservations(k=250):
    # true emission probabilities for a fair die
    fair_emission = [1.0/6] * 6
    # true emission probabilities for a loaded die that rolls a 6 50% of the time
    loaded_emission = [0.1]*5 + [0.5]

actual_emission_matrix = np.array([fair_emission, loaded_emission])

actual_transition_matrix = [ [0.95, 0.05],
                                 [0.05, 0.95]]

# Using the values above, we can generate some data

actual_states = []
    observations = []

cur_state = 0

# Die rolls should be 1-6, but we'll simplify some code and use 0-5
    vocabulary=list(range(6))

def getCoinObservations(k=50):
    # true emission probabilities for a fair coin
    fair_emission = [0.5, 0.5]
    # true emission probabilities for an unfair coin
    unfair_emission = [0.1, 0.9]

actual_emission_matrix = np.array([fair_emission, unfair_emission])

actual_transition_matrix = [ [0.8, 0.2],
                                 [0.2, 0.8]]

# Using the values above, we can generate some data

actual_states = []
    observations = []

cur_state = 0

# 0 for tail, 1 for heads
    vocabulary=list(range(2))

###############################################################
# Now see if we can learn what happened

# X are the observed variables, Z are the latent variables (the hidden states)

# The parameters of a hidden markov model are \Theta = (\pi, A, B)
# \pi(i) = p(z_1=i) (the initial state distribution)
# A(i,j) = p(z_t=j|z_{t-1}=i) (the transition matrix from the previous state to this one)
# B is the emission matrix for the observed variables, p(x_t|z_t|j)

# There are two hidden states
# We'll represent transitions between them with a transition matrix (often called A)
# We'll populate the values by assuming that all probabilities are equal
# The rows must sum to 1; these are all the ways we can leave that state.
# transition_matrix = np.full((2,2), 0.5)

# Now the emission matrix, B
# We would normally begin by counting the possible symbols,
# but we know that they are numbers from 1 to 6
# Begin by assuming some random values
# emission_matrix = np.random.rand(2,6)
# The rows must sum to 1
# emission_matrix = emission_matrix / emission_matrix.sum(axis=1, keepdims=True)

# Now we need pi, the initial state distribution
# Let's just assume that all hidden states are equally likely
# pi = np.full((2), 0.5)

# Now we need to adjust \pi, A, and B to maximize the probability of the observed die rolls
# If that sounds similar our steps to maximize mixture models, good. It is.
# We can use the same expectation-maximization steps.
# As with mixture models, we will adjust our parameters to match the data iteratively.
# With mixture models
# The EM algorithm for Markov models is called the Baum-Welch algorithm
# Estimate the data log likelihood
#   \Sum_{k=1}^K\mathbb{E}[N_k^1]log\pi_k + \Sum_{j=1}^K\Sum_{k=1}^K\mathbb{E}[N_{jk}]logA_{jk} + \Sum_{i=1}^N\Sum_{t=1}^{T_i}\Sum_{k=1}^K p(z_t=k|x_i,\Theta^{old})logp(x_{i,t|\theta_k)
# Those three terms are not so crazy
# \Sum_{k=1}^K\mathbb{E}[N_k^1] = \Sum_{i=1}^N p(z_{i1} = k| x_i, \Theta^{old})
# That's the sum of probabilities that we were in hidden state z_i = k given that we saw x_i.
# In gmm, that was the sum of responsibility of a cluster.
# We'll use this value in a moment to recompute the prior of this hidden state.

# We have all of the data to look at, which means that we can compute the probability of the observation at time i given t=1 to i-1 give the different hidden states.
# Then we go backwards from T to t=i to see if future evidence supports being in each hidden state.

# As usual, all math will be done in log space

# What's the probability of a hidden state being the active one given an observation, x?
# p(z_t = j|x_{1:t}) uses all past information.
# p(z_t = j|x_{1:T}) uses all information. We'll break this into two parts
# p(z_t = j|x_{1:T}) = p(z_t = j|x_{1:t})p(x_{t+1:T}|z_t = j)

# Obviously if we switch from hidden states with high probability, more data doesn't tell us much. We hope for some state stability.
# This is equivalent to saying that mixture models with an infinite number of clusters aren't very helpful.

def forward_pass_log(log_pi, log_A, log_B, observations):
    """Computes the forward probabilities (alpha) in log space.
    This is the first part of the equation below:
      p(z_t = j|x_{1:T}) = p(z_t = j|x_{1:t})p(x_{t+1:T}|z_t = j)
    """
    T = len(observations)
    # The number of hidden states
    N = log_A.shape[0]
    # Initialize with negative infinity (log(0))
    log_alpha = np.full((T, N), -np.inf)

# Initialization (t=0)
    # Begin with the probabilities of seeing the first observation for the emissions of each hidden state
    log_alpha[0, :] = log_pi + log_B[:, observations[0]]

# Termination: Log-likelihood of the observation sequence
    log_likelihood = np.logaddexp.reduce(log_alpha[T-1, :])

return log_alpha, log_likelihood

def backward_pass_log(log_A, log_B, observations):
    """Computes the backward probabilities (beta) in log space.
    This is the second part of the equation below:
      p(z_t = j|x_{1:T}) = p(z_t = j|x_{1:t})p(x_{t+1:T}|z_t = j)
    """
    T = len(observations)
    # The number of hidden states
    N = log_A.shape[0]
    # Initialize with negative infinity (log(0))
    log_beta = np.full((T, N), -np.inf)

# Initialization (t=T-1): P(nothing | S_T-1) = 1, so log(1) = 0
    log_beta[T-1, :] = 0.0

return log_beta

def compute_expectations_log(log_alpha, log_beta, log_A, log_B, observations, log_likelihood):
    """
    Computes gamma (single-state expectation) and xi (two-state transition expectation)
    in log space.
    These are from Murphy 17.101 and 17.102
      gamma is p(z_t = j|x_{i,1:T_i},\Theta)
      xi is p(z_{t-1} = j,z_{t}=k|x_{i,1:T_i},\Theta)
    Keep in mind that T is actually just as long as our history, although we'll
    loop over all values in this function.
    """
    # T is the number of observations, N is the number of hidden states
    T, N = log_alpha.shape

# Log P(O) is the log_likelihood from the forward pass
    log_PO = log_likelihood

# log_xi[t, i, j] = log P(St=i, St+1=j | O)
    log_xi = np.full((T - 1, N, N), -np.inf)

# Divide (or subtract in log space) by the probability of this observation, normalizing this observation
        log_xi[t, :, :] = log_terms - log_PO

return log_gamma, log_xi

# Helper function for numerical stability: log(sum(exp(a)))
def logsumexp(a):
    """Computes log(sum(exp(a))) robustly. Wrapper for np.logaddexp.reduce.
    Sometimes we would compute a sum, such as when building an expectation.
    In log space we cannot sum. This will convert out of log, sum, and convert back,
    e.g. log(exp(x1) + exp(x2) + ... + exp(xN))
    No conversions actually occur, so this is numerically safe even with small values.
    """
    return np.logaddexp.reduce(a)

def maximization_step(log_gamma, log_xi, observations, vocabulary):
    """
    Compute the new HMM parameters (pi_new, A_new, B_new) from the expectations.
    """
    T = len(observations)
    # Shape is T, N, N
    N = log_xi.shape[1]

# 1. Update Initial Probabilities (pi)
    # pi'_i = gamma_0(i)
    # From 17.103 in Murphy
    log_pi = log_gamma[0, :]

# For A, we take the number of expected transitions from j to k,
    # and divide by the expected transitions from j to any other state.
    # This guarantees that each row sums to 1.

# Population the emission matrix. Initialize with 0 probabilities.
    log_B = np.zeros((N, len(vocabulary))) - np.inf

for k in vocabulary:
        # Find indices where observation O_t equals k
        indices_k = (observations == k)

return log_pi, log_A, log_B

def train_hmm(log_A, log_B, log_pi, observations, vocabulary, n_iterations):
    """
    Train a hiddem Markov model with Baum-Welch
    """

# Keep track of the likelihoods during training
    past_log_likelihoods = []

while (len(past_log_likelihoods) < n_iterations):
        ##################
        # Expectation step

# Use the forward-backward algorithm to solve for
        #    p(z_t = j|x_{1:T}) = p(z_t = j|x_{1:t})p(x_{t+1:T}|z_t = j)
        # using our current parameters.
        # This takes O(K^2T), which isn't too bad if our number of hidden states is small.
        # We also need O(KT) space, which can be worse.
        log_alpha, log_likelihood = forward_pass_log(log_pi, log_A, log_B, observations)
        log_beta = backward_pass_log(log_A, log_B, observations)
        past_log_likelihoods.append(log_likelihood)

# See 17.5.2.1 from Murphy
        log_gamma, log_xi = compute_expectations_log(log_alpha, log_beta, log_A, log_B, observations, log_likelihood)

##################
        # Maximization step
        # Maximize the probability of the observations by adjusting our parameters
        log_pi, log_A, log_B = maximization_step(log_gamma, log_xi, observations, vocabulary)

# See if we should stop training
        if len(past_log_likelihoods) > 2:
            if np.abs(past_log_likelihoods[-1] - past_log_likelihoods[-2]) < 10e-5:
                print(f"Converged after {len(past_log_likelihoods)} steps")
                break

return log_pi, log_A, log_B, log_gamma, log_xi, log_alpha, log_beta, past_log_likelihoods

if __name__ == "__main__":
    parser = argparse.ArgumentParser()
    parser.add_argument(
        "test",
        type=str,
        choices=["coin", "die"],
        default="coin",
        help="type of test to run")
    parser.add_argument(
        "--samples",
        required=False,
        type=int,
        default=200,
        help="The number of samples to use.")

args = parser.parse_args()

if args.test == "die":
        actual_states, observations, vocabulary = getDieObservations(args.samples)
    else:
        actual_states, observations, vocabulary = getCoinObservations(args.samples)

# There are two hidden states
    # We'll represent transitions between them with a transition matrix (often called A)
    # We'll populate the values by assuming that all probabilities are equal
    # The rows must sum to 1; these are all the ways we can leave that state.
    transition_matrix = np.full((2,2), 0.5)

# Now the emission matrix, B
    # We would normally begin by counting the possible symbols,
    # but we know that they are numbers from 1 to 6
    # Begin by assuming some random values
    emission_matrix = np.random.rand(transition_matrix.shape[0],len(vocabulary))
    # The rows must sum to 1
    emission_matrix = emission_matrix / emission_matrix.sum(axis=1, keepdims=True)

# Now we need pi, the initial state distribution
    # Let's just assume that all hidden states are equally likely
    pi = np.full((2), 0.5)

# Make log versions of everything
    log_A = np.log(transition_matrix)
    log_B = np.log(emission_matrix)
    log_pi = np.log(pi)

log_pi, log_A, log_B, log_gamma, log_xi, log_alpha, log_beta, past_log_likelihoods = train_hmm(log_A, log_B, log_pi, observations, vocabulary, n_iterations=1000)

A = np.exp(log_A)
    B = np.exp(log_B)

# Figure out which index got the cheating hidden state
    fixed_idx = 1
    if B[fixed_idx,-1] < B[0, -1]:
        fixed_idx = 0
    fair_idx = 1 - fixed_idx

print(f"Transition probabilities for fixed {args.test}: {np.exp(log_A[fixed_idx])}")
    print(f"Transition probabilities for fair {args.test}: {np.exp(log_A[fair_idx])}")
    print(f"Emission probabilities for fixed {args.test}: {np.exp(log_B[fixed_idx])}")
    print(f"Emission probabilities for fair {args.test}: {np.exp(log_B[fair_idx])}")

# Print results
    print(f"Log-Likelihood P(O): {past_log_likelihoods[-1]:.4f}")
    print(f"Total Expected Time in Fair State ({fair_idx}): {np.sum(np.exp(log_gamma[:, fair_idx])):.2f} / {len(observations)} steps")
    print(f"Total Expected Time in Fixed State ({fixed_idx}): {np.sum(np.exp(log_gamma[:, fixed_idx])):.2f} / {len(observations)} steps")

## Check expected transitions (should sum to T-1)
    #expected_total_transitions = np.sum(np.exp(log_xi))
    #print(f"Total Expected Transitions (should be T-1={len(observations)-1}): {expected_total_transitions:.2f}")

## Expected Transitions from Fair
    #expected_fair_transitions = np.sum(np.exp(log_xi[:, fair_idx, :]), axis=0)
    #print(f"\nExpected Transitions FROM Fair (State 0):")
    #print(f"  Fair -> Fair: {expected_fair_transitions[fair_idx]:.2f}")
    #print(f"  Fair -> Fixed: {expected_fair_transitions[fixed_idx]:.2f}")

## Expected Transitions from Fixed
    #expected_fixed_transitions = np.sum(np.exp(log_xi[:, fixed_idx, :]), axis=0)
    #print(f"\nExpected Transitions FROM Fair (State 0):")
    #print(f"  Fixed -> Fixed: {expected_fixed_transitions[fixed_idx]:.2f}")
    #print(f"  Fixed -> Fair: {expected_fixed_transitions[fair_idx]:.2f}")

#print(actual_states)
    likely_state = (log_gamma.T[fixed_idx] > log_gamma.T[fair_idx]).astype(int)
    #print(likely_state)
    accuracy = np.mean(actual_states == likely_state)
    print(f"State prediction accuracy {accuracy}")
```

---

## Predictions

<pre>
$ python3 occassionally_unfair_casino.py coin --samples 300
Converged after 242 steps
Transition probabilities for fixed coin: [0.40704574 0.59295426]
Transition probabilities for fair coin: [0.68086784 0.31913216]
Emission probabilities for fixed coin: [0.00214286 0.99785714]
Emission probabilities for fair coin: [0.53567622 0.46432378]
Log-Likelihood P(O): -181.5306
Total Expected Time in Fair State (0): 167.48 / 300 steps
Total Expected Time in Fixed State (1): 132.52 / 300 steps
State prediction accuracy 0.6733333333333333
</pre>

---

## Reliability

* It is sometimes very close to correct, but sometimes doesn't converge
  * This has to do with the transition probabilities
* Since we are using sequences to improve statistics, there needs to be correlation
  * If the transition matrices all have the same probability, sequences help less

---

## Unfair Die

* Let's do something with more complicated observations

```python
def getDieObservations(k=250):
    # true emission probabilities for a fair die
    fair_emission = [1.0/6] * 6
    # true emission probabilities for a loaded die that rolls a 6 50% of the time
    loaded_emission = [0.1]*5 + [0.5]

actual_emission_matrix = np.array([fair_emission, loaded_emission])

actual_transition_matrix = [ [0.95, 0.05],
                                 [0.05, 0.95]]

# Using the values above, we can generate some data

actual_states = []
    observations = []

cur_state = 0

# Die rolls should be 1-6, but we'll simplify some code and use 0-5
    vocabulary=list(range(6))

---

## Die results

<pre>
$ python3 occassionally_unfair_casino.py die --samples 500
Converged after 95 steps
Transition probabilities for fixed die: [0.14749579 0.85250421]
Transition probabilities for fair die: [0.88183479 0.11816521]
Emission probabilities for fixed die: [0.07026388 0.17342299 0.00586357 0.12337442 0.09226615 0.53480899]
Emission probabilities for fair die: [0.17151531 0.11997255 0.20957302 0.14631032 0.19351346 0.15911533]
Log-Likelihood P(O): -844.1528
Total Expected Time in Fair State (0): 275.24 / 500 steps
Total Expected Time in Fixed State (1): 224.76 / 500 steps
State prediction accuracy 0.752
</pre>

---

## Discussion

* Not only are we predicting the hidden state, we are also predicting A and B
  * And the hidden state has no labels!
* Shortcomings
  * How many hidden states? How should they be initialized?
* But as-is we have something useful

---

## Usefulness

* How is this useful?
* Great for simulations!
  * Used in testing
  * Used in random data generation
  * Used in reinforcement learning

---

## HMM Utility

* HMMs are flexible
  * We can find different uses depending upon our available information
* Many real-world applications are in continuous space
  * But there are plenty in discrete space
  * We'll go over a couple of use-cases next class