import numpy as np
import pandas as pd
On a chessboard, a knight may make (when unconstrained by edges) 8 different moves from his current position, or he can choose to remain where he is. We'll number the moves:
and if the move is '0', remain in place
We'll count positions as 0 through 7 in each direction, starting in the lower left.
For this notebook, we'll assume that each move is equally likely. Later we may try to infer the likelihood of each
In a dynamic model with unobserved state variables, if we want to infer the states, and have a good model for the system, all the known values at all the timesteps influence the probability of the unknown values at each of the timesteps.
Time to learn about hidden markov models
def move(position, move):
#calculate new position
if move == 0: newpos = position
elif move == 1: newpos = position + [1,2]
elif move == 2: newpos = position + [2,1]
elif move == 3: newpos = position + [2,-1]
elif move == 4: newpos = position + [1,-2]
elif move == 5: newpos = position + [-1,-2]
elif move == 6: newpos = position + [-2,-1]
elif move == 7: newpos = position + [-2,1]
elif move == 8: newpos = position + [-1,2]
if newpos.max()<=7 and newpos.min()>=0: return newpos
else: return position
positionlist = pd.DataFrame(index=range(20000), columns=['x','y'])
positionlist.ix[0] = [3,5]
for i in range(len(positionlist)-1):
positionlist.ix[i+1] = move(positionlist.ix[i].values, np.random.randint(0,9))
# I'd like to change the random selection here, as soon as I figure out what's happened to numpy choose.
plt.figure(figsize(4,4))
plt.plot(positionlist['x'], positionlist['y'], 'o-')
plt.xlim(-.5,7.5)
plt.ylim(-.5,7.5)
plt.title("Knight's path")
print 'x values:'
print positionlist['x'].values
x values: [3 1 3 ..., 7 7 6]
positionlist.head(20).T
0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
x | 3 | 2 | 2 | 4 | 4 | 4 | 4 | 4 | 5 | 3 | 1 | 1 | 1 | 1 | 0 | 1 | 1 | 3 | 4 | 2 |
y | 5 | 7 | 7 | 6 | 6 | 6 | 6 | 6 | 4 | 5 | 6 | 6 | 6 | 6 | 4 | 2 | 2 | 3 | 1 | 0 |
Given this list of x coordinates and the above rule for how the knight can move, can we infer the y coordinates, and the distriution of the random move input? How does our accuracy change with the number of samples? What if we only sample every other turn?
Docs: http://scikit-learn.org/stable/modules/generated/sklearn.hmm.MultinomialHMM.html
from sklearn import hmm
def transprob(x1,y1, x2,y2):
count = 0.
for potmove in range(9):
[nx, ny] = move(np.array([x1,y1]), potmove)
if(nx == x2 and ny == y2): count+=1
return count/9
transmat = np.zeros((64,64)) #reshape the original state matrix to (1,1), (1,2) ...
for S1 in range(64):
for(S2) in range(64):
transmat[S1, S2] = transprob(S1/8, S1%8, S2/8, S2%8)
plt.figure(figsize(8,7))
plt.matshow(transmat)
plt.xlabel('State at time t')
plt.ylabel('State at time t+1')
plt.title('Probabilities of transitioning from \none state to the next')
plt.colorbar();
<matplotlib.figure.Figure at 0x115330710>
emissionprob = np.zeros((64,8))
for State in range(64):
for Emission in range(8):
if State/8 == Emission: emissionprob[State, Emission] = 1
plt.matshow(emissionprob.T, origin='lower')
plt.colorbar(orientation='horizontal');
plt.ylabel('Observation')
plt.xlabel('State')
<matplotlib.text.Text at 0x11526fa90>
startprob = np.ones(64)/64
model = hmm.MultinomialHMM(n_components=64, transmat=transmat, startprob=startprob)
model.n_symbols_ = 8
model.emissionprob_ = emissionprob
#model.fit([positionlist['x'].tolist()])
predictions = model.predict(positionlist['x'].tolist())
proba = model.predict_proba(positionlist['x'].tolist())
step=15
plt.matshow(np.reshape(proba[step,:], (8,8), order='F'))
plt.colorbar();
y_estimate = pd.DataFrame(data = 1./8, columns=range(0,8), index=positionlist.index);
for step in y_estimate.index[:-1]:
x = positionlist['x'].ix[step]
y_estimate.ix[step] = proba[step,x*8:(x+1)*8]
plt.figure(figsize(20,6))
plt.subplot(211)
plt.imshow(y_estimate.T, origin='lower',interpolation='nearest')
plt.plot(positionlist['y'], 'wo')
plt.plot(predictions%8, 'r+')
plt.xlim(-.5, 49.5)
plt.ylim(-.5, 7.5)
plt.ylabel('Y position of knight');
plt.colorbar();
plt.subplot(212)
plt.imshow(y_estimate.T, origin='lower',interpolation='nearest')
plt.plot(positionlist['y'], 'wo')
plt.plot(predictions%8, 'r+')
plt.xlim(49.5,99.5)
plt.ylim(-.5, 7.5)
plt.ylabel('Y position of knight');
plt.xlabel('Timestep');
plt.colorbar();
startprob = np.zeros(64)
startprob[positionlist['x'].ix[0]*8+positionlist['y'].ix[0]] = 1
model = hmm.MultinomialHMM(n_components=64, transmat=transmat, startprob=startprob)
model.n_symbols_ = 8
model.emissionprob_ = emissionprob
#model.fit([positionlist['x'].tolist()])
predictions = model.predict(positionlist['x'].tolist())
proba = model.predict_proba(positionlist['x'].tolist())
y_estimate = pd.DataFrame(data = 1./8, columns=range(0,8), index=positionlist.index);
for step in y_estimate.index[:-1]:
x = positionlist['x'].ix[step]
y_estimate.ix[step] = proba[step,x*8:(x+1)*8]
plt.figure(figsize(20,6))
plt.subplot(211)
plt.imshow(y_estimate.T, origin='lower',interpolation='nearest', label='likelihood')
plt.plot(positionlist['y'], 'wo', label='Predictions')
plt.plot(predictions%8, 'r+', label='True Value')
plt.xlim(-.5, 49.5)
plt.ylim(-.5, 7.5)
plt.ylabel('Y position of knight');
plt.colorbar();
plt.legend(loc=(.37,1.05), ncol=2, frameon=False)
plt.subplot(212)
plt.imshow(y_estimate.T, origin='lower',interpolation='nearest')
plt.plot(positionlist['y'], 'wo')
plt.plot(predictions%8, 'r+')
plt.xlim(49.5,99.5)
plt.ylim(-.5, 7.5)
plt.ylabel('Y position of knight');
plt.xlabel('Timestep');
plt.colorbar();
So it looks like we won't ever get perfect prediction, because of the symmetry of the problem, and because it is inherently random - the most likely path isn't likely to be the one we take! The next problem arises because once you get off track, you aren't likely to get back on. So, can we periodically sample the y value to see if we're off track and get back on? Is there an optimal strategy for sampling y such that we minimize the cost while still meeting some performance target?
truthmatrix = pd.DataFrame(data=0, index=y_estimate.index, columns=y_estimate.columns.values)
for timestep, y in zip(positionlist.index, positionlist['y']):
truthmatrix[y].ix[timestep] = 1
delta = .02
accuracy = pd.Series(data=0., index=np.arange(0,1.01, delta))
numpredictions = pd.Series(data=0., index=np.arange(0,1.01, delta))
for val in accuracy.index:
numpredictions[val] = ((y_estimate > val-delta/2) & (y_estimate <= val+delta/2)).sum().sum()
numhits = truthmatrix[((y_estimate > val-delta/2) & (y_estimate <= val+delta/2))].sum().sum()
accuracy[val] = 1.*numhits/numpredictions[val]
plt.figure(figsize=(6,4))
plt.plot(accuracy.index, accuracy.values, 'bo')
plt.plot([0,1],[0,1],'r')
plt.title('Prediction Calibration')
plt.xlabel('Confidence / Expected accuracy')
plt.ylabel('Measured Accuracy')
<matplotlib.text.Text at 0x117abcbd0>
plt.figure(figsize=(6,4))
plt.semilogy(numpredictions.index, numpredictions.values, 'bo');
plt.title('Prediction Usefulness')
plt.ylabel('Number of Predictions')
plt.xlabel('Confidence / Expected accuracy')
<matplotlib.text.Text at 0x1179149d0>