Continuous Time Markov Chain

Questions
1

I’m somewhat new to PyMc, but I can’t seem to find anything to help out on this.

For relevant background on CTMC, see the answer to: https://math.stackexchange.com/questions/876789/continuous-time-markov-chains-is-this-step-by-step-example-correct

A toy example would be to consider I have 2 states. I have a generator matrix Q which is parameterised as:

Q = \begin{bmatrix} - \lambda_{1} & \lambda_{1} \\ \lambda_{2} & -\lambda_{2} \\ \end{bmatrix}

My transition matrix for a period of time, t, is given by exp(Qt), where exp is the matrix exponential.

The above link explains how to fit MLE to the generator matrix, from data of transition times. If I wanted a prior over \lambda_{1}, \ \lambda_{2}, how would I do this in PyMC3? Obviously this is a toy example, and the parameters may be generally more complex.

I appreciate any help you guys can offer.

2

Defining a matrix Q by a set of random variables is simple enough, using

r_12 = pm.HalfNormal('r12', 5.)
r_13 = pm.HalfNormal('r13', 5.)
r_11 = pm.Deterministic('r11', -r_12 - r_13)
(...)
Q = tt.stacklists([[r_11, r_12, r_13], [r_21, r_22, r_23], [r_31, r_32, r_33]])

and P = tt.slinalg.Expm(Q * s) can give a transition matrix for a particular interval.

From here, it seems like there are two approaches. The first is to discretize the process, letting s be the desired step size, and therefore P be a fixed transition matrix, and treat this as a parameter to a discrete probability distribution like

class DiscreteMarkov(pm.Categorical):  # maybe (pm.Discrete) ?
    def __init__(self, P, init_P):
        super(pm.Discrete, self).__init__(*args, **kwargs)
        self.P = P
        self.init_P = init_P

   def logp(self, x):
        # assume x is encoded as state index [0,0,0,1,1,2,1,3,3,0,0,0, ...]
        trans = self.P[x[:-1]] # no transition from the last state
        lik_chain = tt.sum(pm.Categorical.dist(trans).logp(x[1:]))  # after the initial state
        return lik_chain + pm.Categorical.dist(self.init_P).logp(x[0]) # add initial state

An alternative is to create your own Distribution over time-state matrices. It would look much like the above but with something like

def logp(self, X):
    time = X[:,0]
    state = X[:,1]
    s = time[1:] - time[:-1]
    lik = pm.Categorical.dist(self.init_P).logp(state[0])
    for i in xrange(len(s)):
        P = tt.slinalg.Expm(self.Q * s[i])
        lik += pm.Categorical.dist(P[state[i]]).logp(state[i+1])
    return lik

No guarantee that this will work. It’s a valid likelihood over a time-state matrix; but i’m not sure if the pymc3 machinery is amazing enough to draw from it without knowing that X[:,0] is non-negative ascending, and X[:,1] is bounded non-negative integers. Someone more familiar with the API could comment…

Edit:

It actually looks like you could do this with DensityDist by the above logp, and defining your own random sampler. Since the CTMC has exponential transition times (conditioned on state), the generator would sample the initial state; then given a state, sample the transition time, compute the transition matrix P, sample the ending state, and recurse until the time-window has been exceeded.

3 Likes