# CS 530 - Lecture 06

## More Discrete Space Predictions

Bernhard Firner

2026-02-10

---

## Review

* If you want something with more in-depth details on Markov Models
  * Machine Learning: A Probabilistic Perspective
* We'll be going over 14.3.2 from the book today
* But first, let's wrap up some important concepts

---

## Hidden Markov Models

* HMMs are useful for discrete-space predictions
  * In continuous spaces, they are harder
  * Remember that some algorithms has $O(k^2)$ costs, where $k$ is the number of hidden states
* But there are a few points worth lingering upon

---

## HMM Uses

* Filtering
* Smoothing
* Prediction
* And stationary distributions
  * Which we'll talk about in a moment

---

## Filtering

* Get an estimate of the hidden state at time $T$, given observations $X_{t-n}$ ... $X_{t-1}$
* We used the Forward Algorithm to calculate the probability of any particular chain
* And if we want the most likely states along a chain, we used the Viterbi algorithm

---

## Viterbi Vs Forward-Backward

* When learning HMM parameters, we used forward-backward to estimate $\gamma$
  * $\gamma$ was the hidden state probabilities at each time
* How are Viterbi and Forward-Backward different?

---

## MAP vs MPE

* The Viterbi algorithm computes the *maximum a posteriori* (MAP)
  $z^* = \underset{z_{1:T}}{\mathrm{argmax}}(z_{1:T}|x_{t:T})$
