.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/08_mocap6/plot-03-demo=interpret_hdphmm_params_and_run_viterbi.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note Click :ref:`here ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_08_mocap6_plot-03-demo=interpret_hdphmm_params_and_run_viterbi.py: ================================================================ 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. .. GENERATED FROM PYTHON SOURCE LINES 11-26 .. code-block:: default # 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 .. GENERATED FROM PYTHON SOURCE LINES 27-28 Load dataset from file .. GENERATED FROM PYTHON SOURCE LINES 29-34 .. code-block:: default dataset_path = os.path.join(bnpy.DATASET_PATH, 'mocap6') dataset = bnpy.data.GroupXData.read_npz( os.path.join(dataset_path, 'dataset.npz')) .. GENERATED FROM PYTHON SOURCE LINES 35-37 Setup: Function to make a simple plot of the raw data ----------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 38-95 .. code-block:: default 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])) pylab.subplots_adjust( 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 .. GENERATED FROM PYTHON SOURCE LINES 96-98 Visualization of the first sequence (1 of 6) -------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 99-103 .. code-block:: default show_single_sequence(0) .. image-sg:: /examples/08_mocap6/images/sphx_glr_plot-03-demo=interpret_hdphmm_params_and_run_viterbi_001.png :alt: plot 03 demo=interpret hdphmm params and run viterbi :srcset: /examples/08_mocap6/images/sphx_glr_plot-03-demo=interpret_hdphmm_params_and_run_viterbi_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 104-106 Setup: hyperparameters ---------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 107-120 .. code-block:: default 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' .. GENERATED FROM PYTHON SOURCE LINES 121-130 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). .. GENERATED FROM PYTHON SOURCE LINES 131-146 .. code-block:: default 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)), nLap=100, nTask=5, nBatch=1, convergeThr=0.0001, transAlpha=transAlpha, startAlpha=startAlpha, hmmKappa=hmmKappa, gamma=gamma, sF=sF, ECovMat=ECovMat, K=K, initname='randexamples', printEvery=25, ) .. rst-class:: sphx-glr-script-out .. code-block:: none 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 task_output_path: /tmp/mocap6/showcase-K=5-model=HDPHMM+AutoRegGauss-ECovMat=1*eye/1 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 task_output_path: /tmp/mocap6/showcase-K=5-model=HDPHMM+AutoRegGauss-ECovMat=1*eye/2 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 task_output_path: /tmp/mocap6/showcase-K=5-model=HDPHMM+AutoRegGauss-ECovMat=1*eye/3 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 task_output_path: /tmp/mocap6/showcase-K=5-model=HDPHMM+AutoRegGauss-ECovMat=1*eye/4 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 task_output_path: /tmp/mocap6/showcase-K=5-model=HDPHMM+AutoRegGauss-ECovMat=1*eye/5 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. .. GENERATED FROM PYTHON SOURCE LINES 147-152 Visualize the starting-state probabilities ------------------------------------------ start_prob_K : 1D array, size K start_prob_K[k] = exp( E[log Pr(start state = k)] ) .. GENERATED FROM PYTHON SOURCE LINES 153-158 .. code-block:: default start_prob_K = hdphmm_trained_model.allocModel.get_init_prob_vector() print(start_prob_K) .. rst-class:: sphx-glr-script-out .. code-block:: none [0.249 0.336 0.032 0.025 0.168] .. GENERATED FROM PYTHON SOURCE LINES 159-165 Visualize the transition probabilities -------------------------------------- trans_prob_KK : 2D array, K x K trans_prob_KK[j, k] = exp( E[log Pr(z_t = k | z_t-1 = j)] ) .. GENERATED FROM PYTHON SOURCE LINES 166-172 .. code-block:: default trans_prob_KK = hdphmm_trained_model.allocModel.get_trans_prob_matrix() print(trans_prob_KK) .. rst-class:: sphx-glr-script-out .. code-block:: none [[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 ]] .. GENERATED FROM PYTHON SOURCE LINES 173-178 Compute log likelihood of each timestep for sequence 0 ------------------------------------------------------ log_lik_TK : 2D array, T x K log_lik_TK[t, k] = E[ log Pr( observed data at time t | z_t = k)] .. GENERATED FROM PYTHON SOURCE LINES 179-187 .. code-block:: default log_lik_seq0_TK = hdphmm_trained_model.obsModel.calcLogSoftEvMatrix_FromPost( dataset.make_subset([0]) ) print(log_lik_seq0_TK[:10, :]) .. rst-class:: sphx-glr-script-out .. code-block:: none [[ -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 ]] .. GENERATED FROM PYTHON SOURCE LINES 188-194 Run Viterbi algorithm for sequence 0 ------------------------------------ zhat_T : 1D array, size T MAP state sequence zhat_T[t] = state assigned to timestep t, will be int value in {0, 1, ... K-1} .. GENERATED FROM PYTHON SOURCE LINES 195-202 .. code-block:: default 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]) .. rst-class:: sphx-glr-script-out .. code-block:: none [1 1 1 1 1 1 1 1 1 1] .. GENERATED FROM PYTHON SOURCE LINES 203-206 Visualize the segmentation for sequence 0 ----------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 207-211 .. code-block:: default show_single_sequence(0, zhat_T=zhat_seq0_T, K=K) .. image-sg:: /examples/08_mocap6/images/sphx_glr_plot-03-demo=interpret_hdphmm_params_and_run_viterbi_002.png :alt: plot 03 demo=interpret hdphmm params and run viterbi :srcset: /examples/08_mocap6/images/sphx_glr_plot-03-demo=interpret_hdphmm_params_and_run_viterbi_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [, ] .. GENERATED FROM PYTHON SOURCE LINES 212-215 Visualize the segmentation for sequence 1 ----------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 216-227 .. code-block:: default 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) .. image-sg:: /examples/08_mocap6/images/sphx_glr_plot-03-demo=interpret_hdphmm_params_and_run_viterbi_003.png :alt: plot 03 demo=interpret hdphmm params and run viterbi :srcset: /examples/08_mocap6/images/sphx_glr_plot-03-demo=interpret_hdphmm_params_and_run_viterbi_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 51.747 seconds) .. _sphx_glr_download_examples_08_mocap6_plot-03-demo=interpret_hdphmm_params_and_run_viterbi.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot-03-demo=interpret_hdphmm_params_and_run_viterbi.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot-03-demo=interpret_hdphmm_params_and_run_viterbi.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_