pyemma.msm.bayesian_markov_model¶
-
pyemma.msm.
bayesian_markov_model
(dtrajs, lag, reversible=True, statdist=None, sparse=False, connectivity='largest', count_mode='effective', nsamples=100, conf=0.95, dt_traj='1 step', show_progress=True, mincount_connectivity='1/n')¶ Bayesian Markov model estimate using Gibbs sampling of the posterior
Returns a
BayesianMSM
that contains the estimated transition matrix and allows to compute a large number of quantities related to Markov models as well as their statistical uncertainties.- 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) – lagtime for the MSM estimation in multiples of trajectory steps
reversible (bool, optional, default = True) – If true compute reversible MSM, else non-reversible MSM
sparse (bool, optional, default = False) – 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.
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)\]
connectivity (str, optional, default = None) –
Defines if the resulting HMM will be defined on all hidden states or on a connected subset. Connectivity is defined by counting only transitions with at least mincount_connectivity counts. If a subset of states is used, all estimated quantities (transition matrix, stationary distribution, etc) are only defined on this subset and are correspondingly smaller than nstates. Following modes are available: * None or ‘all’ : The active set is the full set of states.
Estimation is done on all weakly connected subsets separately. The resulting transition matrix may be disconnected.
’largest’ : The active set is the largest reversibly connected set.
- ’populous’The active set is the reversibly connected set with
most counts.
nsample (int, optional, default=100) – number of transition matrix samples to compute and store
conf (float, optional, default=0.95) – size of confidence intervals
dt_traj (str, optional, default='1 step') –
Description of the physical time corresponding to the trajectory time step. 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*’show_progress (bool, default=True) – Show progressbars for calculation
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.
- Returns
An
BayesianMSM
object containing the Bayesian MSM estimatorand the model.
Example
Note that the following example is only qualitatively and not quantitatively reproducible because it involves random numbers.
We build a Bayesian Markov model for the following two trajectories at lag time 2:
>>> from pyemma import msm >>> 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]] >>> mm = msm.bayesian_markov_model(dtrajs, 2, show_progress=False)
The resulting Model is an MSM just like you get with estimate_markov_model Its transition matrix does also come from a maximum likelihood estimation, but it’s slightly different from the estimate_markov_mode result because bayesian_markov_model uses an effective count matrix with statistically uncorrelated counts:
>>> print(mm.transition_matrix) # doctest: +SKIP [[ 0.70000001 0.16463699 0.135363 ] [ 0.38169055 0. 0.61830945] [ 0.12023989 0.23690297 0.64285714]]
However bayesian_markov_model returns a SampledMSM object which is able to compute the probability distribution and statistical models of all methods that are offered by the MSM object. This works as follows. You can ask for the sample mean and specify the method you wanna evaluate as a string:
>>> print(mm.sample_mean('transition_matrix')) # doctest: +SKIP [[ 0.71108663 0.15947371 0.12943966] [ 0.41076105 0. 0.58923895] [ 0.13079372 0.23005443 0.63915185]]
Likewise, the standard deviation by element:
>>> print(mm.sample_std('transition_matrix')) # doctest: +SKIP [[ 0.13707029 0.09479627 0.09200214] [ 0.15247454 0. 0.15247454] [ 0.07701315 0.09385258 0.1119089 ]]
And this is the 95% (2 sigma) confidence interval. You can control the percentile using the conf argument in this function:
>>> L, R = mm.sample_conf('transition_matrix') >>> print(L) # doctest: +SKIP >>> print(R) # doctest: +SKIP [[ 0.44083423 0.03926518 0.0242113 ] [ 0.14102544 0. 0.30729828] [ 0.02440188 0.07629456 0.43682481]] [[ 0.93571706 0.37522581 0.40180041] [ 0.69307665 0. 0.8649215 ] [ 0.31029752 0.44035732 0.85994006]]
If you wanna compute expectations of functions that require arguments, just pass these arguments as well:
>>> print(mm.sample_std('mfpt', 0, 2)) # doctest: +SKIP 12.9049811296
And if you want to histogram the distribution or compute more complex statistical moment such as the covariance between different quantities, just get the full sample of your quantity of interest and evaluate it at will:
>>> samples = mm.sample_f('mfpt', 0, 2) >>> print(samples[:4]) # doctest: +SKIP [7.9763615793248155, 8.6540958274695701, 26.295326015231058, 17.909895469938899]
Internally, the SampledMSM object has 100 transition matrices (the number can be controlled by nsamples), that were computed by the transition matrix sampling method. All of the above sample functions iterate over these 100 transition matrices and evaluate the requested function with the given parameters on each of them.
-
class
pyemma.msm.estimators.bayesian_msm.
BayesianMSM
(lag=1, nsamples=100, nsteps=None, reversible=True, statdist_constraint=None, count_mode='effective', sparse=False, connectivity='largest', dt_traj='1 step', conf=0.95, show_progress=True, mincount_connectivity='1/n')¶ Bayesian Markov state model estimator
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, **kw)- 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.
sample_conf
(f, *args, **kwargs)Sample confidence interval of numerical method f over all samples
sample_f
(f, *args, **kwargs)Evaluated method f for all samples
sample_mean
(f, *args, **kwargs)Sample mean of numerical method f over all samples
sample_std
(f, *args, **kwargs)Sample standard deviation of numerical method f over all samples
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
([samples, conf, P, pi, …])- param samples
sampled MSMs
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
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
The transition matrix on the active set.
The fraction of counts in the largest connected set.
The active set of states on which all computations and estimations will be done
The fraction of states in the largest connected set.
Ensures that the connected states are indexed and returns the indices
The reversible connected sets of states, sorted by size (descending)
Returns the connectivity mode of the MSM
The count matrix on the active set given the connectivity mode used.
The count matrix on full set of discrete states, irrespective as to whether they are connected or not.
A list of integer arrays with the discrete trajectories mapped to the connectivity mode used.
A list of integer arrays with the original (unmapped) discrete trajectories:
Description of the physical time corresponding to the lag.
A list of integer arrays with the discrete trajectories mapped to the connectivity mode used.
A list of integer arrays with the original (unmapped) discrete trajectories:
Statistically uncorrelated transition counts within the active set of states
Returns whether the MSM is reversible
Returns whether the MSM is sparse
The lag time at which the Markov model was estimated
The lag time at which the Markov model was estimated
The largest reversible connected set of states
The logger for this class instance
Assignment of states to metastable sets using PCCA++
Probability of metastable states to visit an MSM state by PCCA++
Probabilities of MSM states to belong to a metastable state by PCCA++
Metastable sets using PCCA++
The model estimated by this Estimator
Returns number of jobs/threads to use during assignment of data.
Number of states chosen for PCCA++ computation.
The name of this instance
number of eigenvalues to compute.
Number of active states on which all computations and estimations are done
Number of states in discrete trajectories
The stationary distribution on the MSM states
Returns whether the MSM is reversible
whether to show the progress of heavy calculations on this object.
Returns whether the MSM is sparse
The stationary distribution on the MSM states
Physical time corresponding to one transition matrix step, e.g.
The transition matrix on the active set.
-
P
¶ The transition matrix on the active set.
-
active_count_fraction
¶ The fraction of counts in the largest connected set.
-
active_set
¶ The active set of states on which all computations and estimations will be done
-
active_state_fraction
¶ The fraction of states in the largest connected set.
-
active_state_indexes
¶ Ensures that the connected states are indexed and returns the indices
-
cktest
(nsets, memberships=None, mlags=10, conf=0.95, err_est=False, n_jobs=None, show_progress=True)¶ Conducts a Chapman-Kolmogorow test.
- Parameters
nsets (int) – number of sets to test on
memberships (ndarray(nstates, nsets), optional) – optional state memberships. By default (None) will conduct a cktest on PCCA (metastable) sets.
mlags (int or int-array, optional) – multiples of lag times for testing the Model, e.g. range(10). A single int will trigger a range, i.e. mlags=10 maps to mlags=range(10). The setting None will choose mlags automatically according to the longest available trajectory
conf (float, optional) – confidence interval
err_est (bool, optional) – compute errors also for all estimations (computationally expensive) If False, only the prediction will get error bars, which is often sufficient to validate a model.
n_jobs (int, default=None) – how many jobs to use during calculation
show_progress (bool, optional) – Show progress bars for calculation?
- Returns
cktest
- Return type
References
This test was suggested in [1]_ and described in detail in [2]_.
- 1
F. Noe, Ch. Schuette, E. Vanden-Eijnden, L. Reich and T. Weikl: Constructing the Full Ensemble of Folding Pathways from Short Off-Equilibrium Simulations. Proc. Natl. Acad. Sci. USA, 106, 19011-19016 (2009)
- 2
Prinz, J H, H Wu, M Sarich, B Keller, M Senne, M Held, J D Chodera, C Schuette and F Noe. 2011. Markov models of molecular kinetics: Generation and validation. J Chem Phys 134: 174105
-
coarse_grain
(ncoarse, method='hmm')¶ Returns a coarse-grained Markov model.
Currently only the HMM method described in [1]_ is available for coarse-graining MSMs.
- Parameters
ncoarse (int) – number of coarse states
- Returns
hmsm
- Return type
MaximumLikelihoodHMSM
References
- 1
F. Noe, H. Wu, J.-H. Prinz and N. Plattner: Projected and hidden Markov models for calculating kinetics and metastable states of complex molecules J. Chem. Phys. 139, 184114 (2013)
-
committor_backward
(A, B)¶ Backward committor from set A to set B
- Parameters
A (int or int array) – set of starting states
B (int or int array) – set of target states
-
committor_forward
(A, B)¶ Forward committor (also known as p_fold or splitting probability) from set A to set B
- Parameters
A (int or int array) – set of starting states
B (int or int array) – set of target states
-
connected_sets
¶ The reversible connected sets of states, sorted by size (descending)
-
connectivity
¶ Returns the connectivity mode of the MSM
-
correlation
(a, b=None, maxtime=None, k=None, ncv=None)¶ Time-correlation for equilibrium experiment.
In order to simulate a time-correlation experiment (e.g. fluorescence correlation spectroscopy [NDD11], dynamical neutron scattering [LYP13], …), first compute the mean values of your experimental observable \(a\) by MSM state:
\[a_i = \frac{1}{N_i} \sum_{x_t \in S_i} f(x_t)\]where \(S_i\) is the set of configurations belonging to MSM state \(i\) and \(f()\) is a function that computes the experimental observable of interest for configuration \(x_t\). If a cross-correlation function is wanted, also apply the above computation to a second experimental observable \(b\).
Then the accurate (i.e. without statistical error) autocorrelation function of \(f(x_t)\) given the Markov model is computed by correlation(a), and the accurate cross-correlation function is computed by correlation(a,b). This is done by evaluating the equation
\[\begin{split}acf_a(k\tau) &= \mathbf{a}^\top \mathrm{diag}(\boldsymbol{\pi}) \mathbf{P(\tau)}^k \mathbf{a} \\ ccf_{a,b}(k\tau) &= \mathbf{a}^\top \mathrm{diag}(\boldsymbol{\pi}) \mathbf{P(\tau)}^k \mathbf{b}\end{split}\]where \(acf\) stands for autocorrelation function and \(ccf\) stands for cross-correlation function, \(\mathbf{P(\tau)}\) is the transition matrix at lag time \(\tau\), \(\boldsymbol{\pi}\) is the equilibrium distribution of \(\mathbf{P}\), and \(k\) is the time index.
Note that instead of using this method you could generate a long synthetic trajectory from the MSM and then estimating the time-correlation of your observable(s) directly from this trajectory. However, there is no reason to do this because the present method does that calculation without any sampling, and only in the limit of an infinitely long synthetic trajectory the two results will agree exactly. The correlation function computed by the present method still has statistical uncertainty from the fact that the underlying MSM transition matrix has statistical uncertainty when being estimated from data, but there is no additional (and unnecessary) uncertainty due to synthetic trajectory generation.
- Parameters
a ((n,) ndarray) – Observable, represented as vector on state space
maxtime (int or float) – Maximum time (in units of the input trajectory time step) until which the correlation function will be evaluated. Internally, the correlation function can only be computed in integer multiples of the Markov model lag time, and therefore the actual last time point will be computed at \(\mathrm{ceil}(\mathrm{maxtime} / \tau)\) By default (None), the maxtime will be set equal to the 5 times the slowest relaxation time of the MSM, because after this time the signal is almost constant.
b ((n,) ndarray (optional)) – Second observable, for cross-correlations
k (int (optional)) – Number of eigenvalues and eigenvectors to use for computation. This option is only relevant for sparse matrices and long times for which an eigenvalue decomposition will be done instead of using the matrix power.
ncv (int (optional)) – Only relevant for sparse matrices and large lag times where the relaxation will be computed using an eigenvalue decomposition. The number of Lanczos vectors generated, ncv must be greater than k; it is recommended that ncv > 2*k.
- Returns
times (ndarray (N)) – Time points (in units of the input trajectory time step) at which the correlation has been computed
correlations (ndarray (N)) – Correlation values at given times
Examples
This example computes the autocorrelation function of a simple observable on a three-state Markov model and plots the result using matplotlib:
>>> import numpy as np >>> import pyemma.msm as msm >>> >>> P = np.array([[0.99, 0.01, 0], [0.01, 0.9, 0.09], [0, 0.1, 0.9]]) >>> a = np.array([0.0, 0.5, 1.0]) >>> M = msm.markov_model(P) >>> times, acf = M.correlation(a) >>> >>> import matplotlib.pylab as plt >>> plt.plot(times, acf) # doctest: +SKIP
References
- NDD11
Noe, F., S. Doose, I. Daidone, M. Loellmann, J. D. Chodera, M. Sauer and J. C. Smith. 2011 Dynamical fingerprints for probing individual relaxation processes in biomolecular dynamics with simulations and kinetic experiments. Proc. Natl. Acad. Sci. USA 108, 4822-4827.
- LYP13
Lindner, B., Z. Yi, J.-H. Prinz, J. C. Smith and F. Noe. 2013. Dynamic Neutron Scattering from Conformational Dynamics I: Theory and Markov models. J. Chem. Phys. 139, 175101.
-
count_matrix_active
¶ The count matrix on the active set given the connectivity mode used.
For example, for connectivity=’largest’, the count matrix is given only on the largest reversibly connected set. Attention: This count matrix has been obtained by sliding a window of length tau across the data. It contains a factor of tau more counts than are statistically uncorrelated. It’s fine to use this matrix for maximum likelihood estimated, but it will give far too small errors if you use it for uncertainty calculations. In order to do uncertainty calculations, use the effective count matrix, see:
effective_count_matrix
See also
effective_count_matrix
For a count matrix with effective (statistically uncorrelated) counts.
-
count_matrix_full
¶ The count matrix on full set of discrete states, irrespective as to whether they are connected or not. Attention: This count matrix has been obtained by sliding a window of length tau across the data. It contains a factor of tau more counts than are statistically uncorrelated. It’s fine to use this matrix for maximum likelihood estimated, but it will give far too small errors if you use it for uncertainty calculations. In order to do uncertainty calculations, use the effective count matrix, see:
effective_count_matrix
(only implemented on the active set), or divide this count matrix by tau.See also
effective_count_matrix
For a active-set count matrix with effective (statistically uncorrelated) counts.
-
discrete_trajectories_active
¶ A list of integer arrays with the discrete trajectories mapped to the connectivity mode used. For example, for connectivity=’largest’, the indexes will be given within the connected set. Frames that are not in the connected set will be -1.
-
discrete_trajectories_full
¶ A list of integer arrays with the original (unmapped) discrete trajectories:
-
dt_model
¶ Description of the physical time corresponding to the lag.
-
dt_traj
¶
-
dtrajs_active
¶ A list of integer arrays with the discrete trajectories mapped to the connectivity mode used. For example, for connectivity=’largest’, the indexes will be given within the connected set. Frames that are not in the connected set will be -1.
-
dtrajs_full
¶ A list of integer arrays with the original (unmapped) discrete trajectories:
-
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
-
eigenvalues
(k=None)¶ Compute the transition matrix eigenvalues
- Parameters
k (int) – number of eigenvalues to be returned. By default will return all available eigenvalues
- Returns
ts – transition matrix eigenvalues \(\lambda_i, i = 1, ..., k\)., sorted by descending norm.
- Return type
ndarray(k,)
-
eigenvectors_left
(k=None)¶ Compute the left transition matrix eigenvectors
- Parameters
k (int) – number of eigenvectors to be returned. By default all available eigenvectors.
- Returns
L – left eigenvectors in a row matrix. l_ij is the j’th component of the i’th left eigenvector
- Return type
ndarray(k,n)
-
eigenvectors_right
(k=None)¶ Compute the right transition matrix eigenvectors
- Parameters
k (int) – number of eigenvectors to be computed. By default all available eigenvectors.
- Returns
R – right eigenvectors in a column matrix. r_ij is the i’th component of the j’th right eigenvector
- Return type
ndarray(n,k)
-
estimate
(dtrajs, **kw)¶ - 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.
kw (dict) – Other Parameters. See documentation of class.
- Returns
msm – Estimated Hidden Markov state model
- Return type
-
expectation
(a)¶ Equilibrium expectation value of a given observable.
- Parameters
a ((n,) ndarray) – Observable vector on the MSM state space
- Returns
val – Equilibrium expectation value fo the given observable
- Return type
float
Notes
The equilibrium expectation value of an observable \(a\) is defined as follows
\[\mathbb{E}_{\mu}[a] = \sum_i \pi_i a_i\]\(\pi=(\pi_i)\) is the stationary vector of the transition matrix \(P\).
-
fingerprint_correlation
(a, b=None, k=None, ncv=None)¶ Dynamical fingerprint for equilibrium time-correlation experiment.
- Parameters
a ((n,) ndarray) – Observable, represented as vector on MSM state space
b ((n,) ndarray, optional) – Second observable, for cross-correlations
k (int, optional) – Number of eigenvalues and eigenvectors to use for computation. This option is only relevant for sparse matrices and long times for which an eigenvalue decomposition will be done instead of using the matrix power
ncv (int, optional) – Only relevant for sparse matrices and large lag times, where the relaxation will be computed using an eigenvalue decomposition. The number of Lanczos vectors generated, ncv must be greater than k; it is recommended that ncv > 2*k
- Returns
timescales ((k,) ndarray) – Time-scales (in units of the input trajectory time step) of the transition matrix
amplitudes ((k,) ndarray) – Amplitudes for the correlation experiment
References
Spectral densities are commonly used in spectroscopy. Dynamical fingerprints are a useful representation for computational spectroscopy results and have been introduced in [NDD11].
- NDD11
Noe, F, S Doose, I Daidone, M Loellmann, M Sauer, J D Chodera and J Smith. 2010. Dynamical fingerprints for probing individual relaxation processes in biomolecular dynamics with simulations and kinetic experiments. PNAS 108 (12): 4822-4827.
-
fingerprint_relaxation
(p0, a, k=None, ncv=None)¶ Dynamical fingerprint for perturbation/relaxation experiment.
- Parameters
p0 ((n,) ndarray) – Initial distribution for a relaxation experiment
a ((n,) ndarray) – Observable, represented as vector on state space
k (int, optional) – Number of eigenvalues and eigenvectors to use for computation
ncv (int, optional) – Only relevant for sparse matrices and large lag times, where the relaxation will be computes using an eigenvalue decomposition. The number of Lanczos vectors generated, ncv must be greater than k; it is recommended that ncv > 2*k
- Returns
timescales ((k,) ndarray) – Time-scales (in units of the input trajectory time step) of the transition matrix
amplitudes ((k,) ndarray) – Amplitudes for the relaxation experiment
References
Spectral densities are commonly used in spectroscopy. Dynamical fingerprints are a useful representation for computational spectroscopy results and have been introduced in [NDD11].
- NDD11
Noe, F, S Doose, I Daidone, M Loellmann, M Sauer, J D Chodera and J Smith. 2010. Dynamical fingerprints for probing individual relaxation processes in biomolecular dynamics with simulations and kinetic experiments. PNAS 108 (12): 4822-4827.
-
fit
(X, y=None)¶ Estimates parameters - for compatibility with sklearn.
- Parameters
X (object) – A reference to the data from which the model will be estimated
- Returns
estimator – The estimator (self) with estimated model.
- Return type
object
-
generate_traj
(N, start=None, stop=None, stride=1)¶ Generates a synthetic discrete trajectory of length N and simulation time stride * lag time * N
This information can be used in order to generate a synthetic molecular dynamics trajectory - see
pyemma.coordinates.save_traj()
Note that the time different between two samples is the Markov model lag time tau. When comparing quantities computing from this synthetic trajectory and from the input trajectories, the time points of this trajectory must be scaled by the lag time in order to have them on the same time scale.
- Parameters
N (int) – Number of time steps in the output trajectory. The total simulation time is stride * lag time * N
start (int, optional, default = None) – starting state. If not given, will sample from the stationary distribution of P
stop (int or int-array-like, optional, default = None) – stopping set. If given, the trajectory will be stopped before N steps once a state of the stop set is reached
stride (int, optional, default = 1) – Multiple of lag time used as a time step. By default, the time step is equal to the lag time
- Returns
indexes – trajectory and time indexes of the simulated trajectory. Each row consist of a tuple (i, t), where i is the index of the trajectory and t is the time index within the trajectory. Note that the time different between two samples is the Markov model lag time tau
- Return type
ndarray( (N, 2) )
See also
pyemma.coordinates.save_traj()
in order to save this synthetic trajectory as a trajectory file with molecular structures
-
get_model_params
(deep=True)¶ Get parameters for this model.
- Parameters
deep (boolean, optional) – If True, will return the parameters for this estimator and contained subobjects that are estimators.
- Returns
params – Parameter names mapped to their values.
- Return type
mapping of string to any
-
get_params
(deep=True)¶ Get parameters for this estimator.
- Parameters
deep (boolean, optional) – If True, will return the parameters for this estimator and contained subobjects that are estimators.
- Returns
params – Parameter names mapped to their values.
- Return type
mapping of string to any
-
hmm
(nhidden)¶ Estimates a hidden Markov state model as described in [1]_
- Parameters
nhidden (int) – number of hidden (metastable) states
- Returns
hmsm
- Return type
MaximumLikelihoodHMSM
References
- 1
F. Noe, H. Wu, J.-H. Prinz and N. Plattner: Projected and hidden Markov models for calculating kinetics and metastable states of complex molecules J. Chem. Phys. 139, 184114 (2013)
-
is_reversible
¶ Returns whether the MSM is reversible
-
is_sparse
¶ Returns whether the MSM is sparse
-
lag
¶ The lag time at which the Markov model was estimated
-
lagtime
¶ The lag time at which the Markov model was estimated
-
largest_connected_set
¶ The largest reversible connected set of states
-
classmethod
load
(file_name, model_name='default')¶ Loads a previously saved PyEMMA object from disk.
- Parameters
file_name (str or file like object (has to provide read method)) – The file like object tried to be read for a serialized object.
model_name (str, default='default') – if multiple models are contained in the file, these can be accessed by their name. Use
pyemma.list_models()
to get a representation of all stored models.
- Returns
obj
- Return type
the de-serialized object
-
logger
¶ The logger for this class instance
-
metastable_assignments
¶ Assignment of states to metastable sets using PCCA++
Computes the assignment to metastable sets for active set states using the PCCA++ method [1]_.
pcca()
needs to be called first to make this attribute available.This is only recommended for visualization purposes. You cannot compute any actual quantity of the coarse-grained kinetics without employing the fuzzy memberships!
- Returns
assignments – For each MSM state, the metastable state it is located in.
- Return type
ndarray (n,)
See also
pcca
to compute the metastable decomposition
References
- 1
Roeblitz, S and M Weber. 2013. Fuzzy spectral clustering by PCCA+: application to Markov state models and data classification. Advances in Data Analysis and Classification 7 (2): 147-179
-
metastable_distributions
¶ Probability of metastable states to visit an MSM state by PCCA++
Computes the probability distributions of active set states within each metastable set by combining the the PCCA++ method [1]_ with Bayesian inversion as described in [2]_.
pcca()
needs to be called first to make this attribute available.- Returns
p_out – A matrix containing the probability distribution of each active set state, given that we are in one of the m metastable sets, i.e. p(state | metastable). The row sums of p_out are 1.
- Return type
ndarray (m,n)
See also
pcca
to compute the metastable decomposition
References
- 1
Roeblitz, S and M Weber. 2013. Fuzzy spectral clustering by PCCA+: application to Markov state models and data classification. Advances in Data Analysis and Classification 7, 147-179.
- 2
F. Noe, H. Wu, J.-H. Prinz and N. Plattner. 2013. Projected and hidden Markov models for calculating kinetics and metastable states of complex molecules J. Chem. Phys. 139, 184114.
-
metastable_memberships
¶ Probabilities of MSM states to belong to a metastable state by PCCA++
Computes the memberships of active set states to metastable sets with the PCCA++ method [1]_.
pcca()
needs to be called first to make this attribute available.- Returns
M – A matrix containing the probability or membership of each state to be assigned to each metastable set, i.e. p(metastable | state). The row sums of M are 1.
- Return type
ndarray((n,m))
See also
pcca
to compute the metastable decomposition
References
- 1
Roeblitz, S and M Weber. 2013. Fuzzy spectral clustering by PCCA+: application to Markov state models and data classification. Advances in Data Analysis and Classification 7 (2): 147-179
-
metastable_sets
¶ Metastable sets using PCCA++
Computes the metastable sets of active set states within each metastable set using the PCCA++ method [1]_.
pcca()
needs to be called first to make this attribute available.This is only recommended for visualization purposes. You cannot compute any actual quantity of the coarse-grained kinetics without employing the fuzzy memberships!
- Returns
sets – A list of length equal to metastable states. Each element is an array with microstate indexes contained in it
- Return type
list of ndarray
See also
pcca
to compute the metastable decomposition
References
- 1
Roeblitz, S and M Weber. 2013. Fuzzy spectral clustering by PCCA+: application to Markov state models and data classification. Advances in Data Analysis and Classification 7 (2): 147-179
-
mfpt
(A, B)¶ Mean first passage times from set A to set B, in units of the input trajectory time step
- Parameters
A (int or int array) – set of starting states
B (int or int array) – set of target states
-
model
¶ The model estimated by this Estimator
-
n_jobs
¶ Returns number of jobs/threads to use during assignment of data.
- Returns
If None it will return the setting of ‘PYEMMA_NJOBS’ or
’SLURM_CPUS_ON_NODE’ environment variable. If none of these environment variables exist,
the number of processors /or cores is returned.
Notes
This setting will effectively be multiplied by the the number of threads used by NumPy for algorithms which use multiple processes. So take care if you choose this manually.
-
n_metastable
¶ Number of states chosen for PCCA++ computation.
-
name
¶ The name of this instance
-
neig
¶ number of eigenvalues to compute.
-
nstates
¶ Number of active states on which all computations and estimations are done
-
nstates_full
¶ Number of states in discrete trajectories
-
pcca
(m)¶ Runs PCCA++ [1]_ to compute a metastable decomposition of MSM states
After calling this method you can access
metastable_memberships()
,metastable_distributions()
,metastable_sets()
andmetastable_assignments()
.- Parameters
m (int) – Number of metastable sets
- Returns
pcca_obj – An object containing all PCCA quantities. However, you can also ignore this return value and instead retrieve the quantities of your interest with the following MSM functions:
metastable_memberships()
,metastable_distributions()
,metastable_sets()
andmetastable_assignments()
.- Return type
Notes
If you coarse grain with PCCA++, the order of the obtained memberships might not be preserved. This also applies for
metastable_memberships()
,metastable_distributions()
,metastable_sets()
,metastable_assignments()
References
- 1
Roeblitz, S and M Weber. 2013. Fuzzy spectral clustering by PCCA+: application to Markov state models and data classification. Advances in Data Analysis and Classification 7 (2): 147-179
-
pi
¶ The stationary distribution on the MSM states
-
propagate
(p0, k)¶ Propagates the initial distribution p0 k times
Computes the product
\[p_k = p_0^T P^k\]If the lag time of transition matrix \(P\) is \(\tau\), this will provide the probability distribution at time \(k \tau\).
- Parameters
p0 (ndarray(n,)) – Initial distribution. Vector of size of the active set.
k (int) – Number of time steps
- Returns
pk – Distribution after k steps. Vector of size of the active set.
- Return type
ndarray(n,)
-
relaxation
(p0, a, maxtime=None, k=None, ncv=None)¶ Simulates a perturbation-relaxation experiment.
In perturbation-relaxation experiments such as temperature-jump, pH-jump, pressure jump or rapid mixing experiments, an ensemble of molecules is initially prepared in an off-equilibrium distribution and the expectation value of some experimental observable is then followed over time as the ensemble relaxes towards equilibrium.
In order to simulate such an experiment, first determine the distribution of states at which the experiment is started, \(p_0\) and compute the mean values of your experimental observable \(a\) by MSM state:
\[a_i = \frac{1}{N_i} \sum_{x_t \in S_i} f(x_t)\]where \(S_i\) is the set of configurations belonging to MSM state \(i\) and \(f()\) is a function that computes the experimental observable of interest for configuration \(x_t\).
Then the accurate (i.e. without statistical error) time-dependent expectation value of \(f(x_t)\) given the Markov model is computed by relaxation(p0, a). This is done by evaluating the equation
\[E_a(k\tau) = \mathbf{p_0}^{\top} \mathbf{P(\tau)}^k \mathbf{a}\]where \(E\) stands for the expectation value that relaxes to its equilibrium value that is identical to expectation(a), \(\mathbf{P(\tau)}\) is the transition matrix at lag time \(\tau\), \(\boldsymbol{\pi}\) is the equilibrium distribution of \(\mathbf{P}\), and \(k\) is the time index.
Note that instead of using this method you could generate many synthetic trajectory from the MSM with starting points drawn from the initial distribution and then estimating the time-dependent expectation value by an ensemble average. However, there is no reason to do this because the present method does that calculation without any sampling, and only in the limit of an infinitely many trajectories the two results will agree exactly. The relaxation function computed by the present method still has statistical uncertainty from the fact that the underlying MSM transition matrix has statistical uncertainty when being estimated from data, but there is no additional (and unnecessary) uncertainty due to synthetic trajectory generation.
- Parameters
p0 ((n,) ndarray) – Initial distribution for a relaxation experiment
a ((n,) ndarray) – Observable, represented as vector on state space
maxtime (int or float, optional) – Maximum time (in units of the input trajectory time step) until which the correlation function will be evaluated. Internally, the correlation function can only be computed in integer multiples of the Markov model lag time, and therefore the actual last time point will be computed at \(\mathrm{ceil}(\mathrm{maxtime} / \tau)\). By default (None), the maxtime will be set equal to the 5 times the slowest relaxation time of the MSM, because after this time the signal is constant.
k (int (optional)) – Number of eigenvalues and eigenvectors to use for computation
ncv (int (optional)) – Only relevant for sparse matrices and large lag times, where the relaxation will be computes using an eigenvalue decomposition. The number of Lanczos vectors generated, ncv must be greater than k; it is recommended that ncv > 2*k
- Returns
times (ndarray (N)) – Time points (in units of the input trajectory time step) at which the relaxation has been computed
res (ndarray) – Array of expectation value at given times
-
reversible
¶ Returns whether the MSM is reversible
-
sample_by_distributions
(distributions, nsample)¶ Generates samples according to given probability distributions
- Parameters
distributions (list or array of ndarray ( (n) )) – m distributions over states. Each distribution must be of length n and must sum up to 1.0
nsample (int) – Number of samples per distribution. If replace = False, the number of returned samples per state could be smaller if less than nsample indexes are available for a state.
- Returns
indexes – List of the sampled indices by distribution. Each element is an index array with a number of rows equal to nsample, with rows consisting of a tuple (i, t), where i is the index of the trajectory and t is the time index within the trajectory.
- Return type
length m list of ndarray( (nsample, 2) )
-
sample_by_state
(nsample, subset=None, replace=True)¶ Generates samples of the connected states.
For each state in the active set of states, generates nsample samples with trajectory/time indexes. This information can be used in order to generate a trajectory of length nsample * nconnected using
pyemma.coordinates.save_traj()
or nconnected trajectories of length nsample each usingpyemma.coordinates.save_traj()
- Parameters
nsample (int) – Number of samples per state. If replace = False, the number of returned samples per state could be smaller if less than nsample indexes are available for a state.
subset (ndarray((n)), optional, default = None) – array of states to be indexed. By default all states in the connected set will be used
replace (boolean, optional) – Whether the sample is with or without replacement
- Returns
indexes – list of trajectory/time index arrays with an array for each state. Within each index array, each row consist of a tuple (i, t), where i is the index of the trajectory and t is the time index within the trajectory.
- Return type
list of ndarray( (N, 2) )
See also
pyemma.coordinates.save_traj()
in order to save the sampled frames sequentially in a trajectory file with molecular structures
pyemma.coordinates.save_trajs()
in order to save the sampled frames in nconnected trajectory files with molecular structures
-
sample_conf
(f, *args, **kwargs)¶ Sample confidence interval of numerical method f over all samples
Calls f(*args, **kwargs) on all samples and computes the confidence interval. Size of confidence interval is given in the construction of the SampledModel. f must return a numerical value or an ndarray.
- Parameters
f (method reference or name (str)) – Model method to be evaluated for each model sample
args (arguments) – Non-keyword arguments to be passed to the method in each call
kwargs (keyword-argments) – Keyword arguments to be passed to the method in each call
- Returns
L (float or ndarray) – lower value or array of confidence interval
R (float or ndarray) – upper value or array of confidence interval
-
sample_f
(f, *args, **kwargs)¶ Evaluated method f for all samples
Calls f(*args, **kwargs) on all samples.
- Parameters
f (method reference or name (str)) – Model method to be evaluated for each model sample
args (arguments) – Non-keyword arguments to be passed to the method in each call
kwargs (keyword-argments) – Keyword arguments to be passed to the method in each call
- Returns
vals – list of results of the method calls
- Return type
list
-
sample_mean
(f, *args, **kwargs)¶ Sample mean of numerical method f over all samples
Calls f(*args, **kwargs) on all samples and computes the mean. f must return a numerical value or an ndarray.
- Parameters
f (method reference or name (str)) – Model method to be evaluated for each model sample
args (arguments) – Non-keyword arguments to be passed to the method in each call
kwargs (keyword-argments) – Keyword arguments to be passed to the method in each call
- Returns
mean – mean value or mean array
- Return type
float or ndarray
-
sample_std
(f, *args, **kwargs)¶ Sample standard deviation of numerical method f over all samples
Calls f(*args, **kwargs) on all samples and computes the standard deviation. f must return a numerical value or an ndarray.
- Parameters
f (method reference or name (str)) – Model method to be evaluated for each model sample
args (arguments) – Non-keyword arguments to be passed to the method in each call
kwargs (keyword-argments) – Keyword arguments to be passed to the method in each call
- Returns
std – standard deviation or array of standard deviations
- Return type
float or ndarray
-
samples
¶
-
save
(file_name, model_name='default', overwrite=False, save_streaming_chain=False)¶ saves the current state of this object to given file and name.
- Parameters
file_name (str) – path to desired output file
model_name (str, default='default') – creates a group named ‘model_name’ in the given file, which will contain all of the data. If the name already exists, and overwrite is False (default) will raise a RuntimeError.
overwrite (bool, default=False) – Should overwrite existing model names?
save_streaming_chain (boolean, default=False) – if True, the data_producer(s) of this object will also be saved in the given file.
Examples
>>> import pyemma, numpy as np >>> from pyemma.util.contexts import named_temporary_file >>> m = pyemma.msm.MSM(P=np.array([[0.1, 0.9], [0.9, 0.1]]))
>>> with named_temporary_file() as file: # doctest: +SKIP ... m.save(file, 'simple') # doctest: +SKIP ... inst_restored = pyemma.load(file, 'simple') # doctest: +SKIP >>> np.testing.assert_equal(m.P, inst_restored.P) # doctest: +SKIP
-
score
(dtrajs, score_method=None, score_k=None)¶ Scores the MSM using the dtrajs using the variational approach for Markov processes [1]_ [2]_
Currently only implemented using dense matrices - will be slow for large state spaces.
- Parameters
dtrajs (list of arrays) – test data (discrete trajectories).
score_method (str, optional, default='VAMP2') – Overwrite scoring method if desired. If None, the estimators scoring method will be used. See __init__ for documentation.
score_k (int or None) – Overwrite scoring rank if desired. If None, the estimators scoring rank will be used. See __init__ for documentation.
score_method –
Overwrite scoring method to be used if desired. If None, the estimators scoring method will be used. Available scores are based on the variational approach for Markov processes [1]_ [2]_ :
score_k – The maximum number of eigenvalues or singular values used in the score. If set to None, all available eigenvalues will be used.
References
- 1
Noe, F. and F. Nueske: A variational approach to modeling slow processes in stochastic dynamical systems. SIAM Multiscale Model. Simul. 11, 635-655 (2013).
- 2
Wu, H and F. Noe: Variational approach for learning Markov processes from time series data (in preparation)
- 3
McGibbon, R and V. S. Pande: Variational cross-validation of slow dynamical modes in molecular kinetics, J. Chem. Phys. 142, 124105 (2015)
- 4
Noe, F. and C. Clementi: Kinetic distance and kinetic maps from molecular dynamics simulation. J. Chem. Theory Comput. 11, 5002-5011 (2015)
-
score_cv
(dtrajs, n=10, score_method=None, score_k=None)¶ Scores the MSM using the variational approach for Markov processes [1]_ [2]_ and crossvalidation [3]_ .
Divides the data into training and test data, fits a MSM using the training data using the parameters of this estimator, and scores is using the test data. Currently only one way of splitting is implemented, where for each n, the data is randomly divided into two approximately equally large sets of discrete trajectory fragments with lengths of at least the lagtime.
Currently only implemented using dense matrices - will be slow for large state spaces.
- Parameters
dtrajs (list of arrays) – Test data (discrete trajectories).
n (number of samples) – Number of repetitions of the cross-validation. Use large n to get solid means of the score.
score_method (str, optional, default='VAMP2') –
Overwrite scoring method to be used if desired. If None, the estimators scoring method will be used. Available scores are based on the variational approach for Markov processes [1]_ [2]_ :
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.
References
- 1
Noe, F. and F. Nueske: A variational approach to modeling slow processes in stochastic dynamical systems. SIAM Multiscale Model. Simul. 11, 635-655 (2013).
- 2
Wu, H and F. Noe: Variational approach for learning Markov processes from time series data (in preparation).
- 3
McGibbon, R and V. S. Pande: Variational cross-validation of slow dynamical modes in molecular kinetics, J. Chem. Phys. 142, 124105 (2015).
- 4
Noe, F. and C. Clementi: Kinetic distance and kinetic maps from molecular dynamics simulation. J. Chem. Theory Comput. 11, 5002-5011 (2015).
-
set_model_params
(samples=None, conf=0.95, P=None, pi=None, reversible=None, dt_model='1 step', neig=None)¶ - Parameters
samples (list of MSM objects) – sampled MSMs
conf (float, optional, default=0.68) – Confidence interval. By default one-sigma (68.3%) is used. Use 95.4% for two sigma or 99.7% for three sigma.
-
set_params
(**params)¶ Set the parameters of this estimator. The method works on simple estimators as well as on nested objects (such as pipelines). The former have parameters of the form
<component>__<parameter>
so that it’s possible to update each component of a nested object. :returns: :rtype: self
-
show_progress
¶ whether to show the progress of heavy calculations on this object.
-
simulate
(N, start=None, stop=None, dt=1)¶ Generates a realization of the Markov Model
- Parameters
N (int) – trajectory length in steps of the lag time
start (int, optional, default = None) – starting hidden state. If not given, will sample from the stationary distribution of the hidden transition matrix.
stop (int or int-array-like, optional, default = None) – stopping hidden set. If given, the trajectory will be stopped before N steps once a hidden state of the stop set is reached
dt (int) – trajectory will be saved every dt time steps. Internally, the dt’th power of P is taken to ensure a more efficient simulation.
- Returns
htraj – The state trajectory with length N/dt
- Return type
(N/dt, ) ndarray
-
sparse
¶ Returns whether the MSM is sparse
-
stationary_distribution
¶ The stationary distribution on the MSM states
-
timescales
(k=None)¶ The relaxation timescales corresponding to the eigenvalues
- Parameters
k (int) – number of timescales to be returned. By default all available eigenvalues, minus 1.
- Returns
ts – relaxation timescales in units of the input trajectory time step, defined by \(-\tau / ln | \lambda_i |, i = 2,...,k+1\).
- Return type
ndarray(m)
-
timestep_model
¶ Physical time corresponding to one transition matrix step, e.g. ‘10 ps’
-
trajectory_weights
()¶ Uses the MSM to assign a probability weight to each trajectory frame.
This is a powerful function for the calculation of arbitrary observables in the trajectories one has started the analysis with. The stationary probability of the MSM will be used to reweigh all states. Returns a list of weight arrays, one for each trajectory, and with a number of elements equal to trajectory frames. Given \(N\) trajectories of lengths \(T_1\) to \(T_N\), this function returns corresponding weights:
\[(w_{1,1}, ..., w_{1,T_1}), (w_{N,1}, ..., w_{N,T_N})\]that are normalized to one:
\[\sum_{i=1}^N \sum_{t=1}^{T_i} w_{i,t} = 1\]Suppose you are interested in computing the expectation value of a function \(a(x)\), where \(x\) are your input configurations. Use this function to compute the weights of all input configurations and obtain the estimated expectation by:
\[\langle a \rangle = \sum_{i=1}^N \sum_{t=1}^{T_i} w_{i,t} a(x_{i,t})\]Or if you are interested in computing the time-lagged correlation between functions \(a(x)\) and \(b(x)\) you could do:
\[\langle a(t) b(t+\tau) \rangle_t = \sum_{i=1}^N \sum_{t=1}^{T_i} w_{i,t} a(x_{i,t}) a(x_{i,t+\tau})\]- Returns
weights – The normalized trajectory weights. Given \(N\) trajectories of lengths \(T_1\) to \(T_N\), returns the corresponding weights:
\[(w_{1,1}, ..., w_{1,T_1}), (w_{N,1}, ..., w_{N,T_N})\]- Return type
list of ndarray
-
transition_matrix
¶ The transition matrix on the active set.
-
update_model_params
(**params)¶ Update given model parameter if they are set to specific values
References
- 1
Trendelkamp-Schroer, B, H. Wu, F. Paul and F. Noe: Estimation and uncertainty of reversible Markov models. http://arxiv.org/abs/1507.05990