discreteMarkovChain 0.22

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 one-dimensional 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[state-1] = self.downrate
        elif state == self.m:
            rates[state+1] = self.uprate
        elif state == self.M:
            rates[state-1] = self.downrate
        return rates

Now we initialize the random walk with some values for m and M and calculate the steady-state 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 multi-dimensional 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.