pyemma.msm.estimate_markov_model

pyemma.msm.estimate_markov_model(dtrajs, lag, reversible=True, statdist=None, count_mode='sliding', weights='empirical', sparse=False, connectivity='largest', dt_traj='1 step', maxiter=1000000, maxerr=1e-08, score_method='VAMP2', score_k=10, mincount_connectivity='1/n', core_set=None, milestoning_method='last_core')

Estimates a Markov model from discrete trajectories

Returns a MaximumLikelihoodMSM that contains the estimated transition matrix and allows to compute a large number of quantities related to Markov models.

Parameters

  • dtrajs (list containing ndarrays(dtype=int) or ndarray(n, dtype=int)) – discrete trajectories, stored as integer ndarrays (arbitrary size) or a single ndarray for only one trajectory.

  • lag (int) – lag time at which transitions are counted and the transition matrix is estimated.

  • reversible (bool, optional) – If true compute reversible MSM, else non-reversible MSM

  • statdist ((M,) ndarray, optional) – Stationary vector on the full state-space. Transition matrix will be estimated such that statdist is its equilibrium distribution.

  • count_mode (str, optional, default='sliding') –

    mode to obtain count matrices from discrete trajectories. Should be one of:

    • ’sliding’ : A trajectory of length T will have \(T-\tau\) counts at time indexes

      \[(0 \rightarrow \tau), (1 \rightarrow \tau+1), ..., (T-\tau-1 \rightarrow T-1)\]

    • ’effective’ : Uses an estimate of the transition counts that are statistically uncorrelated. Recommended when used with a Bayesian MSM.

    • ’sample’ : A trajectory of length T will have \(T/\tau\) counts at time indexes

      \[(0 \rightarrow \tau), (\tau \rightarrow 2 \tau), ..., (((T/\tau)-1) \tau \rightarrow T)\]

  • weights (str, optional) –

    can be used to re-weight non-equilibrium data to equilibrium. Must be one of the following:

    • ’empirical’: Each trajectory frame counts as one. (default)

    • ’oom’: Each transition is re-weighted using OOM theory, see 11.

  • sparse (bool, optional) – If true compute count matrix, transition matrix and all derived quantities using sparse matrix algebra. In this case python sparse matrices will be returned by the corresponding functions instead of numpy arrays. This behavior is suggested for very large numbers of states (e.g. > 4000) because it is likely to be much more efficient.

  • connectivity (str, optional) –

    Connectivity mode. Three methods are intended (currently only ‘largest’ is implemented)

    • ’largest’ : The active set is the largest reversibly connected set. All estimation will be done on this subset and all quantities (transition matrix, stationary distribution, etc) are only defined on this subset and are correspondingly smaller than the full set of states

    • ’all’ : The active set is the full set of states. Estimation will be conducted on each reversibly connected set separately. That means the transition matrix will decompose into disconnected submatrices, the stationary vector is only defined within subsets, etc. Currently not implemented.

    • ’none’ : The active set is the full set of states. Estimation will be conducted on the full set of states without ensuring connectivity. This only permits nonreversible estimation. Currently not implemented.

  • dt_traj (str, optional) –

    Description of the physical time corresponding to the lag. May be used by analysis algorithms such as plotting tools to pretty-print the axes. By default ‘1 step’, i.e. there is no physical time unit. Specify by a number, whitespace and unit. Permitted units are (* is an arbitrary string):

    • ’fs’, ‘femtosecond*’

    • ’ps’, ‘picosecond*’

    • ’ns’, ‘nanosecond*’

    • ’us’, ‘microsecond*’

    • ’ms’, ‘millisecond*’

    • ’s’, ‘second*’

  • maxiter (int, optional) – Optional parameter with reversible = True. maximum number of iterations before the transition matrix estimation method exits

  • maxerr (float, optional) – Optional parameter with reversible = True. convergence tolerance for transition matrix estimation. This specifies the maximum change of the Euclidean norm of relative stationary probabilities (\(x_i = \sum_k x_{ik}\)). The relative stationary probability changes \(e_i = (x_i^{(1)} - x_i^{(2)})/(x_i^{(1)} + x_i^{(2)})\) are used in order to track changes in small probabilities. The Euclidean norm of the change vector, \(|e_i|_2\), is compared to maxerr.

  • score_method (str, optional, default='VAMP2') –

    Score to be used with MSM score function. Available scores are based on the variational approach for Markov processes 13 14:

    • ’VAMP1’ Sum of singular values of the symmetrized transition matrix 14 .

      If the MSM is reversible, this is equal to the sum of transition matrix eigenvalues, also called Rayleigh quotient 13 15 .

    • ’VAMP2’ Sum of squared singular values of the symmetrized transition matrix 14 .

      If the MSM is reversible, this is equal to the kinetic variance 16 .

  • score_k (int or None) – The maximum number of eigenvalues or singular values used in the score. If set to None, all available eigenvalues will be used.

  • mincount_connectivity (float or '1/n') – minimum number of counts to consider a connection between two states. Counts lower than that will count zero in the connectivity check and may thus separate the resulting transition matrix. The default evaluates to 1/nstates.

  • core_set (None (default) or array like, dtype=int) – Definition of core set for milestoning MSMs. If set to None, replaces state -1 (if found in discrete trajectories) and performs milestone counting. No effect for Voronoi-discretized trajectories (default). If a list or np.ndarray is supplied, discrete trajectories will be assigned accordingly.

  • milestoning_method (str) – Method to use for counting transitions in trajectories with unassigned frames. Currently available: | ‘last_core’, assigns unassigned frames to last visited core

