deeptime.markov.tools.estimation.transition_matrix

deeptime.markov.tools.estimation.transition_matrix(C, reversible=False, mu=None, method='auto', maxiter: int = 1000000, maxerr: float = 1e-08, rev_pisym: bool = False, return_statdist: bool = False, warn_not_converged: bool = True)

Estimate the transition matrix from the given countmatrix. [1] [2] [3]

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.

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

• maxerr (float, optional, default=1e-8) – 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, optional, 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, optional, default=False) – Optional parameter with reversible = True. If set to true, the stationary distribution is also returned

• return_conv (bool, optional, 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, optional, default=True) – Prints a warning if not converged.

Returns:

result – The MLE transition matrix and optionally the stationary distribution if return_statist=True. The transition matrix has the same data type (dense or sparse) as the input matrix C.

Return type:

Union[array_like, Tuple[array_like, np.ndarray]]

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]

Jan-Hendrik Prinz, Hao Wu, Marco Sarich, Bettina Keller, Martin Senne, Martin Held, John D Chodera, Christof Schütte, and Frank Noé. Markov models of molecular kinetics: generation and validation. The Journal of chemical physics, 134(17):174105, 2011.

[2]

Gregory R Bowman, Kyle A Beauchamp, George Boxer, and Vijay S Pande. Progress and challenges in the automated construction of markov state models for full protein systems. The Journal of chemical physics, 131(12):124101, 2009.

[3]

Benjamin Trendelkamp-Schroer, Hao Wu, Fabian Paul, and Frank Noé. Estimation and uncertainty of reversible markov models. The Journal of chemical physics, 143(17):11B601_1, 2015.

Examples

>>> import numpy as np
>>> from deeptime.markov.tools.estimation import transition_matrix

>>> C = np.array([[10, 1, 1], [2, 0, 3], [0, 1, 4]])

Non-reversible estimate

>>> T_nrev = transition_matrix(C)
>>> print(np.array_str(T_nrev, precision=3))
[[0.833 0.083 0.083]
 [0.4   0.    0.6  ]
 [0.    0.2   0.8  ]]

Reversible estimate

>>> T_rev = transition_matrix(C, reversible=True)
>>> print(np.array_str(T_rev, precision=3))
[[0.833 0.104 0.063]
 [0.351 0.    0.649]
 [0.049 0.151 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)
>>> print(np.array_str(T_mu, precision=3))
[[0.948 0.006 0.046]
 [0.429 0.    0.571]
 [0.111 0.02  0.869]]