Gillespie Decomposition of a Markov Chain

Goal

We will discuss a way to decompose a continuous time Markov chain (CTMC) into a product of two probability distributions:

1. The probability distribution for the waiting time until the next transition.
2. The probability distribution for which transition will occur when the transition takes place.

Gillespie used this decomposition to formulate his algorithm for simulating reaction kinetics. The concepts for this decomposition exist in most treatments of CTMCs, but are rarely put together so explicitly.We will quickly review the assumptions and properties of CTMCs. Then we will demonstrate how to achieve the decomposition just described.

Continuous Time Markov Chain

A CTMC is a mathematical model of systems that three properties:

1. The system resides in a finite number of discrete states
2. The system transitions from state to state in a particular way such that the probability of a transition in a small time interval $\Delta t$ is proportional to $\lambda \Delta t$ for some constant $\lambda$
3. Once the system is in a given a given state it retains no memory of past states or transitions (or any other aspect of its history). In other words, it is "memoryless"

CTMCs are defined by a set of transition probabilities that form a matrix $p_{i,j}$ where $j$ represents the current state of the system, and $i$ represents the new state after a transition. In a CTMC, the $p_{i,j}$ represent infinitesimal rates. Because a state must transition to somewhere the columns of $p$ must sum to 1. Thus, the rows can be considered conditional probabilities of the form:

$p_{i,j} = P[\text{in state } i \text{ at } t+\Deltat ,| \text{ in state } j \text{ at } t]$

In fact, the $p_{i,j}$ are the state-dependent $\lambda$s referred to in item (2) above. Therefore, the are all $o(\Delta t)$.

Connection to Poisson Process

Because a CTMC is based on the same fundamental assumptions as the Poisson process (which is no coincidence), they share some important features. Actually, a Poisson process is a special case of a CTMC in which the states are integer counting events.

Because of this relationship, the waiting time for every transition $i \to j$ is also an exponential distribution, as we will now derive. The Poisson process has a single parameter $\lambda$. A CTMC has an entire matrix of parameters $p_{i,j}$. We must figure out how to incorporate this into our derivation.

Just as with the Poisson process, we begin by deriving the probability of NO transition in a time interval $\tau$ (just as with the Poisson process, because of memorylessness, we are only concerned with intervals). If we are currently in state $j$, then the applicable rate (equivalent to $\lambda$ for the Poisson process) is $p_{j,j}$. Since we have already derived this result, we can immediately say that

$P[\text{still in } j \text{ at } t+\tau \,|\, \text{was in } j \text{ at } t] = e^{-p_{j,j}\,t}$

Math in compose mode

Test of $\sqrt{\sigma}$ inline math in compose mode.

Environment mode:
\begin{equation}
\frac{1}{\sqrt{sigma}}
\end{equation}

Continuous Time Markov Chains 1.

Note: This contains some mistakes and needs to be rewritten. In general it is correct.

Goal

We will define and show some basic properties of Markov chains (MCs), and then discuss continuous time MCs.

Markov Chains

MCs are descriptions of discrete states linked by transition probabilities. Suppose a system can be in three states: $A$, $B$, and $C$. As time goes on, this system transitions from one state to another. For instance, $A \to C$ followed by $C \to A$ followed by $A \to B$. By observing this system we can approximately measure the probabilities or rates at which the system transitions occur.

If the system evolves deterministically, then it will spend exactly a specific amount of time in one state before transitioning to another state. The inverse of this time is the rate for the transition.

If the system evolves randomly, then it may remain in a state for a random amount of time before transitioning. In this case, we have transition probabilities rather than rates. In this respect, there are still important differences between discrete time and continuous time MCs.

We will continue to discuss networks and discrete time Markov Chains (DTMCs) to some degree in order to highlight the particular properties of continuous time Markov Chains (CTMCs). However, this note will focus on CTMCs.

Memory

For many important applications, systems that can be modeled by MCs evolve in a memoryless fashion. Actually, few MCs that are not memoryless can easily be analyzed in closed form, which is probably why we like model systems amenable to memoryless MCs.

Discrete Time Markov Chains

In a deterministic system there can only be one possible transition from each state. Otherwise, how would it choose which way to go? The next step in randomness is to have a number of possible transitions from each state occurring randomly with measurable probabilities. For example, the transition $A \to B$ might occur with a probability $p_{A \to B} = 0.2$ and the transition $A \to C$ occur with probability $p_{A \to C} = 0.7$. Thus, in a CTMC, we have no sense of physical time advancing. Instead, we index successive states of the system by some integer, say $i$ or $n$ or $\alpha$. Notice that the probabilities in my example add to $1$. When transitioning out of the state $A$, the system has to go somewhere! We will see an interesting modification of this when we consider events in continuous time.

