pyemma.msm.estimate_hidden_markov_model¶
-
pyemma.msm.
estimate_hidden_markov_model
(dtrajs, nstates, lag, reversible=True, stationary=False, connectivity=None, mincount_connectivity='1/n', separate=None, observe_nonempty=True, stride=1, dt_traj='1 step', accuracy=0.001, maxit=1000)¶ Estimates a Hidden Markov state model from discrete trajectories
Returns a
MaximumLikelihoodHMSM
that contains a transition matrix between a few (hidden) metastable states. Each metastable state has a probability distribution of visiting the discrete ‘microstates’ contained in the input trajectories. The resulting object is a hidden Markov model that allows to compute a large number of quantities.- 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
nstates (int) – the number of metastable states in the resulting HMM
reversible (bool, optional, default = True) – If true compute reversible MSM, else non-reversible MSM
stationary (bool, optional, default=False) – If True, the initial distribution of hidden states is self-consistently computed as the stationary distribution of the transition matrix. If False, it will be estimated from the starting states. Only set this to true if you’re sure that the observation trajectories are initiated from a global equilibrium distribution.
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.
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.
separate (None or iterable of int) – Force the given set of observed states to stay in a separate hidden state. The remaining nstates-1 states will be assigned by a metastable decomposition.
observe_nonempty (bool) – If True, will restricted the observed states to the states that have at least one observation in the lagged input trajectories.
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*’accuracy (float) – convergence threshold for EM iteration. When two the likelihood does not increase by more than accuracy, the iteration is stopped successfully.
maxit (int) – stopping criterion for EM iteration. When so many iterations are performed without reaching the requested accuracy, the iteration is stopped without convergence (a warning is given)
- Returns
hmsm – Estimator object containing the HMSM and estimation information.
- Return type
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_hidden_markov_model(dtrajs, 2, 2)
We have estimated a 2x2 hidden transition matrix between the metastable states:
>>> print(mm.transition_matrix) [[ 0.684 0.316] [ 0.242 0.758]]
With the equilibrium distribution:
>>> print(mm.stationary_distribution) # doctest: +ELLIPSIS [ 0.43... 0.56...]
The observed states are the three discrete clusters that we have in our discrete trajectory:
>>> print(mm.observable_set) [0 1 2]
The metastable distributions (mm.metastable_distributions), or equivalently the observation probabilities are the probability to be in a given cluster (‘microstate’) if we are in one of the hidden metastable states. So it’s a 2 x 3 matrix:
>>> print(mm.observation_probabilities) # doctest: +SKIP [[ 0.9620883 0.0379117 0. ] [ 0. 0.28014352 0.71985648]]
The first metastable state ist mostly in cluster 0, and a little bit in the transition state cluster 1. The second metastable state is less well defined, but mostly in cluster 2 and less prominently in the transition state cluster 1.
We can print the lifetimes of the metastable states:
>>> print(mm.lifetimes) # doctest: +ELLIPSIS [ 5... 7...]
And the timescale of the hidden transition matrix - now we only have one relaxation timescale:
>>> print(mm.timescales()) # doctest: +ELLIPSIS [ 2.4...]
The mean first passage times can also be computed between metastable states:
>>> print(mm.mfpt(0, 1)) # doctest: +ELLIPSIS 6.3...
-
class
pyemma.msm.estimators.maximum_likelihood_hmsm.
MaximumLikelihoodHMSM
(nstates=2, lag=1, stride=1, msm_init='largest-strong', reversible=True, stationary=False, connectivity=None, mincount_connectivity='1/n', observe_nonempty=True, separate=None, dt_traj='1 step', accuracy=0.001, maxit=1000)¶ Maximum likelihood estimator for a Hidden MSM given a MSM
Methods
cktest
([mlags, conf, err_est, n_jobs, …])Conducts a Chapman-Kolmogorow test.
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
(X, **params)Estimates the model given the data X
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.
get_model_params
([deep])Get parameters for this model.
get_params
([deep])Get parameters for this estimator.
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_observation_probabilities
(nsample)Generates samples according to the current observation probability distribution
save
(file_name[, model_name, overwrite, …])saves the current state of this object to given file and name.
set_model_params
([P, pobs, pi, reversible, …])- param P
coarse-grained or hidden transition matrix
set_params
(**params)Set the parameters of this estimator.
simulate
(N[, start, stop, dt])Generates a realization of the Hidden Markov Model
submodel
([states, obs, …])Returns a HMM with restricted state space
submodel_disconnect
([mincount_connectivity])Disconnects sets of hidden states that are barely connected
submodel_largest
([strong, mincount_connectivity])Returns the largest connected sub-HMM (convenience function)
submodel_populous
([strong, …])Returns the most populous connected sub-HMM (convenience function)
timescales
([k])The relaxation timescales corresponding to the eigenvalues
Uses the HMSM to assign a probability weight to each trajectory frame.
Computes the transition matrix between observed states
update_model_params
(**params)Update given model parameter if they are set to specific values
Attributes
The transition matrix on the active set.
The active set of hidden states on which all hidden state computations are done
A list of integer arrays with the original trajectories.
Transformed original trajectories that are used as an input into the HMM estimation
A list of integer arrays with the discrete trajectories mapped to the observation mode used.
Description of the physical time corresponding to the lag.
A list of integer arrays with the original trajectories.
Transformed original trajectories that are used as an input into the HMM estimation
A list of integer arrays with the discrete trajectories mapped to the observation mode used.
Returns whether the MSM is reversible
Returns whether the MSM is sparse
The lag time in steps
Lifetimes of states of the hidden transition matrix
The logger for this class instance
Computes the assignment to metastable sets for observable states
Returns the output probability distributions. Identical to
Computes the memberships of observable states to metastable sets by
Computes the metastable sets of observable states within each
The model estimated by this Estimator
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 active set of states on which all computations and estimations will be done
Ensures that the observable states are indexed and returns the indices
returns the output probability matrix
The stationary distribution on the MSM states
Returns whether the MSM is reversible
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_set
¶ The active set of hidden states on which all hidden state computations are done
-
cktest
(mlags=10, conf=0.95, err_est=False, n_jobs=None, show_progress=True)¶ Conducts a Chapman-Kolmogorow test.
- Parameters
mlags (int or int-array, default=10) – 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, default = 0.95) – confidence interval
err_est (bool, default=False) – 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, default=True) – Show progressbars for calculation?
- Returns
cktest
- Return type
References
This is an adaption of the Chapman-Kolmogorov Test described in detail in [1]_ to Hidden MSMs as described in [2]_.
- 1
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
- 2
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
-
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.
-
discrete_trajectories_full
¶ A list of integer arrays with the original trajectories.
-
discrete_trajectories_lagged
¶ Transformed original trajectories that are used as an input into the HMM estimation
-
discrete_trajectories_obs
¶ A list of integer arrays with the discrete trajectories mapped to the observation mode used. When using observe_active = True, the indexes will be given on the MSM active set. Frames that are not in the observation set will be -1. When observe_active = False, this attribute is identical to discrete_trajectories_full
-
dt_model
¶ Description of the physical time corresponding to the lag.
-
dt_traj
¶
-
dtrajs_full
¶ A list of integer arrays with the original trajectories.
-
dtrajs_lagged
¶ Transformed original trajectories that are used as an input into the HMM estimation
-
dtrajs_obs
¶ A list of integer arrays with the discrete trajectories mapped to the observation mode used. When using observe_active = True, the indexes will be given on the MSM active set. Frames that are not in the observation set will be -1. When observe_active = False, this attribute is identical to discrete_trajectories_full
-
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_left_obs
¶
-
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)
-
eigenvectors_right_obs
¶
-
estimate
(X, **params)¶ Estimates the model given the data X
- Parameters
X (object) – A reference to the data from which the model will be estimated
params (dict) – New estimation parameter values. The parameters must that have been announced in the __init__ method of this estimator. The present settings will overwrite the settings of parameters given in the __init__ method, i.e. the parameter values after this call will be those that have been used for this estimation. Use this option if only one or a few parameters change with respect to the __init__ settings for this run, and if you don’t need to remember the original settings of these changed parameters.
- Returns
estimator – The estimated estimator with the model being available.
- Return type
object
-
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
-
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
-
is_reversible
¶ Returns whether the MSM is reversible
-
is_sparse
¶ Returns whether the MSM is sparse
-
lagtime
¶ The lag time in steps
-
lifetimes
¶ Lifetimes of states of the hidden transition matrix
- Returns
l – state lifetimes in units of the input trajectory time step, defined by \(-\tau / ln \mid p_{ii} \mid, i = 1,...,nstates\), where \(p_{ii}\) are the diagonal entries of the hidden transition matrix.
- Return type
ndarray(nstates)
-
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
¶ Computes the assignment to metastable sets for observable states
Notes
This is only recommended for visualization purposes. You cannot compute any actual quantity of the coarse-grained kinetics without employing the fuzzy memberships!
- Returns
For each observable state, the metastable state it is located in.
- Return type
ndarray((n) ,dtype=int)
-
metastable_distributions
¶ - Returns the output probability distributions. Identical to
- Returns
Pout – output probability matrix from hidden to observable discrete states
- Return type
ndarray (m,n)
See also
-
metastable_memberships
¶ - Computes the memberships of observable states to metastable sets by
Bayesian inversion as described in [1]_.
- Returns
M – A matrix containing the probability or membership of each observable state to be assigned to each metastable or hidden state. The row sums of M are 1.
- Return type
ndarray((n,m))
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)
-
metastable_sets
¶ - Computes the metastable sets of observable states within each
metastable set
Notes
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 observable state indexes contained in it
- Return type
list of int-arrays
-
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
-
msm_init
¶
-
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_obs
¶ Number of states in discrete trajectories
-
observable_set
¶ The active set of states on which all computations and estimations will be done
-
observable_state_indexes
¶ Ensures that the observable states are indexed and returns the indices
-
observation_probabilities
¶ returns the output probability matrix
- Returns
Pout – output probability matrix from hidden to observable discrete states
- Return type
ndarray (m,n)
-
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_observation_probabilities
(nsample)¶ Generates samples according to the current observation probability distribution
- Parameters
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) )
-
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
-
set_model_params
(P=None, pobs=None, pi=None, reversible=None, dt_model='1 step', neig=None)¶ - Parameters
P (ndarray(m,m)) – coarse-grained or hidden transition matrix
pobs (ndarray (m,n)) – observation probability matrix from hidden to observable discrete states
pi (ndarray(m), optional, default=None) – stationary or distribution. Can be optionally given in case if it was already computed, e.g. by the estimator.
reversible (bool, optional, default=None) – whether P is reversible with respect to its stationary distribution. If None (default), will be determined from P
dt_model (str, optional, default='1 step') – time step of the model
neig (int or None) – The number of eigenvalues / eigenvectors to be kept. If set to None, all eigenvalues will be used.
-
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
-
simulate
(N, start=None, stop=None, dt=1)¶ Generates a realization of the Hidden 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 ((N/dt, ) ndarray) – The hidden state trajectory with length N/dt
otraj ((N/dt, ) ndarray) – The observable state discrete trajectory with length N/dt
-
sparse
¶ Returns whether the MSM is sparse
-
stationary_distribution
¶ The stationary distribution on the MSM states
-
stationary_distribution_obs
¶
-
submodel
(states=None, obs=None, mincount_connectivity='1/n', inplace=False)¶ Returns a HMM with restricted state space
- Parameters
states (None, str or int-array) – Hidden states to restrict the model to. In addition to specifying the subset, possible options are: * None : all states - don’t restrict * ‘populous-strong’ : strongly connected subset with maximum counts * ‘populous-weak’ : weakly connected subset with maximum counts * ‘largest-strong’ : strongly connected subset with maximum size * ‘largest-weak’ : weakly connected subset with maximum size
obs (None, str or int-array) – Observed states to restrict the model to. In addition to specifying an array with the state labels to be observed, possible options are: * None : all states - don’t restrict * ‘nonempty’ : all states with at least one observation in the estimator
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. Default value: 1/nstates.
inplace (Bool) – if True, submodel is estimated in-place, overwriting the original estimator and possibly discarding information. Default value: False
- Returns
hmm – The restricted HMM.
- Return type
HMM
-
submodel_disconnect
(mincount_connectivity='1/n')¶ Disconnects sets of hidden states that are barely connected
Runs a connectivity check excluding all transition counts below mincount_connectivity. The transition matrix and stationary distribution will be re-estimated. Note that the resulting transition matrix may have both strongly and weakly connected subsets.
- Parameters
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
hmm – The restricted HMM.
- Return type
HMM
-
submodel_largest
(strong=True, mincount_connectivity='1/n')¶ Returns the largest connected sub-HMM (convenience function)
- Returns
hmm – The restricted HMM.
- Return type
HMM
-
submodel_populous
(strong=True, mincount_connectivity='1/n')¶ Returns the most populous connected sub-HMM (convenience function)
- Returns
hmm – The restricted HMM.
- Return type
HMM
-
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 HMSM 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
The normalized trajectory weights. Given \(N\) trajectories of lengths \(T_1\) to \(T_N\),
returns the corresponding weights
.. math:: – (w_{1,1}, …, w_{1,T_1}), (w_{N,1}, …, w_{N,T_N})
-
transition_matrix
¶ The transition matrix on the active set.
-
transition_matrix_obs
(k=1)¶ Computes the transition matrix between observed states
Transition matrices for longer lag times than the one used to parametrize this HMSM can be obtained by setting the k option. Note that a HMSM is not Markovian, thus we cannot compute transition matrices at longer lag times using the Chapman-Kolmogorow equality. I.e.:
\[P (k \tau) \neq P^k (\tau)\]This function computes the correct transition matrix using the metastable (coarse) transition matrix \(P_c\) as:
\[P (k \tau) = {\Pi}^-1 \chi^{\top} ({\Pi}_c) P_c^k (\tau) \chi\]where \(\chi\) is the output probability matrix, \(\Pi_c\) is a diagonal matrix with the metastable-state (coarse) stationary distribution and \(\Pi\) is a diagonal matrix with the observable-state stationary distribution.
- Parameters
k (int, optional, default=1) – Multiple of the lag time for which the By default (k=1), the transition matrix at the lag time used to construct this HMSM will be returned. If a higher power is given,
-
update_model_params
(**params)¶ Update given model parameter if they are set to specific values
References
[1]_ is an excellent review of estimation algorithms for discrete Hidden Markov Models. This function estimates a discrete HMM on the discrete input states using the Baum-Welch algorithm [2]_. We use a maximum-likelihood Markov state model to initialize the HMM estimation as described in 3.
- 1
L. R. Rabiner: A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition. Proc. IEEE 77, 257-286 (1989)
- 2
L. Baum, T. Petrie, G. Soules and N. Weiss N: A maximization technique occurring in the statistical analysis of probabilistic functions of Markov chains. Ann. Math. Statist. 41, 164-171 (1970)
- 3
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)