Visualizing learned state sequences and transition probabilities¶

Train a sticky HDP-HMM model on small motion capture data, then visualize the MAP state sequences under the estimated model parameters by running Viterbi.

Also has some info on how to inspect the learned HMM parameters of a sticky HDP-HMM model trained on small motion capture data.

```# sphinx_gallery_thumbnail_number = 3

import bnpy
import numpy as np
import os

import matplotlib
from matplotlib import pylab
import seaborn as sns

np.set_printoptions(suppress=1, precision=3)

FIG_SIZE = (10, 5)
pylab.rcParams['figure.figsize'] = FIG_SIZE
```

```dataset_path = os.path.join(bnpy.DATASET_PATH, 'mocap6')
os.path.join(dataset_path, 'dataset.npz'))
```

Setup: Function to make a simple plot of the raw data¶

```def show_single_sequence(
seq_id,
zhat_T=None,
z_img_cmap=None,
ylim=[-120, 120],
K=5,
left=0.2, bottom=0.2, right=0.8, top=0.95):
if z_img_cmap is None:
z_img_cmap = matplotlib.cm.get_cmap('Set1', K)

if zhat_T is None:
nrows = 1
else:
nrows = 2
fig_h, ax_handles = pylab.subplots(
nrows=nrows, ncols=1, sharex=True, sharey=False)
ax_handles = np.atleast_1d(ax_handles).flatten().tolist()

start = dataset.doc_range[seq_id]
stop = dataset.doc_range[seq_id + 1]
# Extract current sequence
# as a 2D array : T x D (n_timesteps x n_dims)
curX_TD = dataset.X[start:stop]
for dim in range(12):
ax_handles[0].plot(curX_TD[:, dim], '.-')
ax_handles[0].set_ylabel('angle')
ax_handles[0].set_ylim(ylim)
z_img_height = int(np.ceil(ylim[1] - ylim[0]))
wspace=0.1,
hspace=0.1,
left=left, right=right,
bottom=bottom, top=top)
if zhat_T is not None:
img_TD = np.tile(zhat_T, (z_img_height, 1))
ax_handles[1].imshow(
img_TD,
interpolation='nearest',
vmin=-0.5, vmax=(K-1)+0.5,
cmap=z_img_cmap)
ax_handles[1].set_ylim(0, z_img_height)
ax_handles[1].set_yticks([])

bbox = ax_handles[1].get_position()
width = (1.0 - bbox.x1) / 3
height = bbox.y1 - bbox.y0
cax = fig_h.add_axes([right + 0.01, bottom, width, height])
cbax_h = fig_h.colorbar(
ax_handles[1].images[0], cax=cax, orientation='vertical')
cbax_h.set_ticks(np.arange(K))
cbax_h.set_ticklabels(np.arange(K))
cbax_h.ax.tick_params(labelsize=9)

ax_handles[-1].set_xlabel('time')
return ax_handles
```

Visualization of the first sequence (1 of 6)¶

```show_single_sequence(0)
```
```[<AxesSubplot: xlabel='time', ylabel='angle'>]
```

Setup: hyperparameters¶

```K = 5            # Number of clusters/states

# Allocation model (HDP)
gamma      =   5.0  # top-level Dirichlet concentration parameter
transAlpha =   0.5  # trans-level Dirichlet concentration parameter
startAlpha =  10.0  # starting-state Dirichlet concentration parameter
hmmKappa   =  50.0  # set sticky self-transition weight

# Observation model (1st-order Auto-regressive Gaussian)
sF = 1.0          # Set observation model prior so E[covariance] = identity
ECovMat = 'eye'
```

Train HDP-HMM with AutoRegGauss observation model¶

Train single model for all 6 sequences.

Do small number of clusters jut to make visualization easy.

Take the best of 5 random initializations (in terms of evidence lower bound).

