Getting Started

Phasic is a library for modelling and inference on Markov jump processes, which move through multiple states on the the way to a final (absorbing) state. Using algorithms that translate existing matrix-based results for phase-type distributions into a graph framework, Phasic efficiently computes exact moments and distributions of waiting times.

Installing

Phasic is most easily installed using either Pixi, Conda or Pip:

pixi workspace channel add munch-group
pixi install phasic
conda install -c munch-group -c conda-forge phasic
pip install phasic

Your First Model

The simplest possible model has a single state in where it stays with process—you start, wait some random amount of time, then finish. Each states is represented by a vector (list/array) of integers each representing a state property. In this simple case, we only need a single integer.

The easiest way to create that (and any other model) is to write a “callback” function. Think of the function as providing the answer to this question: “Assuming some state in your Markov process, what other states can the process jump to and with what rates?”. When passed a state vector, it must return a (possibly empty) list of child states that the process can jump to, paired with the rates at which these occur.

If you can think about your favorite Markov model in this way, phasic will do the rest for you.

from phasic import Graph
import numpy as np

def model(state):
    transitions = []               # List for transitions from state
    if state[0] < 2:               # Control transitions allowed
        child = state.copy()       # Make a copy to modify
        child[0] += 1              # Make the child state
        trans = [child, 2.0]       # Pair the child with a transition rate (in a list not tuple)
        transitions.append(trans)  # Append transition 
    return transitions

You can now create the graph by passing the callback function to the Graph constructor, along the ipv keyword argument specifying in which state the process begins. Then plot the graph to make sure your callback function works as expected:

graph = Graph(model, ipv=[1])
graph.plot()

This models an exponential distribution with rate 2.0. Once we’ve constructed the graph, we can compute moments directly. The expected waiting time is 1/rate = 0.5, which we verify by calling expectation() on the graph:

graph.expectation(), graph.variance()
(0.5, 0.25)
graph.moments(10)
[0.5, 0.5, 0.75, 1.5, 3.75, 11.25, 39.375, 157.5, 708.75, 3543.75]

You can compute probability distributions. The pdf() method returns the probability density at a given time, while cdf() returns the probability that the waiting time is less than or equal to that time:

import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 2, figsize=(7, 3))
vals = np.linspace(0, 3, 100)
axes[0].plot(vals, graph.pdf(vals))
axes[1].plot(vals, graph.cdf(vals))
plt.show()

You can also draw random samples from the distribution:

samples = graph.sample(10000)
samples[:10]
array([0.09573634, 0.74686978, 0.09943768, 0.4057115 , 0.15184864,
       0.22272146, 0.51297259, 0.04698131, 1.2134511 , 0.02203194])

Sample mean vs. expectation:

np.mean(samples).item(), graph.expectation()
(0.5012805569042454, 0.5)

Using rewards, you can compute expectations and PDFs of marginal waiting times, time spent in each state and much more - all as easily as for this simple example. You an learn more about this, as well as about Bayesian inference in the tutorials. For now, just appreciate how easily complex models can be constructed:

from phasic import Graph
import numpy as np

def mesh(state, max_val=None):
    transitions = []
    for i in range(state.size):
        if state[i] <= max_val:
            child = state.copy()
            child[i] += 1
            trans = [child, state.sum()]
            transitions.append(trans)
    return transitions

graph = Graph(mesh, ipv=[1, 1], max_val=2)
graph.plot()

Note how any keyword arguments are passed onto the callback function:

graph = Graph(mesh, ipv=[1, 1], max_val=8)
print(graph.expectation())
graph.plot()
1.3225441834186733