Bayesian estimation of changepoints

A common introductory problem in Bayesian changepoint detection is the record of UK coal mining disasters from 1851 to 1962. More information can be found in Carlin, Gelfand and Smith (1992). As we can see from the plot below, the number of yearly disasters ranges from 0 to 6 and we will assume that at some point within this time range a change in the accident rate has occured.

The number of yearly disasters can be modelled as a Poisson with a unknown rate depending on the changepoint \(k\):

$$ y_t \sim \text{Po}\left(\rho\right),\qquad \rho = \begin{cases} \mu, & \text{if}\ t=1,2,\dots,k \\ \lambda, & \text{if}\ t = k +1, k + 2, \dots,m \end{cases} $$

Our objective is to estimate in which year the change occurs (the changepoint \(k)\) and the acident rate before (\(\mu\)) and after (\(\lambda\)) the changepoint amounting to the parameter set \(\Phi = \left\lbrace\mu,\lambda,k\right\rbrace\). We will use Crystal (with crystal-gsl) to perform the estimation.

We start by placing independent priors on the parameters:

  • \(k \sim \mathcal{U}\left(0, m\right)\)
  • \(\mu \sim \mathcal{G}\left(a_1, b_1\right)\)
  • \(\lambda \sim \mathcal{G}\left(a_2, b_2\right)\)

For the remainder we'll set \(a_1=a_2=0.5\), \(c_1=c_2=0\) and \(d_1=d_2=1\). The joint posterior of \(\Phi\) is then:

$$ \pi\left(\Phi|Y\right) \propto p\left(Y|\Phi\right) \pi\left(k\right) \pi\left(\mu\right) \pi\left(\lambda\right), $$

where the likelihood is

$$ \begin{align} p\left(Y|\Phi\right) &= \prod_{i=1}^{k} p\left(y_i|\mu,k\right) \prod_{i=k+1}^{m} p\left(y_i|\lambda,k\right) \\ &= \prod_{i=1}^{k} \frac{\mu^{y_i}e^{-\mu}}{y_i!} \prod_{i=k+1}^{m} \frac{\lambda^{y_i}e^{-\lambda}}{y_i!}. \end{align} $$

As such, the full joint posterior can be written as:

$$ \begin{align} \pi\left(\Phi|Y\right) &\propto \prod_{i=1}^{k} \frac{\mu^{y_i}e^{-\mu}}{y_i!} \prod_{i=k+1}^{m} \frac{\lambda^{y_i}e^{-\lambda}}{y_i!} \left(\mu^{a_1-1} e^{-\mu b_1}\right) \left(\lambda^{a_2-1} e^{-\lambda b_2}\right) \frac{1}{m} \\ &= \mu^{a_1 + \sum_{1}^{k}y_i - 1}e^{-\mu\left(k+b_1\right)} \lambda^{a_2 + \sum_{k+1}^{m}y_i - 1}e^{-\lambda\left(m-k+b_2\right)} \end{align}. $$

It follows that the full conditionals are, for \(\mu\):

$$ \begin{align} \pi\left(\mu|\lambda,k,Y\right) &\propto \mu^{a_1 + \sum_{i=1}^{k}y_i-1}e^{-\mu\left(k+b_1\right)} \\ &= \mathcal{G}\left(a_1+\sum_{i=1}^{k}y_i, k + b_1\right) \end{align} $$

We can define the \(\mu\) update as:

def mu_update(data : Array(Int), k : Int, b1 : Float64) : Float64
  Gamma.sample(0.5 + data[0..k].sum, k + b1)

The full conditional for \(\lambda\) is:

$$ \begin{align} \pi\left(\lambda|\mu,k,Y\right) &\propto \lambda^{a_2 + \sum_{i=k+1}^{m}y_i-1}e^{-\lambda\left(m-k+b_2\right)} \\ &= \mathcal{G}\left(a_2+\sum_{i=k+1}^{m}y_i, m - k + b_2\right), \end{align} $$

which we implement as:

def lambda_update(data : Array(Int), k : Int, b2 : Float64) : Float64
  Gamma.sample(0.5 + data[(k+1)..M].sum, M - k + b2)
The next step is to take $$ \begin{align} b_1 &\sim \mathcal{G}\left(a_1 + c_1,\mu + d_1\right) \\ b_2 &\sim \mathcal{G}\left(a_2 + c_2,\lambda + d_2\right), \end{align} $$ which we will implement as:
def b1_update(mu : Float64) : Float64
  Gamma.sample(0.5, mu + 1.0)

