Zen Coding: A Markov chain in Clojure

I’m currently in the middle of, among learning many new things, getting back into Clojure. I never really could imagine deploying production-grade Clojure programs, but writing code in a functional paradigm is always a fascinating, challenging, and rewarding exercise.

Plus, any Lisp looks pretty cool syntactically, and such languages are always extremely terse.

Here’s a Clojure adaptation of the HIV/AIDS Markov chain scenario model that I wrote about and implemented in Python yesterday. Look at how short it is!

(require '[clojupyter.misc.helper :as helper])
(helper/add-dependencies '[net.mikera/core.matrix "0.62.0"])
(use '[clojure.core.matrix :as m])

(def hiv [[0.9 0.07 0.02 0.01] [0 0.93 0.05 0.02] [0 0 0.85 0.15] [0 0 0 1]])
(def p-0 [[0.85, 0.1, 0.05, 0]])

(defn evolutions [p-matrix initial steps]
  (reductions m/mmul
              initial
              (repeat steps p-matrix)))

(defn future-distribution [p-matrix initial steps]
  (last
    (evolutions p-matrix
                initial
                steps)))

(defn watch-evolve [column p-matrix initial steps]
  (-> (evolutions p-matrix
                  initial
                  steps)
      (get-column 0)
      (get-column column)))

Note: It’s written in a way that works in Hydrogen for Atom (i.e., via a hidden Jupyter Notebook). I’m not sure how you’d go about importing core.matrix in a one-off REPL session (not a Leiningen project), but now that I’ve discovered Hydrogen, I doubt there will be many terminal REPL sessions in my future!

Coding in Clojure means non-object-oriented and no mutable state. You can write a function to compute the probability distribution at time $t$, but the system won’t stay in that state.

In fact, there is no system—only evolution. In that sense, you quite literally end up with a memoryless Markov chain. The design philosophy of Clojure is all about impermanence, which is very Zen, when you think about it. (Would that make each program a koan?)

Let’s take a quick tour of the functions.

First, there is evolutions [p-matrix initial steps]. That’s a base function that works as a building block for the other functions. In pseudo-math, it would look like this:

$$ E(P, \pi^{(0)}, t) = {\pi^{(1)}, \cdots, \pi^{(t)}} $$

It’s easy to see that computing some $\pi^{(n)}$, the probability distribution (state of the system) at a given time, can be accomplished with a single reduce.

But Clojure has an interesting function reductions (docs), which stores all of the intermediate outputs in a list, with the expected output of reduce as the last entry. So it makes sense to write two functions here:

one, evolutions, that gives us all $\{\pi^{(1)}, \cdots, \pi^{(t)}\}$, and
another, future-distribution, that takes the same inputs and outputs only the last evolution.

In pseudo-math:

$$ F(P, \pi^{(0)}, t) = \pi^{(t)} $$

In action:

(future-distribution hiv p-0 20)
; => [[0.10334015640198398 0.24687062185910383 0.12037223090781371 0.5294169908310989]]

What else can we do with the output of evolutions? Well, we can choose one state of the system and watch its distribution evolve over time by picking out the $k$-th item of each evolution (a row vector).

In this case, we can observe how the proportion of a segment of the population (asymptomatic, symptomatic, has AIDS, dead) changes over $t$ time-steps.

In pseudo-math:

$$ W(s, P, \pi^{(0)}, t) = \Big\{ P(s_0), P(s_1 | s_0), \cdots, P(s_t | s_{t-1}) \Big\} $$

In action:

(watch-evolve 0 hiv p-0 5)
; => [0.85 0.765 0.6885 0.61965 0.5576850000000001 0.5019165000000001]

What a cool language.

Go Top