Solve Markov chains with a discrete state space.
Project description
While for statistical and scientific programming languages such as R various packages are available for analyzing Markov chains, equivalent packages in Python are rather scarce. This discreteMarkovChain package for Python addresses the problem of obtaining the steady state distribution of a Markov chain, also known as the stationary distribution, limiting distribution or invariant measure. The package is for Markov chains with discrete and finite state spaces, which are most commonly encountered in practical applications.
This package is based on numpy and scipy for efficient computations and limited use of resources. Markov chains with several million states can be solved. The package introduces the markovChain class which has the following features.
 States can be either integers or vectors of integers.
 Steady state distributions can be calculated for continous time Markov chains (CTMC) as well as discrete time Markov chains (DTMC).
 The user can either manually specify the probability/rate matrix of the Markov chain, or let the program do this automatically using an indirect or direct method.
 The indirect method requires the user to specify an initial state and transition function (giving for each state the reachable states and their probabilities/rates).
 By repeatedly calling the transition function on unvisited states, the state space and the probability matrix are built up automatically.
 This makes it easy to implement your own Markov chains!
 The direct method requires the user to specify a transition function and a function that gives the complete state space.
 While the implementation is typically more complex, this may have some computational advantage over the indirect method for large state spaces with vector states.
 The steady state distribution can be calculated by a method of choice:
 The power method,
 Solving a system of linear equations,
 Determing the first left eigenvector,
 Searching in Krylov subspace.
 Checks are included to see whether all states in the Markov chain are connected.
 Memory consumption is reduced by using sparse matrices.
When the user calls a certain solution method, the markovChain object gets the attribute pi which specifies the steady state probability of each state. When the user uses the direct or indirect method, the object gets the mapping attribute which is a dictionary that links each index of pi with a corresponding state. Using the mapping and pi, it becomes simple to calculate performance measures for your Markov chain, such as the average cost per time unit or the number of blocked customers in a queue with blocking.
Installation
The package can be installed with the command
pip install discreteMarkovChain
or by downloading the source distribution and installing manually with
python setup.py install
Examples
The markovChain class can be used to initialize your own Markov chains. We import it by using
from discreteMarkovChain import markovChain
First, lets consider a simple Markov chain with two states, where we already know the probability matrix P.
P = np.array([[0.5,0.5],[0.6,0.4]]) mc = markovChain(P) mc.computePi('linear') #We can also use 'power', 'krylov' or 'eigen' print(mc.pi)
We get the following steady state probabilities:
[ 0.54545455 0.45454545]
Now we show an example of a onedimensional random walk in continuous time between integers m and M. We move up and down with rates 1. We will use the indirect method to determine the rate matrix for us automatically. The indirect method is rather flexible, and allows the transition function to return a dictionary with reachable states and rates. We first introduce our randomWalk class.
class randomWalk(markovChain): #A random walk where we move up and down with rate 1.0 in each state between bounds m and M. #For the transition function to work well, we define some class variables in the __init__ function. def __init__(self,m,M): super(randomWalk, self).__init__() #always use this as first line when creating your own __init__ self.initialState = m self.m = m self.M = M self.uprate = 1.0 self.downrate = 1.0 def transition(self,state): #Specify the reachable states from state and their rates. #A dictionary is extremely easy here! rates = {} if self.m < state < self.M: rates[state+1] = self.uprate rates[state1] = self.downrate elif state == self.m: rates[state+1] = self.uprate elif state == self.M: rates[state1] = self.downrate return rates
Now we initialize the random walk with some values for m and M and calculate the steadystate vector pi.
mc = randomWalk(0,5) mc.computePi() mc.printPi()
The stationary probabilities are given below.
0 0.166666666667 1 0.166666666667 2 0.166666666667 3 0.166666666667 4 0.166666666667 5 0.166666666667
Not unexpectedly, they are the same for each state. We can repeat this for a multidimensional random walk. Now we use the direct method. Here, we need to use a transition function returning numpy arrays and we need to define a function that calculates the state space.
from discreteMarkovChain import partition class randomWalkNumpy(markovChain): #Now we do the same thing with a transition function that returns a 2d numpy array. #We also specify the statespace function so we can use the direct method. #This one is defined immediately for general n. def __init__(self,m,M,n,direct=True): super(randomWalkNumpy, self).__init__(direct=direct) self.initialState = m*np.ones(n,dtype=int) self.n = n self.m = m self.M = M self.uprate = 1.0 self.downrate = 1.0 #It is useful to define the variable 'events' for the the transition function. #The possible events are 'move up' or 'move down' in one of the random walks. #The rates of these events are given in 'eventRates'. self.events = np.vstack((np.eye(n,dtype=int),np.eye(n,dtype=int))) self.eventRates = np.array([self.uprate]*n+[self.downrate]*n) def transition(self,state): #First check for the current state which of the 'move up' and 'move down' events are possible. up = state < self.M down = state > self.m possibleEvents = np.concatenate((up,down)) #Combine into one boolean array. #The possible states after the transition follow by adding the possible 'move up'/'move down' events to the current state. newstates = state+self.events[possibleEvents] rates = self.eventRates[possibleEvents] return newstates,rates def statespace(self): #Each random walk can be in a state between m and M. #The function partition() gives all partitions of integers between min_range and max_range. min_range = [self.m]*self.n max_range = [self.M]*self.n return partition(min_range,max_range)
Now we initialize n=2 random walks between m=0 and M=2 and print the stationary distribution.
mc = randomWalkNumpy(0,2,n=2) mc.computePi('linear') mc.printPi() [0 0] 0.111111111111 [1 0] 0.111111111111 [2 0] 0.111111111111 [0 1] 0.111111111111 [1 1] 0.111111111111 [2 1] 0.111111111111 [0 2] 0.111111111111 [1 2] 0.111111111111 [2 2] 0.111111111111
We could also solve much larger models. The example below has random walks in 5 dimensions with 100.000 states. For these larger models, it is often better to use the power method. The linear algebra solver may run into memory problems.
mc = randomWalkNumpy(0,9,n=5) mc.computePi('power')
On a dual core computer from 2006, the rate matrix and pi can be calculated within 10 seconds.
Changes in v0.22
 Added documentation for the markovChain class and all its methods, including examples.
 Added the function partition that can be used to determine the state space when states are consists of all integers between ranges. The optional parameter max_sum can be specified if the state vectors should sum up to less than max_sum (useful in some queueing and inventory applications).
 Fixed an error when calling krylovMethod(), linearMethod() and eigenMethod() on Markov chains with one state.
 Included a workaround for an error when calling eigenMethod() on Markov chains with two states.
Project details
Release history Release notifications  RSS feed
Download files
Download the file for your platform. If you're not sure which to choose, learn more about installing packages.
Filename, size  File type  Python version  Upload date  Hashes 

Filename, size discreteMarkovChain0.22py2.py3noneany.whl (20.0 kB)  File type Wheel  Python version 2.7  Upload date  Hashes View 
Filename, size discreteMarkovChain0.22.tar.gz (14.9 kB)  File type Source  Python version None  Upload date  Hashes View 
Filename, size discreteMarkovChain0.22.zip (24.2 kB)  File type Source  Python version None  Upload date  Hashes View 
Hashes for discreteMarkovChain0.22py2.py3noneany.whl
Algorithm  Hash digest  

SHA256  6893a74275535baaeabda074a5dec030e916e39b1b5c588580617b218a396302 

MD5  f23fde03fa4dd7eb73b26aadfd01b8d9 

BLAKE2256  7785d1d2e71296f9cd4f5d8d1e2c198dcb44bbc056722351ca24e70b3c45b090 