```hdphmm_trained_model, hmmar_info_dict = bnpy.run(
dataset, 'HDPHMM', 'AutoRegGauss', 'memoVB',
output_path=(
'/tmp/mocap6/showcase-K=%d-model=HDPHMM+AutoRegGauss-ECovMat=1*eye/'
% (K)),
transAlpha=transAlpha, startAlpha=startAlpha, hmmKappa=hmmKappa,
gamma=gamma,
sF=sF, ECovMat=ECovMat,
K=K, initname='randexamples',
printEvery=25,
)
```
```Dataset Summary:
GroupXData
total size: 6 units
batch size: 6 units
num. batches: 1
Allocation Model:  None
Obs. Data  Model:  Auto-Regressive Gaussian with full covariance.
Obs. Data  Prior:  MatrixNormal-Wishart on each mean/prec matrix pair: A, Lam
E[ A ] =
[[1. 0.]
[0. 1.]] ...
E[ Sigma ] =
[[1. 0.]
[0. 1.]] ...
Initialization:
initname = randexamples
K = 5 (number of clusters)
seed = 1607680
elapsed_time: 0.0 sec
Learn Alg: memoVB | task  1/5 | alg. seed: 1607680 | data order seed: 8541952
1.000/100 after      0 sec. |    229.7 MiB | K    5 | loss  2.772882295e+00 |
2.000/100 after      0 sec. |    229.7 MiB | K    5 | loss  2.502603548e+00 | Ndiff  470.393
3.000/100 after      0 sec. |    229.7 MiB | K    5 | loss  2.426864576e+00 | Ndiff  119.274
25.000/100 after      3 sec. |    229.7 MiB | K    5 | loss  2.302556401e+00 | Ndiff   18.913
50.000/100 after      6 sec. |    229.7 MiB | K    5 | loss  2.299789291e+00 | Ndiff    1.326
75.000/100 after      9 sec. |    229.7 MiB | K    5 | loss  2.285556973e+00 | Ndiff    1.744
100.000/100 after     12 sec. |    229.7 MiB | K    5 | loss  2.268502057e+00 | Ndiff    1.794
... done. not converged. max laps thru data exceeded.
Learn Alg: memoVB | task  2/5 | alg. seed: 6454144 | data order seed: 7673856
1.000/100 after      0 sec. |    229.7 MiB | K    5 | loss  2.789290311e+00 |
2.000/100 after      0 sec. |    229.7 MiB | K    5 | loss  2.547246672e+00 | Ndiff  332.182
3.000/100 after      0 sec. |    229.7 MiB | K    5 | loss  2.446687466e+00 | Ndiff   96.506
25.000/100 after      3 sec. |    229.7 MiB | K    5 | loss  2.282761453e+00 | Ndiff    2.677
50.000/100 after      6 sec. |    229.7 MiB | K    5 | loss  2.279866842e+00 | Ndiff    0.028
75.000/100 after      9 sec. |    229.7 MiB | K    5 | loss  2.279866670e+00 | Ndiff    0.003
100.000/100 after     12 sec. |    229.7 MiB | K    5 | loss  2.279866668e+00 | Ndiff    0.000
... done. not converged. max laps thru data exceeded.
Learn Alg: memoVB | task  3/5 | alg. seed: 6168832 | data order seed: 7360256
1.000/100 after      0 sec. |    229.7 MiB | K    5 | loss  2.838795721e+00 |
2.000/100 after      0 sec. |    229.7 MiB | K    5 | loss  2.614908627e+00 | Ndiff  188.922
3.000/100 after      0 sec. |    229.7 MiB | K    5 | loss  2.526132362e+00 | Ndiff  143.993
25.000/100 after      3 sec. |    229.7 MiB | K    5 | loss  2.319895680e+00 | Ndiff    1.371
50.000/100 after      6 sec. |    229.7 MiB | K    5 | loss  2.318585849e+00 | Ndiff    0.184
75.000/100 after      9 sec. |    229.7 MiB | K    5 | loss  2.316143371e+00 | Ndiff    0.030
100.000/100 after     12 sec. |    229.7 MiB | K    5 | loss  2.316143320e+00 | Ndiff    0.001
... done. not converged. max laps thru data exceeded.
Learn Alg: memoVB | task  4/5 | alg. seed: 8205184 | data order seed: 900864
1.000/100 after      0 sec. |    229.7 MiB | K    5 | loss  2.805840898e+00 |
2.000/100 after      0 sec. |    229.7 MiB | K    5 | loss  2.587110039e+00 | Ndiff  428.347
3.000/100 after      0 sec. |    229.7 MiB | K    5 | loss  2.495055939e+00 | Ndiff  176.781
25.000/100 after      3 sec. |    229.7 MiB | K    5 | loss  2.338830039e+00 | Ndiff    0.020
47.000/100 after      6 sec. |    229.7 MiB | K    5 | loss  2.338829969e+00 | Ndiff    0.000
... done. converged.
Learn Alg: memoVB | task  5/5 | alg. seed: 2986368 | data order seed: 6479872
1.000/100 after      0 sec. |    229.7 MiB | K    5 | loss  2.828188354e+00 |
2.000/100 after      0 sec. |    229.7 MiB | K    5 | loss  2.592647224e+00 | Ndiff  435.868
3.000/100 after      0 sec. |    229.7 MiB | K    5 | loss  2.478529769e+00 | Ndiff  171.138
25.000/100 after      3 sec. |    229.7 MiB | K    5 | loss  2.302103693e+00 | Ndiff    0.362
50.000/100 after      6 sec. |    229.7 MiB | K    5 | loss  2.301974529e+00 | Ndiff    0.194
75.000/100 after      9 sec. |    229.7 MiB | K    5 | loss  2.301969041e+00 | Ndiff    0.000
... done. converged.
```