def b2_update(lambda : Float64) : Float64
  Gamma.sample(0.5, lambda + 1.0)

And finally we choose the next year, \(k\), according to

$$ p\left(k|Y,\Phi\right)=\frac{L\left(Y|\Phi\right)}{\sum_{k^{\prime}} L\left(Y|\Phi^{\prime}\right)} $$


$$ L\left(Y|\Phi\right) = e^{\left(\lambda-\mu\right)k}\left(\frac{\mu}{\lambda}\right)^{\sum_i^k y_i} $$

implemented as

def l(data : Array(Int), k : Int, lambda : Float64, mu : Float64) : Float64
  Math::E**((lambda - mu)*k) * (mu / lambda)**(data[0..k].sum)

So, let's start by writing our initials conditions:

iterations = 100000

b1 = 1.0
b2 = 1.0

M = data.size # number of data points

# parameter storage
mus = Array(Float64).new(iterations, 0.0)
lambdas = Array(Float64).new(iterations, 0.0)
ks = Array(Int32).new(iterations, 0)

We can then cast the priors:

mus[0] = Gamma.sample(0.5, b1)
lambdas[0] = Gamma.sample(0.5, b2)
ks[0] =

And define the main body of our Gibbs sampler:

(1...iterations).map { |i|

  k = ks[i-1]

  mus[i] = mu_update(data, k, b1)
  lambdas[i] = lambda_update(data, k, b2)

  b1 = b1_update(mus[i])
  b2 = b2_update(lambdas[i])

  ks[i] = Multinomial.sample((0...M).map { |kk| 
    l(data, kk, lambdas[i], mus[i])


Looking at the results, we see that the mean value of \(k\) is 38.761, which seems to indicate that the change in accident rates occured somewhere near \(1850+38.761\approx 1889\). We can visually check this by looking at the graph below. Also plotted are the density for the accident rates before (\(\mu\)) and after (\(\lambda\)) the change.

Of course, one the main advantages of implementing the solution in Crystal is not only the boilerplate-free code, but the execution speed. Compared to an equivalent implementation in R the Crystal code executed roughly 17 times faster.

Language Time (s)

t as mixture of Normals

(Based on Rasmus Bååth's post)

A scaled \(t\) distribution, with \(\mu\) mean, \(s\) scale and \(\nu\) degrees of freedom, can be simulated from a mixture of Normals with \(\mu\) mean and precisions following a Gamma distribution:

$$ y \sim \mathcal{N}\left(\mu,\sigma\right) \\ \sigma^2 \sim \mathcal{IG}\left(\frac{\nu}{2},s^2\frac{\nu}{2}\right) $$

Since I've recently pickep up again the crystal-gsl in my spare time, I've decided to replicate the previously mentioned post using a Crystal one-liner. To simulate 10,000 samples from \(t_2\left(0,3\right)\) using the mixture, we can then write:

samples = (0..10000).map { |x| 
   Normal.sample 0.0, 1.0/Math.sqrt(Gamma.sample 1.0, 9.0) 

We can see the mixture distribution (histogram) converging nicely to the \(t_2(0,3)\) (red):

Langton's Ant

Last week, at the North East Functional Programming meet up, we were given a code Kata consisting of the Langton's ant algorithm. I've had a go at Scala but decided later on to put a live version in this blog. I considered several implementation options, such as scala.js and Elm, but in the end decided to implement it in plain Javascript.

A Gibbs sampler in Crystal

Recently, I've been following with interest the development of the Crystal language.

Crystal is a statically type language with a syntax resembling Ruby's. The main features which drawn me to it were its simple boilerplate-free syntax (which is ideal for quick prototyping), tied with the ability to compile directly to native code along with a dead simple way of creating bindings to existing C code. These features make it quite attractive, in my opinion, for scientific computing.

To test it against more popular languages, I've decided to run the Gibbs sampling examples created in Darren Wilkinson's blog. I recommend reading this post, and in fact, if you are interested in Mathematics and scientific computing in general, I strongly recommend you follow the blog.

As explained in the linked post, I will make a Gibbs sampler for $$f\left(x,y\right)=kx^2\exp\left\lbrace-xy^2-y^2+2y-4x\right\rbrace$$ with \[ x|y\sim Ga\left(3,y^2+4\right) \\ y|x\sim N\left(\frac{1}{1+x},\frac{1}{2\left(1+x\right)}\right) \]

The original examples were ran again, without any code alterations. I've just added the Crystal version. This implementation uses a very simple wrapper I wrote to the famous GNU Scientific Language (GSL).

require "../libs/gsl/"
require "math"

def gibbs(n : Int = 50000, thin : Int = 1000)
  x = 0.0
  y = 0.0
  puts "Iter  x  y"
  (0..n).each do |i|
    (0..thin).each do |j|
      x = Statistics::Gamma.sample(3.0, y*y+4.0)
      y = Statistics::Normal.sample(1.0/(x+1.0), 1.0/Math.sqrt(2.0*x+2.0))
    puts "#{i} #{x} #{y}"


(As you can see, the Crystal code is quite similar to the Python one).

To make sure it's a fair comparison, I ran it in compiled (and optimised) mode build using

crystal build --release
time ./gibbs > gibbs_crystal.csv

Looking at the results, you can see that they are consistent with the other implementations:

The timings for each of the different versions (ran in a 1.7 GHz Intel Core i7 Macbook Air) were

Language Time (s)

So there you have it. A Ruby-like language which can easily compete with C performance-wise.

I sincerely hope that Crystal gets some traction in the scientific community. That of course won't depend solely on its merits but rather on an active community along with a strong library ecosystem. This is lacking at the moment, simply because it is relatively new language with the specs and standard library still being finalised.

Scala types map

Scala has a strong static type system which arguably can be considered complex.

For my own benefit, this is an attempt at mapping Scala's 2.11.7 type hierarchy interactively.

Unlike Java where you have primitive types which step out of the hierarchy, Scala has a unified type system where all types descend from Any (a "top type", ⊤, in type theory), which in the standard library has two direct subtypes (AnyRef and AnyVal). At the other end of type hierarchy1 (commonly called the "bottom types", ⊥) we have Null$ and Nothing$. In essence all of Scala's remaining library types will live between ⊤ and ⊥. The nodes in the graph are coloured according to its category, i.e. class, trait, object or case class, aliases in the standard library are kept with their nominal value and there is no distinction between class and abstract class.

When a type is selected, the tree paths will represent all the upwards reachable nodes, this is, its immediate successors (and successors' successors) will be highlighted up the hierarchy.

You can select types using the URL. To pre-select a type and it's ascendants go to, for instance:

If you have any suggestions or comments, please leave them at @ruimvieira on Twitter.

A map of Scala types

1 - One of the existential reasons for ⊥ is closely related to the Turing Halting Problem and how to fit it within a consistent strong static type system.

Without going into too much detail (perhaps a future blog post), these types exist partly to represent the result of computations for which the Scala compiler can't easily deduce the result type. Let's assume we have a never ending computation represented by the function foo:

def foo():T = ???
What should T be? While in many other languages (notoriously, Java) we could get away with typing foo with just about anything (for instance, String or TransactionAwarePersistenceManagerFactoryProxy), Scala is more strict, since we can't verify that assertion. The solution within our type system is simply to declare
def foo():Nothing = ???
Of course that in Java you could write void interpreted as "return nothing". However, to a functional programmer, having a function that returns nothing (which is quite different from Unit) is surely a code smell.


There is now a new section called Ponies.

Ponies is a deliberately simple snippet and recipe repository motivated by the fact that I usually don't remember most of the API details and idiosyncrasies, and, even worse, could never bring myself to create a system of quickly finding those recurring little tasks which end up being googled/stack overflowed/random favourite forum-ed.

Admittedly, the first question that springs to mind is "WYARR (Why Yet Another Recipe Repository)?". No special reason, apart from being able to cherry pick my favourite methods, whenever an API doesn't enforce TOOWTDI ("There's Only One Way To Do It")

The second question might be "Why is it called Ponies?".
Simply because it's a kind of string of one-trick-ponies and "string" would eventually be a (even) more confusing name.

It is just starting and, as mentioned, I'll use it as personal a note holder. Nevertheless, please feel free to use it if you find it useful.
It supports basic title search as well as language tag searching. Using it as

list lang:scala
will search for (you've guessed it) all the recipes mentioning "list" in the title and tagged with "scala".

Full-text search is not implemented but, time permitting, will be eventually.

Happy coding.

MCMC notifications

It is said that patience is a virtue but the truth is that no one likes waiting (especially waiting around: this interesting article explores why people prefer walking 8 minutes to the airport’s baggage claim and having the bags ready rather than waiting the same amount of time entirely in the claim area).

Anyone performing computationally heavy work, such as Monte Carlo methods, will know that these are usually computationally expensive algorithms which, even in modern hardware, can result in waiting times in the magnitude of hours, days and even weeks. These long running times coupled with the fact that in certain cases it is not easy to accurately predict how long a certain number of iterations will take, usually leads to a tiresome behaviour of constantly checking for good (or bad) news. Although it is perfectly possible to specify that your simulation should stop after a certain amount of time (especially valid for very long simulations), this doesn’t seem to be the standard practice.

In this post I’ll detail my current setup for being notified exactly of when simulations are finished. To implement this setup, the following stack is required:

  • A JDK
  • Apache Maven
  • A messaging service (Pushbullet will be used)
  • A smartphone, tablet, smartwatch (or any other internet enabled device)

To start, we can create an account in Pushbullet, which will involve, in the simplest case, signing up using some authentication service such as Google.

Next, we will install the client application (available for Android, iOS and most modern browsers) after which we can enable notifications (at least in the Android client, I’m not familiar with the iPhone version).

Since my current work started as a plain Java project which in time evolved mainly to Scala, it consists of an unholy mixture of Maven as a build tool for Scala code. This shouldn’t be a problem for other setups, but I’ll just go through my specific setup (i.e. using Maven dependencies to a Scala project).

To implement communication between the code and the messaging service, we can use a simple library such as jpushbullet.
The library works well enough, although at the time of writing it only supports Pushbullet’s v1 API but not the newer v2 API. Since the project, unfortunately, is not in Maven central, you should build it from scratch. Fortunately, in a sensibly configured machine, this is trivial.

In the machine where you plan to perform the simulations, clone and build jpushbullet.

  git clone
  mvn clean install

Once the build is complete, you can add it as a local dependency in your project’s pom.xml:


For the purpose of this example, lets assume that you have the following Object as the entry point of your simulation. The next step is to add a call to the Pushbullet service before the exit point. Please keep in mind that it is very bad practice to include your personal API key in your committed code. I strongly suggest you keep this information in a separate file (e.g. in resources), read it at runtime and add it to .gitignore. That being said, place the messaging code as such:

package benchmarks

import com.shakethat.jpushbullet.PushbulletClient

object MCMC {

  def main(args:Array[String]):Unit = {

    // Your MCMC code

    val client = new PushbulletClient(api_key)
    val devices = client.getDevices
    val title = "MCMC simulation finished"
    val body = "Some summary can be included"
    // n is the preferred device number
    client.sendNote(true, devices.getDevices.get(n).getIden(), title, body)


Usually, I would call this code via ssh into my main machine from Maven (using Scala Maven) as:

  nohup mvn scala:run -DmainClass=benchmarks.MCMC &

Finally, when the simulation is completed, you will be notified in the client devices (you can select which ones by issuing separate sendNote calls) and include a result summary, as an example.

My current setup generates an R script from a template which is run by Rscript in order to produce a PDF report. However, be careful, since file quotas in Pushbullet are limited, so text notifications should be used without worry of going over the free usage tier.

Keep in mind, that there are other alternatives to jpushbullet, such as send-notification, a general notification library for Java for which the setup is quite similar.

Hope this was helpful.


I have been working in a new library called gulp which you can on
On the project\'s page there are some usage examples but I will try to summarise the main points here.
'The purpose of this library is to facilitate the parallel development of R and Java code, using rJava as the bridge. Creating bindings in rJava is quite simple, the tricky part of the process (in my opinion) being the maintenance of the bindings (usually done by hand) when refactoring your code.
As an example, let's assume you have the following Java class
public class Test {
  // Java code
That you wish to call from R.

New version

Finally the rewriting of the website is going at a good speed and basic functional version is available.
This site is being rewritten from scratch, and there a few thing you might notice
In this blog I don't intend to implement any tagging system, but try to implement any context handling by implementing an above average search system.
There is also a RESTful API which can be consulted here. It's in a very early stage, so please expect things to change.
And if you know your web frameworks, why not try to find out in which this site was written? ;).
You can drop your answer on twitter (check the about page).
Stay tuned for further developments.)