Representations of Markov Chains

Two convenient representations of MCs are employed. The first is graphical. Essentially, a dependency graph. Graph elements representing individual states are connected with arrows representing possible transitions. I'll make one to show as soon as I install and learn to use canviz.

The second representation is in matrix form and lends itself to mathematical analysis. Both rows and columns represent the different states of the system, thus the matrix is square. Each entry describes the transition from the row state to the column state. This description will be different depending on whether the the system evolves deterministically, randomly, and/or in discrete or continuous time (more on that later).

Here is a matrix that represents all the transitions in a hypothetical CTMC:

\begin{bmatrix} 0.1 & 0.7 & 0.2 \\ 0.8 & 0.1 $ 0.1 \\ 0.4 & 0.3 & 0.3 \\ \end{bmatrix}

If our example were deterministic network, then the matrix could instead just have ones and zeros (for discrete time) or it could have rates which would represent approximately how long you need to wait or the system to change state. Also in the case of a deterministic network, or a random but discrete time MC, there would be no numbers on the diagonal, since we would again only be concerned with changes in state, and no strict passage of time. However, in the latter case, the rows would add to 1. This is because the numbers in the rows would represent probabilities for transitioning to the state in the respective column. Since the systems transitions to somewhere with probability 1, the rows must sum to 1.

A CTMC has a very important difference, which is reflected in the fact that it does have entries on the diagonal, as in our example. A CTMC represents a system in which time is advancing continuously. The matrix entries represent infinitesimal transition rates for a very tiny time step $\delta t$. It is understood that the system will eventually go to a new state.

However, in any particular infinitesimal interval the system may NOT transition. We represent this by a probability or rate for the system to stay in its current state. In the diagrammatic representation this is expressed as an arrow that loops back on the same state. In the matrix representation this is expressed by transition rates on the diagonal, since these represent the probabilities that the system will stay in its current state during the next $\Delta t$.

Examples

Here is a short list of MC applications that is both incomplete insufficient in detail. It does give a sense of the broad range of applications as well as the different ways one can think of states and transitions.

As mentioned above, CTMCs can model systems that are or can be thought of as being approximately memoryless, which again means that the system does not remember what states it has been in leading up to being in the current state -- systems that live entirely in the moment.

Queues

One of the important practical applications of MCs is in modeling queuing systems. A classic example is machines which may break down randomly at any time, and one or more technicians who service the machines. Perhaps the servicing time also has some probability distribution. As will be discussed in another note, if the probability of a machine breaking down at any moment is fairly small, then the machine breakdowns are essentially memoryless. In this example, the "state" of the system is the list of which machines are active and which are in repair. This classic example also has numerous equivalent variants: people in line to get service at some counter, requests to an internet server, etc.

Population and Counting

Another classic application is in various "counting processes". In a counting process, the various states of the system are the numbers of something, such as the number of molecules of a certain chemical in a reaction, the number of individuals in a population, etc. An important class of these is called "birth and death processes", which have application in biology including population dynamics of species and the spread of epidemics.

Real and Hypothetical Networks

Google uses a Markov type model to rank pages. An world wide web user remains on a certain page for a certain amount of time before going to another. Thus, there is a probabilistic rate associated with page links. The world-wide-web is somewhere between a physical network of interconnected physical entities and a virtual network of hypothetical entities or states. The internet, which is the physical network of servers, is a physical network. Markov models can also be used to model network traffic. For instance, the path of a single packet through the internet can be a Markov chain. The packet travels to a server and remains there for a certain amount of time before being routed to another server.

Climate

Markov chains are often used in modeling climate. By discretizing the values of a climate variable such as temperature or pressuer, these values can be thought of as the states of a Markov chain. Because climate has natural periodicities, one can use DTMCs for modeling successive days, for instance. Or, for events that occur intermittently, such as rain or extreme winds, one can use CTMCs.

Summary and Conclusion

CTMCs are a useful tool for modeling many real world processes that can be described as a system that will reside for random periods of time in one of many states before transitioning to another, randomly chosen state. The additional requirement is that the time between transitions as well as the choice of transition are completely random and unaffected by the history of the system, including even the immediately preceding transition.

CTMCs also have many useful and interesting properties, some of which we will discuss in other notes. These include the Poisson Process, which governs the probability distribution of waiting times between transitions between states.

jsMath problems in blogger

Inline math

This is an example of $\frac{1}{\sqrt{2 \pi \sigma^2}}$ inline math

gives

This is an example of $\frac{1}{\sqrt{2 \pi \sigma^2}}$ inline math

Equation Environment

\begin{equation}
p(x) = \frac{1}{\sqrt{2 \pi \sigma^2}} e^-x^2/{2 \sigma^2}
\end{equation}