Returns

msm – Estimator object containing the MSM and estimation information.

Return type

MaximumLikelihoodMSM

See also

MaximumLikelihoodMSM

An MSM object that has been estimated from data

class pyemma.msm.estimators.maximum_likelihood_msm.MaximumLikelihoodMSM(*args, **kwargs)

Maximum likelihood estimator for MSMs given discrete trajectory statistics

Methods

cktest(nsets[, memberships, mlags, conf, …])

Conducts a Chapman-Kolmogorow test.

coarse_grain(ncoarse[, method])

Returns a coarse-grained Markov model.

committor_backward(A, B)

Backward committor from set A to set B

committor_forward(A, B)

Forward committor (also known as p_fold or splitting probability) from set A to set B

correlation(a[, b, maxtime, k, ncv])

Time-correlation for equilibrium experiment.

eigenvalues([k])

Compute the transition matrix eigenvalues

eigenvectors_left([k])

Compute the left transition matrix eigenvectors

eigenvectors_right([k])

Compute the right transition matrix eigenvectors

estimate(dtrajs, **kwargs)

param dtrajs

discrete trajectories, stored as integer ndarrays (arbitrary size)

expectation(a)

Equilibrium expectation value of a given observable.

fingerprint_correlation(a[, b, k, ncv])

Dynamical fingerprint for equilibrium time-correlation experiment.

fingerprint_relaxation(p0, a[, k, ncv])

Dynamical fingerprint for perturbation/relaxation experiment.

fit(X[, y])

Estimates parameters - for compatibility with sklearn.

generate_traj(N[, start, stop, stride])

Generates a synthetic discrete trajectory of length N and simulation time stride * lag time * N

get_model_params([deep])

Get parameters for this model.

get_params([deep])

Get parameters for this estimator.

hmm(nhidden)

Estimates a hidden Markov state model as described in 1

load(file_name[, model_name])

Loads a previously saved PyEMMA object from disk.

mfpt(A, B)

Mean first passage times from set A to set B, in units of the input trajectory time step

pcca(m)

Runs PCCA++ 1 to compute a metastable decomposition of MSM states

propagate(p0, k)

Propagates the initial distribution p0 k times

relaxation(p0, a[, maxtime, k, ncv])

Simulates a perturbation-relaxation experiment.

sample_by_distributions(distributions, nsample)

Generates samples according to given probability distributions

sample_by_state(nsample[, subset, replace])

Generates samples of the connected states.

save(file_name[, model_name, overwrite, …])

saves the current state of this object to given file and name.

score(dtrajs[, score_method, score_k])

Scores the MSM using the dtrajs using the variational approach for Markov processes 1 2

score_cv(dtrajs[, n, score_method, score_k])

Scores the MSM using the variational approach for Markov processes 1 2 and crossvalidation 3 .

set_model_params(P[, pi, reversible, …])

Call to set all basic model parameters.

set_params(**params)

Set the parameters of this estimator.

simulate(N[, start, stop, dt])

Generates a realization of the Markov Model

timescales([k])

The relaxation timescales corresponding to the eigenvalues

trajectory_weights()

Uses the MSM to assign a probability weight to each trajectory frame.

update_model_params(**params)

Update given model parameter if they are set to specific values

Attributes

property effective_count_matrix

Statistically uncorrelated transition counts within the active set of states

You can use this count matrix for Bayesian estimation or error perturbation.

References

[1] Noe, F. (2015) Statistical inefficiency of Markov model count matrices

http://publications.mi.fu-berlin.de/1699/1/autocorrelation_counts.pdf

References

