elynx-seq-0.0.1: Handle molecular sequences

Copyright(c) Dominik Schrempf 2017
Portabilitynon-portable (not tested)
Safe HaskellNone



Some helper functions that come handy when working with rate matrices of continuous-time discrete-state Markov processes.

  • Changelog

To be imported qualified.



type RateMatrix = Matrix R Source #

A rate matrix is just a real matrix.

type ExchangeabilityMatrix = Matrix R Source #

A matrix of exchangeabilities, we have q = e * pi, where q is a rate matrix, e is the exchangeability matrix and pi is the diagonal matrix containing the stationary frequency distribution.

type StationaryDistribution = Vector R Source #

Stationary distribution of a rate matrix.

totalRate :: StationaryDistribution -> RateMatrix -> Double Source #

Get average number of substitutions per unit time.

normalize :: RateMatrix -> RateMatrix Source #

Normalizes a Markov process generator such that one event happens per unit time. Calculates stationary distribution from rate matrix.

normalizeWith :: StationaryDistribution -> RateMatrix -> RateMatrix Source #

Normalizes a Markov process generator such that one event happens per unit time. Stationary distribution has to be given.

setDiagonal :: RateMatrix -> RateMatrix Source #

Set the diagonal entries of a matrix such that the rows sum to 0.

toExchangeabilityMatrix :: RateMatrix -> StationaryDistribution -> ExchangeabilityMatrix Source #

Extract the exchangeability matrix from a rate matrix.

fromExchangeabilityMatrix :: ExchangeabilityMatrix -> StationaryDistribution -> RateMatrix Source #

Convert exchangeability matrix to rate matrix.

getStationaryDistribution :: RateMatrix -> StationaryDistribution Source #

Get stationary distribution from RateMatrix. Involves eigendecomposition. If the given matrix does not satisfy the required properties of transition rate matrices and no eigenvector with an eigenvalue nearly equal to 0 is found, an error is thrown. Is there an easier way to calculate the stationary distribution or a better way to handle errors (of course I could use the Maybe monad, but then the error report is just delayed to the calling function)?