.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/08_mocap6/plot-04-demo_ar_gaussian_prior_hypers.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-04-demo_ar_gaussian_prior_hypers.py: ================================================================ Visualizing learned state sequences and transition probabilities ================================================================ Train a sticky HMM model with auto-regressive Gaussian likelihoods on small motion capture data. Discover how the likelihood hyperparameters might impact performance. .. GENERATED FROM PYTHON SOURCE LINES 13-28 .. 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 29-30 Load dataset from file .. GENERATED FROM PYTHON SOURCE LINES 31-40 .. code-block:: default dataset_path = os.path.join(bnpy.DATASET_PATH, 'mocap6', 'dataset.npz') dataset = bnpy.data.GroupXData.read_npz(dataset_path) dataset_biasfeature = bnpy.data.GroupXData( X=dataset.X, Xprev=np.hstack([dataset.Xprev, np.ones((dataset.X.shape[0], 1))]), doc_range=dataset.doc_range) .. GENERATED FROM PYTHON SOURCE LINES 41-43 Setup: Function to make a simple plot of the raw data ----------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 44-101 .. 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 102-104 Visualization of the first sequence (1 of 6) -------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 105-109 .. code-block:: default show_single_sequence(0) .. image-sg:: /examples/08_mocap6/images/sphx_glr_plot-04-demo_ar_gaussian_prior_hypers_001.png :alt: plot 04 demo ar gaussian prior hypers :srcset: /examples/08_mocap6/images/sphx_glr_plot-04-demo_ar_gaussian_prior_hypers_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 110-112 Setup training: hyperparameters ---------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 113-127 .. code-block:: default K = 5 # Number of clusters/states # Allocation model (Standard Bayesian HMM with sticky self-transitions) 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' nTask = 1 .. GENERATED FROM PYTHON SOURCE LINES 128-140 Train with EM an HMM with *AutoRegGauss* observation model ---------------------------------------------------------- Train single model using likelihood without any free parameter $ x_t ~ \mbox{Normal}( A_k x_t-1, \Sigma_k) $ Take the best of 10 random initializations (in terms of evidence lower bound). .. GENERATED FROM PYTHON SOURCE LINES 141-156 .. code-block:: default bias0_trained_model, bias0_info_dict = bnpy.run( dataset, 'FiniteHMM', 'AutoRegGauss', 'EM', output_path=( '/tmp/mocap6/bias0-K=%d-model=HMM+AutoRegGauss-ECovMat=1*eye/' % (K)), nLap=100, nTask=nTask, convergeThr=0.0001, transAlpha=transAlpha, startAlpha=startAlpha, hmmKappa=hmmKappa, sF=sF, ECovMat=ECovMat, MMat='eye', K=K, initname='randexamples', printEvery=25, ) .. rst-class:: sphx-glr-script-out .. code-block:: none Dataset Summary: GroupXData size: 6 units (documents) dimension: 12 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: EM | task 1/1 | alg. seed: 1607680 | data order seed: 8541952 task_output_path: /tmp/mocap6/bias0-K=5-model=HMM+AutoRegGauss-ECovMat=1*eye/1 1/100 after 0 sec. | 229.7 MiB | K 5 | loss 1.029847757e+10 | 2/100 after 0 sec. | 229.7 MiB | K 5 | loss 2.401709596e+00 | Ndiff 383.432 25/100 after 2 sec. | 229.7 MiB | K 5 | loss 2.054112455e+00 | Ndiff 0.371 50/100 after 4 sec. | 229.7 MiB | K 5 | loss 2.046961483e+00 | Ndiff 5.413 75/100 after 7 sec. | 229.7 MiB | K 5 | loss 2.021134069e+00 | Ndiff 0.864 100/100 after 9 sec. | 229.7 MiB | K 5 | loss 2.012157346e+00 | Ndiff 1.147 ... done. not converged. max laps thru data exceeded. .. GENERATED FROM PYTHON SOURCE LINES 157-171 Train with EM an HMM with *AutoRegGauss* observation model ---------------------------------------------------------- Train single model using likelihood WITH free mean parameter $ x_t ~ \mbox{Normal}( A_k x_t-1 + \mu_k, \Sigma_k) $ This is equivalent to including column of all ones into the x-previous, which is how we do it in practice... Take the best of 10 random initializations (in terms of evidence lower bound). .. GENERATED FROM PYTHON SOURCE LINES 172-187 .. code-block:: default bias1_trained_model, bias1_info_dict = bnpy.run( dataset_biasfeature, 'FiniteHMM', 'AutoRegGauss', 'EM', output_path=( '/tmp/mocap6/bias1-K=%d-model=HMM+AutoRegGauss-ECovMat=1*eye/' % (K)), nLap=100, nTask=nTask, convergeThr=0.0001, transAlpha=transAlpha, startAlpha=startAlpha, hmmKappa=hmmKappa, sF=sF, ECovMat=ECovMat, MMat='eye', K=K, initname='randexamples', printEvery=25, ) .. rst-class:: sphx-glr-script-out .. code-block:: none Dataset Summary: GroupXData size: 6 units (documents) dimension: 12 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: EM | task 1/1 | alg. seed: 1607680 | data order seed: 8541952 task_output_path: /tmp/mocap6/bias1-K=5-model=HMM+AutoRegGauss-ECovMat=1*eye/1 1/100 after 0 sec. | 229.7 MiB | K 5 | loss 1.029794586e+10 | 2/100 after 0 sec. | 229.7 MiB | K 5 | loss 2.400629439e+00 | Ndiff 352.956 25/100 after 2 sec. | 229.7 MiB | K 5 | loss 2.022326047e+00 | Ndiff 3.077 50/100 after 4 sec. | 229.7 MiB | K 5 | loss 2.004862867e+00 | Ndiff 0.832 75/100 after 7 sec. | 229.7 MiB | K 5 | loss 1.987434495e+00 | Ndiff 0.006 90/100 after 8 sec. | 229.7 MiB | K 5 | loss 1.987434487e+00 | Ndiff 0.000 ... done. converged. .. GENERATED FROM PYTHON SOURCE LINES 188-201 Train with VB an HMM with *AutoRegGauss* observation model ---------------------------------------------------------- Train single model using likelihood with free mean parameter $ x_t ~ \mbox{Normal}( A_k x_t-1 + \mu_k, \Sigma_k) $ This is equivalent to including column of all ones into the x-previous. Take the best of 10 random initializations (in terms of evidence lower bound). .. GENERATED FROM PYTHON SOURCE LINES 202-231 .. code-block:: default bias1vb_trained_model, bias1vb_info_dict = bnpy.run( dataset_biasfeature, 'FiniteHMM', 'AutoRegGauss', 'VB', output_path=( '/tmp/mocap6/bias1vb-K=%d-model=HMM+AutoRegGauss-ECovMat=1*eye/' % (K)), nLap=100, nTask=nTask, convergeThr=0.0001, transAlpha=transAlpha, startAlpha=startAlpha, hmmKappa=hmmKappa, sF=sF, ECovMat=ECovMat, MMat='eye', K=K, initname=bias1_info_dict['task_output_path'], printEvery=25, ) bias0vb_trained_model, bias0vb_info_dict = bnpy.run( dataset, 'FiniteHMM', 'AutoRegGauss', 'VB', output_path=( '/tmp/mocap6/bias1vb-K=%d-model=HMM+AutoRegGauss-ECovMat=1*eye/' % (K)), nLap=100, nTask=nTask, convergeThr=0.0001, transAlpha=transAlpha, startAlpha=startAlpha, hmmKappa=hmmKappa, sF=sF, ECovMat=ECovMat, MMat='eye', K=K, initname=bias0_info_dict['task_output_path'], printEvery=25, ) .. rst-class:: sphx-glr-script-out .. code-block:: none Dataset Summary: GroupXData size: 6 units (documents) dimension: 12 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 = /tmp/mocap6/bias1-K=5-model=HMM+AutoRegGauss-ECovMat=1*eye/1 K = 5 (number of clusters) seed = 1607680 elapsed_time: 0.1 sec Learn Alg: VB | task 1/1 | alg. seed: 1607680 | data order seed: 8541952 task_output_path: /tmp/mocap6/bias1vb-K=5-model=HMM+AutoRegGauss-ECovMat=1*eye/1 1/100 after 0 sec. | 229.7 MiB | K 5 | loss 2.279007765e+00 | 2/100 after 0 sec. | 229.7 MiB | K 5 | loss 2.278850534e+00 | Ndiff 2.120 25/100 after 3 sec. | 229.7 MiB | K 5 | loss 2.277292722e+00 | Ndiff 1.431 50/100 after 6 sec. | 229.7 MiB | K 5 | loss 2.275682910e+00 | Ndiff 0.113 72/100 after 8 sec. | 229.7 MiB | K 5 | loss 2.275561264e+00 | Ndiff 0.000 ... done. converged. Dataset Summary: GroupXData size: 6 units (documents) dimension: 12 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 = /tmp/mocap6/bias0-K=5-model=HMM+AutoRegGauss-ECovMat=1*eye/1 K = 5 (number of clusters) seed = 1607680 elapsed_time: 0.1 sec Learn Alg: VB | task 1/1 | alg. seed: 1607680 | data order seed: 8541952 task_output_path: /tmp/mocap6/bias1vb-K=5-model=HMM+AutoRegGauss-ECovMat=1*eye/1 1/100 after 0 sec. | 229.7 MiB | K 5 | loss 2.274308619e+00 | 2/100 after 0 sec. | 229.7 MiB | K 5 | loss 2.274112727e+00 | Ndiff 1.193 25/100 after 3 sec. | 229.7 MiB | K 5 | loss 2.273779830e+00 | Ndiff 0.318 50/100 after 5 sec. | 229.7 MiB | K 5 | loss 2.273609130e+00 | Ndiff 0.003 75/100 after 8 sec. | 229.7 MiB | K 5 | loss 2.273609126e+00 | Ndiff 0.001 100/100 after 11 sec. | 229.7 MiB | K 5 | loss 2.273609125e+00 | Ndiff 0.000 ... done. not converged. max laps thru data exceeded. .. GENERATED FROM PYTHON SOURCE LINES 232-235 Compare loss function traces for all methods -------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 236-269 .. code-block:: default pylab.figure() markersize = 5 pylab.plot( bias0_info_dict['lap_history'], bias0_info_dict['loss_history'], 'bs-', markersize=markersize, label='EM without bias feature') pylab.plot( bias1_info_dict['lap_history'], bias1_info_dict['loss_history'], 'ks-', markersize=markersize, label='EM WITH bias feature') pylab.plot( bias0vb_info_dict['lap_history'], bias0vb_info_dict['loss_history'], 'bd--', markersize=markersize, label='VB without bias feature') pylab.plot( bias1vb_info_dict['lap_history'], bias1vb_info_dict['loss_history'], 'kd--', markersize=markersize, label='VB WITH bias feature') pylab.legend(loc='upper right') pylab.ylim([1.5, 3.0]) pylab.xlabel('num. laps') pylab.ylabel('loss: - log p(x)') pylab.draw() pylab.tight_layout() pylab.show(block=False) .. image-sg:: /examples/08_mocap6/images/sphx_glr_plot-04-demo_ar_gaussian_prior_hypers_002.png :alt: plot 04 demo ar gaussian prior hypers :srcset: /examples/08_mocap6/images/sphx_glr_plot-04-demo_ar_gaussian_prior_hypers_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 36.140 seconds) .. _sphx_glr_download_examples_08_mocap6_plot-04-demo_ar_gaussian_prior_hypers.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-04-demo_ar_gaussian_prior_hypers.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot-04-demo_ar_gaussian_prior_hypers.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_