gives

\begin{equation} p(x) = \frac{1}{\sqrt{2 \pi \sigma^2}} e^-x^2/{2 \sigma^2} \end{equation}

Shorthand Equation Environment

$$p(x) = \frac{1}{\sqrt{2 \pi \sigma^2}} e^-x^2/{2 \sigma^2}$$

gives

$$p(x) = \frac{1}{\sqrt{2 \pi \sigma^2}} e^-x^2/{2 \sigma^2}$$

<myserver>/jsMath/easy/load.js

Poisson Processes

Goal

Probability distribution, at time $t$ for $N$ events in a subsequent time interval $[t, t+\tau]$: $P[N,\tau,t]$. As we shall see, as a practical matter, we primarily focus on calculating $P[1, \tau, t]$.

Assumptions

1. For a small time interval $h$, as $h \to 0$, $P[1,h,t] \to \lambda h$ for all $t$ and $P[N>1,h,t] \to 0$. Therefore $P[0,h,t] \to 1 - \lambda h$.
2. $P$ is independent of $t$: $P[N, \tau, t] = P[N, \tau]$ for all $t$. In other words, it is stationary.
3. The probabilities of events in disjoint intervals are independent. If $t_2-t_1>\tau_1$ then,
   $P[N_2,\tau_2,t_2 ; N_1,\tau_1,t_1] = P[N_2,\tau_2,t_2] \cdot P[N_1,\tau_1,t_1]$
   Alternatively, consider four successive time points $t_4 \geq t_3 \geq t_2 \geq t_1$. Then:
   $P[N_2,t_4-t_3,t_3 ; N_1,t_2-t_1,t_1] = P[N_2,t_4-t_3,t_3] \cdot P[N_1,t_2-t_1,t_1]$
4. $P$ is memoryless. We might feel as though having already waited a long time for an event to occur that the next event should happen sooner. However, since $P$ is memoryless, it does not remember how long it took for the last event to occur.
   Note that this is not the same as being stationary.

Derivation of $P[0,\tau]$

We wamt to derive an expression for $P[N,\tau]$ using the assumed properties listed above. To do so, we start by looking at probabilities for 0 events in time intervals. As an aid in understanding, think about the following situation: you start standing at a bus stop at a certain time $t_1$. You wait until a later time $t_2$ and observe that no bus has arrived. For some reason, you stick around until an even later time $t_3$ and notice that a bus has still not arrived. This is strange, because you know that the buses arrive in a way that obeys to assumptions of the Poisson process. Mathematically, consider three successive times $t_3>t_2>t_1$, and the joint probability
$P[0,t_3-t_1 ; 0, t_2 - t_1]$
By the independence of events in disjoint intervals,
$P[0,t_3-t_1 ; 0, t_2 - t_1] = P[0,t_3-t_2] \cdot P[0,t_2-t_1]$

Alternatively, due to the memoryless assumption, we have

$P[0,t_3-t_1 ; 0, t_2 - t_1] = P[0,t_3-t_1]$
This may seem counterintuitive, but it can be understood if we think abou the fact that we have arbitrarily divide the interval $t_3-t_1$ into two by inserting $t_2$. Doing so should not alter the probability! Therefore, we have the above expression.
Combining the two above, we have
$P[0,t_3-t_1] = P[0,t_3-t_2] \cdot P[0,t_2-t_1]$
This expression in hand, we now derive the analytical form for $P[0,\tau]$. To do so, we differentiate the equation by $t_2$. Along the way, we adopt some abbreviated notation that should be clear enough
$0 = p_0'(t_3-t_2) \cdot p_0(t_2-t_1) + p_0'(t_3-t_2) \cdot p_0(t_2-t_1)$

From this we conclude that

$\frac{1}{p_0(t_3-t_2)}\frac{dp_0(t_3-t_2)}{d(t_2)} = \frac{1}{p_0(t_2-t_1)}\frac{dp_0(t_2-t_1)}{d(t_2)}$

Using a standard technique in separation of variables, we note that the two sides of the equation are respectively functions of only $p_0(t_3-t_2)$ and $p_0(t_2-t_1)$. Therefore each side is separately equal to a constant, $\lambda$. Writing $t_2-t_1$ as $\tau$,

$\frac{1}{p_0(\tau)}\frac{dp_0(\tau)}{d\tau} = \lambda$

which has the unique solution

$p_0(\tau) = e^{-\lambda \tau + C}$

Since $p_0(0) = 0$,

$p_0(\tau) = e^{-\lambda \tau + C}$

Derivation of $P[N,\tau]$

From $p_0(\tau)$ we can proceed to calculate an analytical expression for $p_n(\tau)$ (we will use this notation from now on).

