msmtools.estimation.tmatrix

msmtools.estimation.tmatrix(C, reversible=False, mu=None, method='auto', **kwargs)

Estimate the transition matrix from the given countmatrix.

Parameters:
  • C (numpy ndarray or scipy.sparse matrix) – Count matrix
  • reversible (bool (optional)) – If True restrict the ensemble of transition matrices to those having a detailed balance symmetry otherwise the likelihood optimization is carried out over the whole space of stochastic matrices.
  • mu (array_like) – The stationary distribution of the MLE transition matrix.
  • method (str) – Select which implementation to use for the estimation. One of ‘auto’, ‘dense’ and ‘sparse’, optional, default=’auto’. ‘dense’ always selects the dense implementation, ‘sparse’ always selects the sparse one. ‘auto’ selects the most efficient implementation according to the sparsity structure of the matrix: if the occupation of the C matrix is less then one third, select sparse. Else select dense. The type of the T matrix returned always matches the type of the C matrix, irrespective of the method that was used to compute it.
  • **kwargs (Optional algorithm-specific parameters. See below for special cases) –
  • Xinit ((M, M) ndarray) – Optional parameter with reversible = True. initial value for the matrix of absolute transition probabilities. Unless set otherwise, will use X = diag(pi) t, where T is a nonreversible transition matrix estimated from C, i.e. T_ij = c_ij / sum_k c_ik, and pi is its stationary distribution.
  • maxiter (1000000 : int) – Optional parameter with reversible = True. maximum number of iterations before the method exits
  • maxerr (1e-8 : float) – Optional parameter with reversible = True. convergence tolerance for transition matrix estimation. This specifies the maximum change of the Euclidean norm of relative stationary probabilities (\(x_i = \sum_k x_{ik}\)). The relative stationary probability changes \(e_i = (x_i^{(1)} - x_i^{(2)})/(x_i^{(1)} + x_i^{(2)})\) are used in order to track changes in small probabilities. The Euclidean norm of the change vector, \(|e_i|_2\), is compared to maxerr.
  • rev_pisym (bool, default=False) – Fast computation of reversible transition matrix by normalizing \(x_{ij} = \pi_i p_{ij} + \pi_j p_{ji}\). \(p_{ij}\) is the direct (nonreversible) estimate and \(pi_i\) is its stationary distribution. This estimator is asympotically unbiased but not maximum likelihood.
  • return_statdist (bool, default=False) – Optional parameter with reversible = True. If set to true, the stationary distribution is also returned
  • return_conv (bool, default=False) – Optional parameter with reversible = True. If set to true, the likelihood history and the pi_change history is returned.
  • warn_not_converged (bool, default=True) – Prints a warning if not converged.
  • sparse_newton (bool, default=False) – If True, use the experimental primal-dual interior-point solver for sparse input/computation method.
Returns:

  • P ((M, M) ndarray or scipy.sparse matrix) – The MLE transition matrix. P has the same data type (dense or sparse) as the input matrix C.
  • The reversible estimator returns by default only P, but may also return
  • (P,pi) or (P,lhist,pi_changes) or (P,pi,lhist,pi_changes) depending on the return settings
  • P (ndarray (n,n)) – transition matrix. This is the only return for return_statdist = False, return_conv = False
  • (pi) (ndarray (n)) – stationary distribution. Only returned if return_statdist = True
  • (lhist) (ndarray (k)) – likelihood history. Has the length of the number of iterations needed. Only returned if return_conv = True
  • (pi_changes) (ndarray (k)) – history of likelihood history. Has the length of the number of iterations needed. Only returned if return_conv = True

Notes

The transition matrix is a maximum likelihood estimate (MLE) of the probability distribution of transition matrices with parameters given by the count matrix.

References

[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]Bowman, G R, K A Beauchamp, G Boxer and V S Pande. 2009. Progress and challenges in the automated construction of Markov state models for full protein systems. J. Chem. Phys. 131: 124101
[3]Trendelkamp-Schroer, B, H Wu, F Paul and F. Noe. 2015 Estimation and uncertainty of reversible Markov models. J. Chem. Phys. 143: 174101

Examples

>>> import numpy as np
>>> from msmtools.estimation import transition_matrix
>>> C = np.array([[10, 1, 1], [2, 0, 3], [0, 1, 4]])

Non-reversible estimate

>>> T_nrev = transition_matrix(C)
>>> T_nrev
array([[ 0.83333333,  0.08333333,  0.08333333],
       [ 0.4       ,  0.        ,  0.6       ],
       [ 0.        ,  0.2       ,  0.8       ]])

Reversible estimate

>>> T_rev = transition_matrix(C, reversible=True)
>>> T_rev
array([[ 0.83333333,  0.10385551,  0.06281115],
       [ 0.35074677,  0.        ,  0.64925323],
       [ 0.04925323,  0.15074677,  0.8       ]])

Reversible estimate with given stationary vector

>>> mu = np.array([0.7, 0.01, 0.29])
>>> T_mu = transition_matrix(C, reversible=True, mu=mu)
>>> T_mu
array([[ 0.94771371,  0.00612645,  0.04615984],
       [ 0.42885157,  0.        ,  0.57114843],
       [ 0.11142031,  0.01969477,  0.86888491]])