# Jim's typical preface
% matplotlib inline
from numpy import *
from matplotlib.pyplot import plot
import matplotlib.pyplot as plt
import numpy as np
Getting started¶
Let's play around with this and see what's what.
First let's draw explicitly the markov states and transitions. There are two hidden states (Rain, not Rain) and two observable states (Umbrella, not Umbrella).
The transition matrix that moves the hidden (Rain, not Rain) states forward is
$$ \left[ \begin{array}{cc} 0.7 & 0.3 \\ 0.3 & 0.7 \\ \end{array} \right] \left[ \begin{array}{c} R_{t} \\ !R_{1} \end{array} \right] = \left[ \begin{array}{c} R_{t+1} \\ !R_{t+1} \end{array} \right] $$and the observation matrix that determines the probabilities of (Umbrella, not Umbrella) given the hidden state is
$$ \left[ \begin{array}{cc} 0.9 & 0.2 \\ 0.1 & 0.8 \\ \end{array} \right] \left[ \begin{array}{c} R_{t} \\ !R_{1} \end{array} \right] = \left[ \begin{array}{c} U_{t} \\ !U_{t+1} \end{array} \right] $$The columns must sum to 1.0. The rows need not. (Explain why.)
Ignoring the observables for a moment, we can plot a possible time evolution.
# First define the transition matrix T
transition = matrix([[0.7, 0.3], [0.3, 0.7]])
# Then set some starting condition (R[0], !R[0]).
# Let's assume that the weather starts sunny, just to see what happens.
r0 = matrix([[0.0], [1.0]])
# Then we can propogate the probabilities forward
# in time by successive multiplications by the matrix.
time_steps = 6
for t in xrange(time_steps):
print transition**t * r0
# Here's the probability of rain over time.
rain = array([ (transition**t * r0)[0,0] for t in range(time_steps) ])
print "P(rain) = {}".format(rain.round(3))
# A plot is always nice.
plt.figure()
ax = plt.subplot(111) # This means plot 1 x 1 plots, number 1.
ax.plot(rain, 'ro') # red (r) circles (o)
plt.xlabel("time")
plt.ylabel("P(rain)")
# ax.set_yscale('log')
# ax.set_xscale('log')
ax.set_xlim(-1, time_steps)
ax.set_ylim(-0.1, 1.1)
# You can use LaTeX in plot annotations ...
# ax.text(10, 1e6, '$f(n) = 10n + 0.1 \, n \, \log(n)$', color='red', fontsize=16)
# ax.text(10, 1e5, '$g(n) = 0.001 \, n^2$', color='blue', fontsize=16)
plt.show()
Observations¶
How does this change if we have observed a specific sequence of umbrellas on days 1 through 5, say U = [True, True, False, False, False] ?
We are given $ P(umbrella=? | rain=?) $. From that we can use Bayes Theorem to find $P(rain=? | umbrella=?)$. If way we see an umbrella at time t, then
$$ P(R_t|U_t) = C P(U_t|R_t) * P(R_t) P(!R_t|U_t) = C P(U_t|!R_t) * P(!R_t) $$where as usual the normalization constant C can be found by requiring that these sum to 1.0.
The conditional probabilities P(U|R) and P(U|!R) are already given. The a priori probability P(R) and P(!R) are what we have without the obvservation.
So for each time step, we weight the probabilities (P(R), P(!R)) that we have by either (P(U|R), P(U|!R)) if we see an umbrella, or by (P(!U|R), P(!U|!R)) if we don't see an umbrella at that time.
This is the FORWARD algorithm.
# This time I'll assume that the starting a priori is P(R)=P(!R)=0.5
observations = ['', 'U', 'U', 'U', '!U', '!U', '!U'] # day 0 is our starting point.
rain = [0.5, None, None, None, None, None] # P(R), what we know so far.
r = matrix([[rain[0]], [1.0-rain[0]]]) # initial hidden state probabilities
for t in range(1, time_steps):
# First propogate the hidden states forward, without the observations.
r = transition * r # (Note that this is a matrix equation.)
(p_R, p_not_R) = map(float, list(r)) # extract the two probabilities
# Now weight P(R) and P(!R) by what we know for this time.
if observations[t] == 'U':
(p_R, p_not_R) = (0.9 * p_R, 0.2 * p_not_R) # weight by P(U|R), P(U|!R)
else:
(p_R, p_not_R) = (0.1 * p_R, 0.8 * p_not_R) # weight by P(!U|R), P(!U|!R)
r[0,0] = p_R/(p_R + p_not_R) # normalize so they sum to 1.0
r[1,0] = p_not_R/(p_R + p_not_R)
rain[t] = r[0,0] # remember this value
rain = array(rain)
print "P(rain) = {}".format(rain.round(3))
# A plot is always nice.
plt.figure()
ax = plt.subplot(111) # This means plot 1 x 1 plots, number 1.
ax.plot(rain, 'ro') # red (r) circles (o)
plt.xlabel("time")
plt.ylabel("P(rain)")
# ax.set_yscale('log')
# ax.set_xscale('log')
ax.set_xlim(-1, time_steps)
ax.set_ylim(-0.1, 1.1)
# You can use LaTeX in plot annotations ...
# ax.text(10, 1e6, '$f(n) = 10n + 0.1 \, n \, \log(n)$', color='red', fontsize=16)
# ax.text(10, 1e5, '$g(n) = 0.001 \, n^2$', color='blue', fontsize=16)
plt.show()
The result is intuitive:
- The starting probabilitiy is assumed to be 0.5
- Then for three days we observe umbrella, and p(rain) increases towards 1.0
- Then for two days we see no umbrella, and p(rain) decreases towards 0.0.
Each of these points is the probability given the previous evidence only - we're not trying to look back after knowing all the observations to see what the early state should have been. (The FORWARD-BACKWARD algorithm does that.)
Next¶
So ... what next to understand the material in this chapter?
Repeat this "filter" calculation for another similar problem, such as the Dow Jones example.
Practice the matrix stuff if that is new to you.
Determine which of the Dow Jones models is a better fit to the data ... and discuss how you might go about making that determination.
Learn about the other variations of this idea: FORWARD-BACKWARD and VITERBI. What exactly do those implement?
Use a similar approach to the FORWARD algorith to implement a "particle filter" algorithm for robot location.
For the robot problem, the state space (all the poossible robot locations) is typically too large for a practical version of the procedure we just followed. So instead we randomly choose some number of "particles", i.e. virtual robots. Then we either add or remove these virtual robots (resample) depending on the observations, in the same way that the observations weighted the state probabilities in the FORWARD algorithm.
More on that to come ...