Visualize the starting-state probabilities¶

start_prob_K1D array, size K

start_prob_K[k] = exp( E[log Pr(start state = k)] )

```start_prob_K = hdphmm_trained_model.allocModel.get_init_prob_vector()

print(start_prob_K)
```
```[0.249 0.336 0.032 0.025 0.168]
```

Visualize the transition probabilities¶

trans_prob_KK2D array, K x K

trans_prob_KK[j, k] = exp( E[log Pr(z_t = k | z_t-1 = j)] )

```trans_prob_KK = hdphmm_trained_model.allocModel.get_trans_prob_matrix()

print(trans_prob_KK)
```
```[[0.917 0.008 0.052 0.009 0.011]
[0.025 0.953 0.    0.    0.019]
[0.09  0.011 0.861 0.035 0.   ]
[0.012 0.    0.077 0.907 0.   ]
[0.021 0.026 0.    0.    0.95 ]]
```

Compute log likelihood of each timestep for sequence 0¶

log_lik_TK2D array, T x K

log_lik_TK[t, k] = E[ log Pr( observed data at time t | z_t = k)]

```log_lik_seq0_TK = hdphmm_trained_model.obsModel.calcLogSoftEvMatrix_FromPost(
dataset.make_subset([0])
)

print(log_lik_seq0_TK[:10, :])
```
```[[ -135.664   -23.757   -42.577  -122.033 -1174.337]
[ -505.011   -28.26    -52.154  -483.826 -4680.957]
[  -82.065   -24.821   -45.635   -58.799  -179.048]
[ -168.085   -28.907   -37.378  -149.684  -678.56 ]
[ -302.      -28.805   -35.963  -282.073 -2687.59 ]
[  -70.994   -25.882   -34.915   -68.371  -535.383]
[ -437.198   -43.6     -50.451  -369.265 -3470.826]
[  -99.618   -26.81    -51.818   -72.34   -142.138]
[ -162.653   -32.128   -48.969  -163.484  -942.637]
[ -582.385   -28.568   -77.439  -573.153 -5090.44 ]]
```

Run Viterbi algorithm for sequence 0¶

zhat_T1D array, size T

MAP state sequence zhat_T[t] = state assigned to timestep t, will be int value in {0, 1, … K-1}

```zhat_seq0_T = bnpy.allocmodel.hmm.HMMUtil.runViterbiAlg(
log_lik_seq0_TK, np.log(start_prob_K), np.log(trans_prob_KK))

print(zhat_seq0_T[:10])
```
```[1 1 1 1 1 1 1 1 1 1]
```

Visualize the segmentation for sequence 0¶

```show_single_sequence(0, zhat_T=zhat_seq0_T, K=K)
```
```[<AxesSubplot: ylabel='angle'>, <AxesSubplot: xlabel='time'>]
```

Visualize the segmentation for sequence 1¶

```log_lik_seq1_TK = hdphmm_trained_model.obsModel.calcLogSoftEvMatrix_FromPost(
dataset.make_subset([1])
)
zhat_seq1_T = bnpy.allocmodel.hmm.HMMUtil.runViterbiAlg(
log_lik_seq1_TK, np.log(start_prob_K), np.log(trans_prob_KK))

show_single_sequence(1, zhat_T=zhat_seq1_T, K=K)

pylab.show(block=False)
```

Total running time of the script: ( 0 minutes 51.747 seconds)

Gallery generated by Sphinx-Gallery