The mathematical theory of Markov (state) model estimation was introduced in 1 . Further theoretical developments were made in 2 . The term Markov state model was coined in 3 . Continuous-time Markov models (Master equation models) were suggested in 4. Reversible Markov model estimation was introduced in 5 , and further developed in 6 7 9 . It was shown in 8 that the quality of Markov state models does in fact not depend on memory loss, but rather on where the discretization is suitable to approximate the eigenfunctions of the Markov operator (the ‘reaction coordinates’). With a suitable choice of discretization and lag time, MSMs can thus become very accurate. 9 introduced a number of methodological improvements and gives a good overview of the methodological basics of Markov state modeling today. 10 is a more extensive review book of theory, methods and applications.

1(1,2,3,4,5)

Schuette, C. , A. Fischer, W. Huisinga and P. Deuflhard: A Direct Approach to Conformational Dynamics based on Hybrid Monte Carlo. J. Comput. Phys., 151, 146-168 (1999)

2(1,2,3)

Swope, W. C., J. W. Pitera and F. Suits: Describing protein folding kinetics by molecular dynamics simulations: 1. Theory J. Phys. Chem. B 108, 6571-6581 (2004)

3(1,2)

Singhal, N., C. D. Snow, V. S. Pande: Using path sampling to build better Markovian state models: Predicting the folding rate and mechanism of a tryptophan zipper beta hairpin. J. Chem. Phys. 121, 415 (2004).

4

Sriraman, S., I. G. Kevrekidis and G. Hummer, G. J. Phys. Chem. B 109, 6479-6484 (2005)

5

Noe, F.: Probability Distributions of Molecular Observables computed from Markov Models. J. Chem. Phys. 128, 244103 (2008)

6

Buchete, N.-V. and Hummer, G.: Coarse master equations for peptide folding dynamics. J. Phys. Chem. B 112, 6057–6069 (2008)

7

Bowman, G. R., K. A. Beauchamp, G. Boxer and V. S. Pande: Progress and challenges in the automated construction of Markov state models for full protein systems. J. Chem. Phys. 131, 124101 (2009)

8

Sarich, M., F. Noe and C. Schuette: On the approximation quality of Markov state models. SIAM Multiscale Model. Simul. 8, 1154-1177 (2010)

9(1,2)

Prinz, J.-H., H. Wu, M. Sarich, B. Keller, M. Senne, M. Held, J. D. Chodera, C. Schuette and F. Noe: Markov models of molecular kinetics: Generation and Validation J. Chem. Phys. 134, 174105 (2011)

10

Bowman, G. R., V. S. Pande and F. Noe: An Introduction to Markov State Models and Their Application to Long Timescale Molecular Simulation. Advances in Experimental Medicine and Biology 797, Springer, Heidelberg (2014)

11

Nueske, F., Wu, H., Prinz, J.-H., Wehmeyer, C., Clementi, C. and Noe, F.: Markov State Models from short non-Equilibrium Simulations - Analysis and

Correction of Estimation Bias J. Chem. Phys. (submitted) (2017)

12

H. Wu and F. Noe: Variational approach for learning Markov processes from time series data (in preparation)

13(1,2)

Noe, F. and F. Nueske: A variational approach to modeling slow processes in stochastic dynamical systems. SIAM Multiscale Model. Simul. 11, 635-655 (2013).

14(1,2,3)

Wu, H and F. Noe: Variational approach for learning Markov processes from time series data (in preparation)

15

McGibbon, R and V. S. Pande: Variational cross-validation of slow dynamical modes in molecular kinetics, J. Chem. Phys. 142, 124105 (2015)

16

Noe, F. and C. Clementi: Kinetic distance and kinetic maps from molecular dynamics simulation. J. Chem. Theory Comput. 11, 5002-5011 (2015)

Example

>>> from pyemma import msm
>>> import numpy as np
>>> np.set_printoptions(precision=3)
>>> dtrajs = [[0,1,2,2,2,2,1,2,2,2,1,0,0,0,0,0,0,0], [0,0,0,0,1,1,2,2,2,2,2,2,2,1,0,0]]  # two trajectories
>>> mm = msm.estimate_markov_model(dtrajs, 2)

Which is the active set of states we are working on?

>>> print(mm.active_set)
[0 1 2]

Show the count matrix

>>> print(mm.count_matrix_active)
[[ 7.  2.  1.]
 [ 2.  0.  4.]
 [ 2.  3.  9.]]

Show the estimated transition matrix

>>> print(mm.transition_matrix)
[[ 0.7    0.167  0.133]
 [ 0.388  0.     0.612]
 [ 0.119  0.238  0.643]]

Is this model reversible (i.e. does it fulfill detailed balance)?

>>> print(mm.is_reversible)
True

What is the equilibrium distribution of states?

>>> print(mm.stationary_distribution)
[ 0.393  0.17   0.437]

Relaxation timescales?

>>> print(mm.timescales())
[ 3.415  1.297]

Mean first passage time from state 0 to 2:

>>> print(mm.mfpt(0, 2))  
9.929...