There are many approaches to deriving an analytical expression for $p_n(\tau)$. One way is to solve the "master equation", or "birth and death" equation by one of many methods. One can use generating functions, induction, or other approaches. I found a much more satisfying approach online notes.

Let $p_n(\tau)$ denote the probability of n events in an interval $\tau$. We have just derived that

$p_0(\tau) = e^{-\lambda \tau}$

Next, we derive $p_1(\tau)$ from $p_0(\tau)$. Consider a time interval $\tau=t_3-t_1$. Let the single event occur at a time $t_2$. We don't need to be explicit about when the event occurs. Therefore we can place the event, which takes place in a vanishingly small time $\delta \tau$ anywhere in the interval. Based on the independence of disjoint intervals, we can immediately write

$p_1(t_3-t_1 | \text{event at } t_2) = p_0(t_2-t_1) \cdot p_1(\delta t) \cdot p_0(t-3-t_2)$

To remove the conditioning, we integrate over $t_2$ from $0 \ldots \tau$.

$p_1(t_3-t_1) = \int_0^\tau p_1(t_3-t_1 | \text{event at } t_2) dt_2$

$=\int_0^\tau e^_{-\lambda t_2} \cdot \lambda \cdot e^{-\lambda(\tau-t_2)} dt_2$

$=\lambda \cdot e^_{-\lambda \tau} \cdot \int_0^\tau dt_2$

$=\lambda \tau e^_{-\lambda \tau}$

We can now extend this process for $N$ events by further subdividing the interval. Let $\tau = t_N-t_0$ and $t_i$ be the time of the $i^{th}$ event.

$p_N(\tau | t_1 \ldots t_N) = p_0(t_1-t_0)\Pi_{n \in 2\ldots N} \lambda \cdot p_0(t_n-t_{n-1})$

$=\lambda^N e^{-\lambda \tau}$

Next we remove the conditioning on the $t_1 \ldots t_N$.

$p_N(\tau) = \lambda^N e^{-\lambda \tau} int_{t_0}^{t_2}dt'_1 \int_{t_1}^{t_3}dt'_2 \cdots \int_{t_{n-1}}^{t_{n+1}}dt'_n \cdots \int_{t_{N-2}}^{\tau}dt'_{n-1}$

This integral is going to be nasty do evaluate. Therefore, we use a trick. Although I numbered the times of the events, $t_n$, consecutively, I will consider the indices to be arbitrary markers. This means that the $t_n$ can actually happen in any order. In order to keep the calculation correct, I will now also need to divide by the $N!$ different orders in which the events can happen. Now all of the $t_n$ can range from $t_0$ to $\tau$. This simplifies the integral so that we have

$p_N(\tau) = \frac{1}{N!}\lambda^N e^{-\lambda \tau} \times \ldots$

$int_{t_0}^{t_2}dt'_1 \int_{t_0}^{\tau}dt'_2 \cdots \int_{t_0}^{\tau}dt'_n \cdots \int_{t_0}^{\tau}dt'_{n-1}

which evaluates to

$ = \frac{1}{N!} \lambda^N \tau^N e^{-\lambda \tau} = \frac{(\lambda \tau)^N}{N!}e^{-\lambda \tau}$

Conclusion

This concludes a derivation of the Poisson process. From the three assumptions stated initially, we derived the probability distribution for $N$ events in a time interval $\tau$:

$p_N(\tau) = frac{(\lambda \tau)^N}{N!}e^{-\lambda \tau}$

The Poisson process has many useful applications and also many interesting and useful properties, some of which we will discuss in later posts.

Appendix A -- quick and dirty $p_0(\tau)$

I found a quicker approach to $p_0(\tau)$ here. Actually, I think it is the more common one.

$p_0(\tau + \Delta \tau) = P[\text{no event in } \tau ; \text{ no event in } \Delta \tau]$

By independence,

$p_0(\tau + \Delta \tau) = P[\text{no event in } \tau \cdot P[\text{ no event in } \Delta \tau]$

$\approx p_0(\tau) \cdot (1 - \lambda \Delta \tau)$

From this we can write

$p_0(\tau + \Delta \tau) - p_0(\tau) = - \lambda \Delta \tau \cdot p_0(\tau)$

which is leads to the differential equation

$\frac{dp_0(\tau)}{d\tau} = - \lambda \cdot p_0(\tau)$

with solution

$p_0(\tau) = e^{-\lambda \tau)$

I guess the main difference between this and what I did is $t_2 \to \tau-\Delta \tau$.

mathtest 3

\begin{equation}
\sqrt{\sigma}
\end{equation}

math test 2

Math test 2
$\sqrt(\sigma)$
\\begin{equation}
\sqrt(\sigma)
\\end{equation}


\begin{equation}
\sqrt(\sigma)
\end{equation}