* $\gamma$ is an estimate that maximizes the posterior marginals (MPM):
  $\hat{z} = (\underset{z_1}{\mathrm{argmax}}p(z_1|x_{1:T}), ..., (\underset{z_T}{\mathrm{argmax}}p(z_T|x_{1:T}))$

---

## Simplified

* Viterbi finds the single most likely sequence of events
* $\gamma$ is the maximized likelihood at any individual time step
* Consider the sentence "recognize speech"
  * That may be the most likely total chain given the phonemes you hear
  * But "wreck a nice beach" could contain the individual most likely words

---

## Another Example

* Let's do something useful, and compare the two approaches
* We will localize a robot based upon some basic sensing
  * Perhaps this is a cleaning robot with short-range ultrasonic sensors

---

## Localization

* We can represent every location as a hidden state
  * And the things that can be sensed will be the emission matrix entries
* In this case, we will sense walls in four directions with ultrasonic sensors

---

## Metrics

* We could care about two different things
* Localization error
  * How close is the estimated final location?
* Path error
  * How wrong is the path, on average?
* Will use Manhattan distances

---

## Errors

* Let's add in errors
  * Independent for each sensor
* Flip 0's to 1's and vice-versa
* This will show how robust the algorithms are to noise

-v-

## HMM Code

```python
##
# Implementation of hidden markov model algorithms.
# See Machine Learning: A Probabilistic Perspective by Kevin P. Murphy, 2012
# for more details on the equations and notation.
# References to specific equations in the book are provided in comments.
##

import numpy
import random
import functools

# Helper function for numerical stability: log(sum(exp(a)))
# This is only here as a comment and explanation, the function already exists in numpy.
def logsumexp(a):
    """Computes log(sum(exp(a))) robustly. Wrapper for numpy.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 numpy.logaddexp.reduce(a)

def multistate_log_prob(log_B, neg_log_B, O):
    """
    This implementation assumes multiple states in the observation, with a
    vector of 0s and 1s in the observations.
    The probability of a 1 at each observation for different states is B[:, i].
    The probability of the observation at time t is (O[t,i] * B[:,i]) + ((1 - O[t,i]) * (1 - B[:,i]))
    The probability over all of the states, assuming that they are uncorrelated,
    is numpy.prod((O[t,i] * B[:,i]) + ((1 - O[t,i]) * (1 - B[:,i])), axis=1)
    We are doing everything in log space, so that turns into
    numpy.sum(numpy.logaddexp.reduce((O * log_B[:,i]) + ((1 - O) * (neg_log_B[:,i]))), axis=1)
    """
    # Continue using nan_to_num to prevent -inf values that turn into nans with more math
    # Any conversion back to non-log probabilities will still yield 0.
    return numpy.nan_to_num(numpy.sum((O * log_B) + ((1 - O) * (neg_log_B)), axis=1))

def log_prob(log_B, O):
    """
    Returns the probability of seeing the given observation with the provided
    emission probabilities.
    This works with a single observation at each unit time and is simply:
        log_B[:, O]
    where O is the index of the observation.
    """
    return log_B[:, observations[0]]

def viterbi_log(log_pi, log_A, log_B, observations):
    """
    Compute the most probably sequence of states:
        z^* = argmax_z_{1|T} p(z_{1:T}|x_{1:T})
        See Murphy eqn 17.68
    This is the MAP (maximum a posteriori), the jointly most probable sequence of states.
    This is the most likely chain, taken as a whole, and is always globally consistent.
    The MPE (maximizer of the posterior marginals) obtained from
    forwards-backwards is the most likely element at each time, but is not
    necessarily globally consistent.
    The example Murphy gives is "recognize speech" vs "wreck a nice beach"
    """
    # Set up the proper function to get the probabilities of observations
    observations = numpy.array(observations)
    if (1 == observations.ndim):
        observe_prob = functools.partial(log_prob, log_B, _)
    else:
        # We didn't require this as an argument to remain consistent with
        # forward_pass_log but we still need it.
        neg_log_B = numpy.nan_to_num(numpy.log(1.0 - numpy.exp(log_B)))
        observe_prob = functools.partial(multistate_log_prob, log_B, neg_log_B)

# Initialize with the initial state probabilities
    max_probs = [log_pi]

# The number of hidden states
    N = log_A.shape[0]

# Now we find the maximally likely states at each time T given the observations
    # z^* = \underset{z_{1|T}}{\mathrm{argmax}} p(z_{1:T}|x_{1:T})
    for t in range(0, len(observations)):
        # We find the probability of transitioning into each state and seeing
        # an observation, then pick the max

# That's p(z_t = j | x_t)
        log_B_t = observe_prob(observations[t])

max_prob = []

# Loop over the hidden states
        for j in range(N):
            # Get the probabilities of transitioning into this state
            # Continue using nan_to_num to prevent -inf values that turn into nans with more math
            # Any conversion back to non-log probabilities will still yield 0.
            prob_transition = numpy.nan_to_num(max_probs[-1] + log_A[j])
            # One of these was the maximally likely path into this state
            max_j = numpy.argmax(prob_transition)
            # Need to multiply (add in log space) the probability of the observation
            max_prob.append(max_j + log_B_t[j])

# Add the most likely probabilities of the next states onto the chain
        max_probs.append(max_prob)

# Return the most likely sequence of states and the probabilities of each
    # hidden state at each time.
    return numpy.argmax(max_probs, axis=1), max_probs

def forward_pass_log(log_pi, log_A, log_B, observations):
    """Computes the forward state 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)
    """
    observations = numpy.array(observations)
    if (1 == observations.ndim):
        observe_prob = functools.partial(log_prob, log_B, _)
    else:
        # We didn't require this as an argument to remain consistent with
        # forward_pass_log but we still need it.
        neg_log_B = numpy.nan_to_num(numpy.log(1.0 - numpy.exp(log_B)))
        observe_prob = functools.partial(multistate_log_prob, log_B, neg_log_B)
    # The observations are being used to select, to they remain 0s and 1s
    T = len(observations)
    # The number of hidden states
    N = log_A.shape[0]
    # Initialize with negative infinity (log(0))
    log_alpha = numpy.full((T, N), -numpy.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 + observe_prob(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 = observe_prob(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 numpy.logaddexp.reduce for numerical stability)
            log_sum = numpy.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 = numpy.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)
    """
    observations = numpy.array(observations)
    if (1 == observations.ndim):
        observe_prob = functools.partial(log_prob, log_B, _)
    else:
        # We didn't require this as an argument to remain consistent with
        # forward_pass_log but we still need it.
        neg_log_B = numpy.nan_to_num(numpy.log(1.0 - numpy.exp(log_B)))
        observe_prob = functools.partial(multistate_log_prob, log_B, neg_log_B)
    T = len(observations)
    # The number of hidden states
    N = log_A.shape[0]
    # Initialize with negative infinity (log(0))
    log_beta = numpy.full((T, N), -numpy.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 = observe_prob(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] = numpy.logaddexp.reduce(log_sum)

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.
    """
    observations = numpy.array(observations)
    if (1 == observations.ndim):
        observe_prob = functools.partial(log_prob, log_B, _)
    else:
        # We didn't require this as an argument to remain consistent with
        # forward_pass_log but we still need it.
        neg_log_B = numpy.nan_to_num(numpy.log(1.0 - numpy.exp(log_B)))
        observe_prob = functools.partial(multistate_log_prob, log_B, neg_log_B)
    # 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 = numpy.full((T - 1, N, N), -numpy.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 = observe_prob(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, :, numpy.newaxis] +   # log alpha_t(n),    shape (N, 1)
            log_A +                         # log A_ij,          shape (N, N)
            log_B_t_plus_1[numpy.newaxis, :] + # log B_j(x_{t+1}),  shape (1, N)
            log_beta[t+1, numpy.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

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 = numpy.logaddexp.reduce(log_xi)
    log_gamma_sum = numpy.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]
    # numpy.newaxis ensures the denominator broadcasts correctly across columns j
    # From 17.103 in Murphy
    log_A = log_xi_sum - log_gamma_sum[:, numpy.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 = numpy.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 = numpy.zeros((N, len(vocabulary))) - numpy.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 = numpy.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] = -numpy.inf
            # We already initialized these to 0
            pass

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 numpy.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

def generate_random_observations(A, B, pi, length):
    """
    Generate a numpy array of random observations based upon the provided information.
    Arguments:
        A (numpy.array[[float]]): The transition matrix.
        B (numpy.array[[float]]): The emission matrix.
        pi  (numpy.array[float]): The initial state probabilities.
        length             (int): The length of the random observations.
    Returns:
        path ([int]): List of traversed states, sized length+1. The first state has no observations.
        observations ([[int]]): List of observations, 0 or 1 per observations.
    """

num_states = A.shape[0]
    vocab_len = B.shape[1]
    begin = random.choices(numpy.arange(num_states), weights=pi)[0]

path = [begin]
    observations = []

# Keep drawing new states and observations until we generate all that were requested
    while len(observations) < length:
        next_step = random.choices(numpy.arange(num_states), weights=A[path[-1]])[0]
        path.append(next_step)
        new_observations = []
        for oidx in range(vocab_len):
            p = B[path[-1],oidx]
            if random.random() < p:
                new_observations.append(1)
            else:
                new_observations.append(0)
        observations.append(new_observations)

return path, observations
```

-v-

## Localization Code

```python
#! /usr/bin/python3

##
# Use the Viterbi algorithm to select the most probable sequence of locations
# given a sequence of readings of adjacent spaces, indicating whether or not
# they are traversable.
# The user must also provide a map, in matrix form, that indicates the
# locations. Locations are treated as hidden states, and the observations are
# whether or not another state is accessible in each of the cardinal
# directions.
#
# There are both discrete and continuous versions of this problem.
##

import hmm_lib as hmm
import numpy

import argparse
import math

if __name__ == "__main__":
    parser = argparse.ArgumentParser()
    parser.add_argument(
        "map",
        type=str,
        help="A txt file with 0s and 1s in a grid representing a rectangular map. 0s are traversable, 1s are barriers.")
    parser.add_argument(
        "--walk_length",
        required=False,
        type=int,
        default=50,
        help="The length of the random walk through the space.")
    parser.add_argument(
        "--trials",
        required=False,
        type=int,
        default=1,
        help="The number of trials to run. Gathers statistics.")
    parser.add_argument(
        "--error_rate",
        required=False,
        type=float,
        default=0.0,
        help="The error rate for adjacent traversable state detection.")
    parser.add_argument(
        "--allow_stay",
        required=False,
        action='store_true',
        default=False,
        help="Allow the agent to remain in the same state without moving.")
    parser.add_argument(
        "--algorithm",
        required=False,
        choices=['map', 'mpe'],
        default='map',
        help="Use either the viterbi algorithm (MAP) or estimate with forward-backward (MPE).")

args = parser.parse_args()

## Read the map file

occupancy = numpy.genfromtxt(args.map, delimiter=1, dtype=int)

## Convert the map into a set of hidden states and emission probabilities
    ## We are making multiple independent observations, so the probabilities can be multiplied.

nhidden = numpy.count_nonzero(occupancy == 0)

# The transition matrix
    A = numpy.zeros((nhidden, nhidden))
    # The emission matrix. There are four possible sensing directions.
    B = numpy.zeros((nhidden, 4))

# Also remember the locations of the hidden states for later plotting
    hidden_to_location = []
    # And remember how to go from locations back to the hidden states
    location_to_hidden = [[-1] * occupancy.shape[1] for _ in range(occupancy.shape[0])]

hidden_idx = 0
    for row in range(occupancy.shape[0]):
        for column in range(occupancy.shape[1]):
            # The hidden state index
            h_idx = len(hidden_to_location)
            # If this is a hidden state, then add it into A and B
            # We can fill in the occupancy now.
            # Transitions will wait until we have the index for each location
            if occupancy[row, column] == 0:
                # Make the emission matrix, which has our expected sensor readings at this location
                # Instead of being binomial, we'll treat this as multinomial, so the total probabilities don't need to sum to 1.
                t_idx = 0
                for y, x in [(row, column-1), (row-1, column), (row, column+1), (row+1, column)]:
                    if y >= 0 and y < occupancy.shape[0] and x >= 0 and x < occupancy.shape[1] and occupancy[y][x] == 0:
                        B[h_idx, t_idx] = 1.0 - args.error_rate
                    else:
                        B[h_idx, t_idx] = args.error_rate
                    t_idx += 1

# Remember the mapping from location to the hidden index and back
                location_to_hidden[row][column] = h_idx
                hidden_to_location.append((row, column))

# Now that all of the hidden states have indices, let's fill in the transition matrix
    for ai in range(nhidden):
        row, column = hidden_to_location[ai]
        # We'll allow a transition to the current location itself if the user permits it.
        # The exact value will be scaled by the number of valid transitions at the end of the loop
        if args.allow_stay:
            A[ai,ai] = 1
        for y, x in [(row, column-1), (row-1, column), (row, column+1), (row+1, column)]:
            if y >= 0 and y < occupancy.shape[0] and x >= 0 and x < occupancy.shape[1] and occupancy[y][x] == 0:
                other_idx = location_to_hidden[y][x]
                A[ai,other_idx] = 1
        # Transitions should sum to 1
        A[ai] = A[ai] / numpy.sum(A[ai])

# pi, the initial state distribution
    # Make all hidden states equally likely
    pi = numpy.full((nhidden), 1.0/nhidden)

accuracies = []
    average_path_errors = []
    localization_errors = []
    for trial in range(args.trials):

# For now, we'll support generating a random path and then localizing with the HMM.
        # TODO We could do more things with the model.
        true_path, observations = hmm.generate_random_observations(A, B, pi, args.walk_length)

# Make log versions of everything
        # Use nan_to_num to convert -inf values to very negative finite numbers instead.
        # This prevents some unexpected headaches, such as -inf * 0 = nan rather than 0.
        log_A = numpy.nan_to_num(numpy.log(A))
        log_B = numpy.nan_to_num(numpy.log(B))
        log_pi = numpy.nan_to_num(numpy.log(pi))

if args.algorithm == 'map':
            # Use only the observations to reconstruct the path
            m_a_posteriori, all_state_probs = hmm.viterbi_log(log_pi, log_A, log_B, observations)
            estimated_path = m_a_posteriori

if args.algorithm == 'mpe':
            # Compute MPE
            log_alpha, log_likelihood = hmm.forward_pass_log(log_pi, log_A, log_B, observations)
            log_beta = hmm.backward_pass_log(log_A, log_B, observations)
            log_gamma, log_xi = hmm.compute_expectations_log(log_alpha, log_beta, log_A, log_B, observations, log_likelihood)
            # Notice that this path may not be consistent
            estimated_path = numpy.argmax(log_gamma, axis=1)
            # We need to pad by 1
            estimated_path = estimated_path
            # This method doesn't estimate the initial state
            true_path = true_path[1:]

manhatten_distances = []
        for actual, estimated in zip(true_path, estimated_path):
            actual_location = hidden_to_location[actual]
            estimated_location = hidden_to_location[estimated]
            distance = math.fabs(actual_location[0] - estimated_location[0]) + math.fabs(actual_location[1] - estimated_location[1])
            manhatten_distances.append(distance)

average_path_error = numpy.mean(manhatten_distances)
        localization_error = manhatten_distances[-1]

# These are rather long to look at.
        #print(f"Actual path: {true_path}")
        #print(f"Predicted path: {estimated_path}")

accuracy = numpy.mean(true_path == estimated_path)
        print(f"State prediction accuracy {accuracy}")
        # Print the localization error, as the Manhattan distance from the true location
        print(f"Localization error {localization_error}")
        # Print the path error, the average Manhattan distance of states on the path from the true path
        print(f"Average path error {average_path_error}")

accuracies.append(accuracy)
        localization_errors.append(localization_error)
        average_path_errors.append(average_path_error)

print("Results:")
    print(f"algorithm, walk length, trials, error rate, accuracy, localization error, average path error")
    print(f"{args.algorithm}, {args.walk_length}, {args.trials}, {args.error_rate}, {numpy.mean(accuracies)}, {numpy.mean(localization_errors)}, {numpy.mean(average_path_errors)}")
```

---

## Localization Error

</div>
<div class="col">

</div>
</div>

---

## Path Error

</div>
<div class="col">

</div>
</div>

---

## Path Accuracy

</div>
<div class="col">

</div>
</div>

---

## Fixed-Lag Smoothing

* We can use MAP in an online system
  * What about MPE? We need future events, so we can't, right?
* We could accept a lag in our estimate
* Wait for additional sensor readings, then find the MPE estimate
* This is called fixed-lag smoothing

---

## Limitations

* Hidden Markov models are great for discrete state problems
  * And for things that can be binned into discrete spaces
* For example, we could model the emissions as parameters of distrbutions rather than categorical outputs

---

## Speech Modelling

* HMMs are often used in speed processing
  * Hidden states could be phonemes
  * There are correlations between them, which means a transition matrix
* The emissions will be probablistic gaussians rather than categorical

---

## Discretization

* We could also discretize a continuous space
  * For example, to localize in a continuous space
* Remember, the cost of HMM algorithms grows with the square of the hidden states
  * The more we quantize a continuous state space, the more expensive this becomes
* Let's go over one last technique with Markov Models before we move on!

---

## Stationary Distributions

* $\pi_i$ is the distribution of state probabilities at time $i$
* If certain conditions are met, there is a special distribution:
  * $\pi = \pi A$
* Meaning that this state, when multiplied by the transition matrix, remains the same
* Only exists if all states are reachable and the states are not time dependent

---

## What is it?

* The stationary distribution tells us the state probabilities over infinite time
  * We could attempt to simulate these, but solving gives us precise values
* Solving is tricky by hand, but can be done with numerical methods

---

## Solving

* Turn the HMM into a system of linear equations
* One equation is from statistics:
  $\sum_{i=1}^{k}\pi(i) = 1$
  * The sum of all state probabilities must be 1

---

## Other Equations

* Expand each column of $\pi = \pi A$ to get an equation
  * e.g. Let's say state Z is transitioned into from states X and Y with probability 1
  * $\pi(Z) = \pi(X) + \pi(Y)$
* Meaning that we should see state Z twice as much as X and Y

---

## Surprising Results

* Let's say you are observing random digits
* Which sequence are you more likely to observe first
  * 111
  * 123
  * 121

-v-

## Brute Force

```python
import math
import random
import numpy

times = []

for _ in range(100000):
    # Time to get to 111
    vocab = [1, 2, 3]
    count = 0
    target = [1, 1, 1]
    time = 0
    while count < len(target):
        v = random.choice(vocab)
        if v == target[count]:
            count += 1
        else:
            count = 0
        time += 1

times.append(time)

print(f"Took {numpy.mean(times)} steps for 111")

times = []

for _ in range(100000):
    # Time to get to 111
    vocab = [1, 2, 3]
    count = 0
    target = [1, 2, 1]
    time = 0
    while count < len(target):
        v = random.choice(vocab)
        if v == target[count]:
            count += 1
        elif v == target[0]:
            # Recovery
            count = 1
        else:
            count = 0
        time += 1

times.append(time)

print(f"Took {numpy.mean(times)} steps for 121")

times = []

for _ in range(100000):
    # Time to get to 111
    vocab = [1, 2, 3]
    count = 0
    target = [1, 2, 3]
    time = 0
    while count < len(target):
        v = random.choice(vocab)
        if v == target[count]:
            count += 1
        elif v == target[0]:
            # Recovery
            count = 1
        else:
            count = 0
        time += 1

times.append(time)

print(f"Took {numpy.mean(times)} steps for 123")
```

---

## 111 Transitions

<p><math display="block" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><mi>A</mi><mo>=</mo><mtable><mtr><mtd columnalign="center" style="text-align: center"><mi>a</mi></mtd></mtr><mtr><mtd columnalign="center" style="text-align: center"><mi>b</mi></mtd></mtr><mtr><mtd columnalign="center" style="text-align: center"><mi>c</mi></mtd></mtr><mtr><mtd columnalign="center" style="text-align: center"><mi>d</mi></mtd></mtr><mtr><mtd columnalign="center" style="text-align: center"><mspace width="1.0em"></mspace></mtd></mtr></mtable><mrow><mo stretchy="true" form="prefix">(</mo><mtable><mtr><mtd columnalign="center" style="text-align: center"><mfrac><mn>2</mn><mn>3</mn></mfrac></mtd><mtd columnalign="center" style="text-align: center"><mfrac><mn>1</mn><mn>3</mn></mfrac></mtd><mtd columnalign="center" style="text-align: center"><mn>0</mn><mi>.</mi></mtd><mtd columnalign="center" style="text-align: center"><mn>0</mn><mi>.</mi></mtd></mtr><mtr><mtd columnalign="center" style="text-align: center"><mfrac><mn>2</mn><mn>3</mn></mfrac></mtd><mtd columnalign="center" style="text-align: center"><mn>0</mn></mtd><mtd columnalign="center" style="text-align: center"><mfrac><mn>1</mn><mn>3</mn></mfrac></mtd><mtd columnalign="center" style="text-align: center"><mn>0</mn><mi>.</mi></mtd></mtr><mtr><mtd columnalign="center" style="text-align: center"><mfrac><mn>2</mn><mn>3</mn></mfrac></mtd><mtd columnalign="center" style="text-align: center"><mn>0</mn></mtd><mtd columnalign="center" style="text-align: center"><mn>0</mn></mtd><mtd columnalign="center" style="text-align: center"><mfrac><mn>1</mn><mn>3</mn></mfrac></mtd></mtr><mtr><mtd columnalign="center" style="text-align: center"><mn>1</mn></mtd><mtd columnalign="center" style="text-align: center"><mn>0</mn></mtd><mtd columnalign="center" style="text-align: center"><mn>0</mn></mtd><mtd columnalign="center" style="text-align: center"><mn>0</mn></mtd></mtr></mtable><mo stretchy="true" form="postfix">)</mo></mrow></mrow><annotation encoding="application/x-tex">
A= 
\begin{matrix}
a \\
b \\
c \\
d \\
\quad
\end{matrix}
\begin{pmatrix}
 \frac{2}{3} &   \frac{1}{3}  & 0.   & 0. \\
 \frac{2}{3} &   0   & \frac{1}{3}  & 0. \\
 \frac{2}{3} &   0   & 0   & \frac{1}{3} \\
 1 &   0   & 0   & 0
\end{pmatrix}</annotation></semantics></math></p>

---

## 111 Equations

* $a + b + c + d = 1$
* $b = \frac{1}{3}a$
* $c = \frac{1}{3}b$
* $d = \frac{1}{3}c$
* $a = \frac{2}{3}a + \frac{2}{3}b + \frac{2}{3}c + d$

---

## Solving 111

* $a + \frac{1}{3}a + \frac{1}{9}a + \frac{1}{27}a = 1$
* $a = \frac{27}{40}$
* $b = \frac{9}{40}$, $c = \frac{3}{40}$, $d = \frac{1}{40}$
* So we expect 40 steps to arrive at state d
  * However, subtract one because our state trasitions include restarting

---

## 121 Transitions

<p><math display="block" xmlns="http://www.w3.org/1998/Math/MathML"><semantics><mrow><mi>A</mi><mo>=</mo><mtable><mtr><mtd columnalign="center" style="text-align: center"><mi>a</mi></mtd></mtr><mtr><mtd columnalign="center" style="text-align: center"><mi>b</mi></mtd></mtr><mtr><mtd columnalign="center" style="text-align: center"><mi>c</mi></mtd></mtr><mtr><mtd columnalign="center" style="text-align: center"><mi>d</mi></mtd></mtr><mtr><mtd columnalign="center" style="text-align: center"><mspace width="1.0em"></mspace></mtd></mtr></mtable><mrow><mo stretchy="true" form="prefix">(</mo><mtable><mtr><mtd columnalign="center" style="text-align: center"><mfrac><mn>2</mn><mn>3</mn></mfrac></mtd><mtd columnalign="center" style="text-align: center"><mfrac><mn>1</mn><mn>3</mn></mfrac></mtd><mtd columnalign="center" style="text-align: center"><mn>0</mn><mi>.</mi></mtd><mtd columnalign="center" style="text-align: center"><mn>0</mn><mi>.</mi></mtd></mtr><mtr><mtd columnalign="center" style="text-align: center"><mfrac><mn>1</mn><mn>3</mn></mfrac></mtd><mtd columnalign="center" style="text-align: center"><mfrac><mn>1</mn><mn>3</mn></mfrac></mtd><mtd columnalign="center" style="text-align: center"><mfrac><mn>1</mn><mn>3</mn></mfrac></mtd><mtd columnalign="center" style="text-align: center"><mn>0</mn><mi>.</mi></mtd></mtr><mtr><mtd columnalign="center" style="text-align: center"><mfrac><mn>2</mn><mn>3</mn></mfrac></mtd><mtd columnalign="center" style="text-align: center"><mn>0</mn></mtd><mtd columnalign="center" style="text-align: center"><mn>0</mn></mtd><mtd columnalign="center" style="text-align: center"><mfrac><mn>1</mn><mn>3</mn></mfrac></mtd></mtr><mtr><mtd columnalign="center" style="text-align: center"><mn>1</mn></mtd><mtd columnalign="center" style="text-align: center"><mn>0</mn></mtd><mtd columnalign="center" style="text-align: center"><mn>0</mn></mtd><mtd columnalign="center" style="text-align: center"><mn>0</mn></mtd></mtr></mtable><mo stretchy="true" form="postfix">)</mo></mrow></mrow><annotation encoding="application/x-tex">
A= 
\begin{matrix}
a \\
b \\
c \\
d \\
\quad
\end{matrix}
\begin{pmatrix}
 \frac{2}{3} &   \frac{1}{3}  & 0.   & 0. \\
 \frac{1}{3} &   \frac{1}{3}   & \frac{1}{3}  & 0. \\
 \frac{2}{3} &   0   & 0   & \frac{1}{3} \\
 1 &   0   & 0   & 0
\end{pmatrix}</annotation></semantics></math></p>

---

## What should you remember?

* A forward pass can filter, predicting the most likely state
* Smoothing can be done with a forward-backward algorithm
* And the parameters of unknown states can be estimated with Baum-Welch
  * This is a version of Expectation-Maximization
* A stationary distribution can reveal exact solutions to a process

---

## Continuous Space

* If we want to filter when our state space is continuous, we have other options
  * They will be called "state space models"
  * The name doesn't make a distinction between continuous or discrete states clear

---

## Filters with Assumptions

* Kalman Filters
  * All noise must be normal
  * There are approximations if some assumptions don't hold: extended or unscented Kalman filters

---

## Heuristic Approaches

* Particle Filters
  * Stochastic approximation algorithm
  * Can almost always be applied
  * Needs more tuning because of that

---

## Why Heuristics?

* Approximation algorithms generally let us trade off resources for accuracy
* The cost of forward-backward in HMMs is $O(K^2T)$ steps and $O(KT)$ space
  * These are pretty bad when our state space is large
* As with many approximation algorithms, we can tune how many resources we use in the algorithm

---

## Continuous Problem

* So what is an example of a continuous space problem?
* Intertial Measurement IMU
  * Gives frequent measurments of acceleration, rotation, and (usually) megnetometer
* Inertial navigation is a very common use case, so let's talk about IMUs

---

## Acceleration

* Acceleration is usually $m/s^2$
* Gravity is usually present, so you need to remove it
* Ideally, we can derive velocity by multiplying by time
  * But remember, we are sampling
  * We can only estimate the true value

---

## Acceleration Troubles

<pre>
UncalAccel,utcTimeMillis,elapsedRealtimeNanos,UncalAccelXMps2,UncalAccelYMps2,UncalAccelZMps2,BiasXMps2,BiasYMps2,BiasZMps2,CalibrationAccuracy
UncalAccel,1770758949098,2298392917820300,0.05746084,1.1971009,9.376891,-0.05231083,0.24177434,-0.031528257,3
UncalAccel,1770758949101,2298392920331914,0.1436521,1.0295067,9.376891,-0.05231083,0.24177434,-0.031528257,3
UncalAccel,1770758949103,2298392922844206,0.19632454,0.8714894,9.38168,-0.05231083,0.24177434,-0.031528257,3
UncalAccel,1770758949106,2298392925357018,0.20590134,0.7517793,9.400833,-0.05231083,0.24177434,-0.031528257,3
</pre>

---

## Gravity

* Gravity is $9.8m/s^2$
  * But we don't see that in the data!
* Our sensor is never perfectly flat
  * We are going to have to continually guess which components are down and remove them

---

## Rotation

* Gyroscope data is in radians/second
* Ideally, these measurements are 0 when you are still
  * But they have a tendency to drift
* If we are estimating a heading, a tiny error leads to huge long-term errors

<pre>
# UncalGyro,utcTimeMillis,elapsedRealtimeNanos,UncalGyroXRadPerSec,UncalGyroYRadPerSec,UncalGyroZRadPerSec,DriftXRadPerSec,DriftYRadPerSec,DriftZRadPerSec,CalibrationAccuracy
UncalAccel,1770758949098,2298392917820300,0.05746084,1.1971009,9.376891,-0.05231083,0.24177434,-0.031528257,3
UncalAccel,1770758949101,2298392920331914,0.1436521,1.0295067,9.376891,-0.05231083,0.24177434,-0.031528257,3
UncalAccel,1770758949103,2298392922844206,0.19632454,0.8714894,9.38168,-0.05231083,0.24177434,-0.031528257,3
UncalAccel,1770758949106,2298392925357018,0.20590134,0.7517793,9.400833,-0.05231083,0.24177434,-0.031528257,3
</pre>

---

## Egomotion and Navigation

* Most navigation requires a rough estimate of movement
* Thus, an agent needs to evaluation its own motion (egomotion)
* This isn't hard over short distances, but errors accumulate quickly
  * So this can end up being one of the more difficult problems in robotics

---

## Next Topics

* We'll go into filtering in continuous spaces
* That will lead us into navigation
* I'll also post the first homework (collecting and using IMU data)
* Also post the second homework, which is to organize into teams and choose a topic