Andrew Birkett's blog

In the previous post, I talked about Bayesian stats and MCMC methods in general. In this post, I’ll work through an example where we try to infer how fair a coin-toss is, based on the results of ten coin flips. Most people use JAGS via an R interface, but I’m going to use JAGS directly to avoid obfuscation.

(Note: a coin-toss is a physical event determined by physics, so the “randomness” arises only through uncertainty of how hard it’s tossed, how fast it spins, where it lands etc, and therefore is open to all sorts of evil)

Firstly, we have to tell JAGS about our problem – eg. how many coin tosses we’ll do, and that we believe each coin toss is effectively a draw from a Bernoulli distribution with unknown proportion theta, and what our prior beliefs about theta are.

To do this, we create “example.model” containing:

model {
  for (i in 1:N){
    x[i] ~ dbern(theta)
  }
  theta ~ dunif(0,1)
}

This says that we’ll have N coin-flips, and each coin flip is assumed to be drawn from the same Bernoulli distribution with unknown proportion theta. We also express our prior belief that all values of theta from zero to one are equally likely.

We can now launch “jags” in interactive mode:

$ jags
Welcome to JAGS 4.2.0 on Sun Feb 26 14:31:57 2017
JAGS is free software and comes with ABSOLUTELY NO WARRANTY
Loading module: basemod: ok
Loading module: bugs: ok

.. and tell it to load our example.model file ..

. model in example.model

If the file doesn’t exist, or the model is syntactically invalid you’ll get an error – silence means everything has gone fine.

Next, we need the data about the coin flip, which corresponds to the x[1] .. x[N] in our model. We create a file called “example.data” containing:

N < - 10
x <- c(0,1,0,1,1,1,0,1,0,0)

The format for this file matches what R’s dump() function spits out. Here we’re saying that we have flipped ten coins (N is 10) and the results were tails/heads/tails/heads/heads etc. I’ve chosen the data so we have the same number of heads and tail, suggesting a fair coin.

We tell JAGS to load this file as data:

. data in example.data
Reading data file example.data

Again, it’ll complain about syntax errors (in an old-school bison parser kinda way) or if you have duplicate bindings. But it won’t complain yet if you set N to 11 but only provided 10 data points.

Next, we tell JAGS to compile everything. This combines your model and your data into an internal graph structure, ready for evaluating. It’s also where JAGS will notice if you’ve got too few data points or any unbound names in your model.

. compile
Reading data file example.data
. compile
Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 10
   Unobserved stochastic nodes: 1
   Total graph size: 14

The graph consists of ten “observed” nodes (one per coin flip) and one unobserved stochastic node (the unknown value of theta). The other nodes presumably include the bernoulli distribution and the uniform prior distribution.

At this stage, we can tell JAGS where it should start its random walk by providing an initial value for theta. To do this, we create a file “example.inits” containing:

theta < - 0.5

.. and tell JAGS about it ..

. parameters in example.inits
Reading parameter file example.inits

Finally, we tell JAGS to initialize everything so we’re ready for our MCMC walk:

. initialize
Initializing model

Now we’re ready to start walking. We need to be a bit careful at first, because we have to choose a starting point for our random walk (we chose theta=0.5) and if that’s not a good choice (ie. it corresponds to a low posterior probability) then it will take a while for the random walk to dig itself out of the metaphorical hole we dropped it in. So, we do a few thousand steps of our random walk, give it a fancy name like “burn-in period” and cross our fingers that our burn-in period was long enough:

. update 4000
Updating 4000
-------------------------------------------------| 4000
************************************************** 100%

(JAGS give some enterprise-level progress bars when in interactive mode, but not in batch mode).

JAGS has happily done 4000 steps in our random walk, but it hasn’t been keeping track of anything. We want to know what values of theta is jumping between, since that sequence (aka “chain”) of values is what we want as output.

To tell JAGS to start tracking where it’s been, we create a sampler for our ‘theta’ variable, before proceeding for another 4000 steps, and then writing the results out to a file:

. monitor theta
. update 4000
-------------------------------------------------| 4000
************************************************** 100%
. coda *

The last command causes two files to be written out – CODAindex.txt and CODAchain1.txt. CODA is a hilariously simple file format, coming originally from the “Convergence Diagnostic and Output Analysis” package in R/S-plus. Each line contains a step number (eg. 4000) and the value of theta at that step (eg. 0.65).

Here’s an interesting thing – why would we need a “Convergence Diagnostic” tool? When we did our “burn-in” phase we crossed our fingers and hoped we’d ran it for long enough. Similarly, when we did the random walk we also used 4000 steps. Is 4000 enough? Too many? We can answer these questions by looking at the results of the random walk – both to get the answer to our original question, but also to gain confidence that our monte-carlo approximation has thrown enough darts to be accurate.

At this point, we’ll take our coda files and load them into R to visualize the results.

$ R
R version 3.0.2 (2013-09-25) -- "Frisbee Sailing"
Copyright (C) 2013 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

> require(coda)
Loading required package: coda

> c < - read.coda(index.file="CODAindex.txt",output.file="CODAchain1.txt")
Abstracting theta ... 5000 valid values

> summary(c)
Iterations = 4001:9000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 5000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

          Mean             SD       Naive SE Time-series SE 
      0.501658       0.139819       0.001977       0.001977 

2. Quantiles for each variable:

  2.5%    25%    50%    75%  97.5% 
0.2436 0.4000 0.5022 0.6017 0.7675 

This is telling us that, given ten coin flips and our prior uniform belief and our bernoulli assumption, the most probably value for theta (the proportion of coin-flip yielding heads) is close to 0.5. Half of the probability mass lies between theta=0.4 and theta=0.6, and 95% of the probability mass lies between theta=0.25 and theta=0.75.

So it’s highly unlikely that the coin flip is extremely biased – ie. theta<0.25 or theta>0.75. Pleasantly, “highly unlikely” means “probability is less than 5%”. That’s a real common-or-garden probability. Not any kind of frequencist null-hypothesis p-value. We can make lots of other statements too – the probability that the bias is greater than 0.75 is about 40%. If we had a second coin (or coin flipper) we could make statement like “the probability that coin2 has a higher bias than coin1 is xx%”.

Let’s briefly revisit the question of convergence. There’s a few ways to determine how well your random walk represents (or “has converged to”) the true posterior distribution. One way, by Rubin and Gelman, is to run several random walks and look at the variance between them. The coda package in R comes with a function gelman.diag() for this purpose. However, in our simple example we only did one chain so we can’t run it on our coda files. (Incidentally, Gelman writes a great blog about stats).

In the next post, I’m will look at the performance characteristics of JAGS – how it scales with the number of data points, and what tools you can use to track this.

JAGS is a useful statistics tool, helping you decide how to generalise experimental results. For example, if you roll a die ten times and it comes up “six” half of the time, is that strong evidence that the die is loaded? Or if you toss a coin ten times and it comes up heads nine times, what the probability that the coin is a “normal” coin?

JAGS is based on the bayesian approach to statistics, which uses Bayes rule to go from your experimental results to a (probabilistic) statement of how loaded the dice is. This is different approach from the frequencist approach to statistics which most textbooks cover – with p values and null hypothesis test. The upside of the bayesian approach is that it answers the kind of questions you want to ask (like, “what is the probability that drug A is better than drug B at treating asthma”) as opposed to the convoluted questions which frequencist statistics answer (“assuming that there’s no difference between drug A and drug B, what’s the probability that you’d get a measured difference at least as large as the one you saw?”). The downside of Bayesian stats is that you have to provide a “prior” probability distribution, which expresses your beliefs of how likely each outcome is prior to seeing any experiment results. That can seem a bit ugly, since it introduces a subjective element to the calculation. Some people find that unacceptable, and indeed that’s why statistics forked in the early 1900s from its bayesian origins to spawn the frequencist school of probability with it’s p-values driven by popular works by Ronald Fisher. But on the other hand, no experiment is run in a vacuum. We do not start each experiment in complete ignorance, nor is our selection of which experiments to run, or which hypotheses to check, determined objectively. The prior allows us to express information from previous knowledge, and we can rerun our analysis with a range of priors to see how sensitive our results are to the choice of prior.

Although Bayes rule is quite simple, only simpler textbook examples can be calculated exactly using algebra. This does include a few useful cases, like the coin-flipping example used earlier (so long as your prior comes from a particular family of probability distribution). But for more real-world examples, we end up using numerical techniques – in particular, “Markov Chain Monte Carlo” methods. “Monte Carlo” methods are anything where you do simple random calculations which, when repeated enough times, converge to the right answer. A nice example is throwing darts towards a circular darts board mounted on a square piece of wood – if the darts land with uniform probability across the square, you can count what fraction land inside the circle and from that get an approximation of Pi. As you throw more and more darts, the approximation gets closer and closer to the right answer. “Markov Chain” is the name given to any approach where the next step in your calculation only depends on the previous step, but not any further back in history. In Snakes and Ladders, your next position depends only on your current position and the roll of the die – it’s irrelevant where you’ve been before that.

When using MCMC methods for Bayesian statistics, we provide our prior (a probability distribution) and some data and a choice of model with some unknown parameters, and the task is to produce probability distributions for these unknown parameters. So our model for a coin toss might be a Bernoilli distribution with unknown proportion theta, our prior might be a uniform distribution from 0 to 1 (saying that we think all values of theta are equally likely) and the data would be a series of 0’s or 1’s (corresponding to heads and tails). We run our MCMC algorithm of choice, and out will pop a probability distribution over possible values of theta (which we call the ‘posterior’ distribution). If our data was equally split between 0’s and 1’s, then the posterior distribution would say that theta=0.5 was pretty likely, theta=0.4 or theta=0.6 fairly likely and theta=0.1 or theta=0.9 much less likely.

There’s several MCMC methods which can be used here. Metropolis-Hasting, created in 1953 during the Teller’s hydrogen bomb project, works by jumping randomly around the parameter space – always happy to jump towards higher probability regions, but will only lump to lower probability regions some of the time. This “skipping around” yields a sequence (or “chain”) of values for theta drawn from the posterior probability distribution. So we don’t ever directly get told what the posterior distribution is, exactly, but we can draw arbitrarily many values from it in order to answer our real-world question to sufficient degree of accuracy.

JAGS uses a slightly smarter technique called Gibbs Sampling which can be faster because, unlike Metropolis-Hasting, it never skips/rejects any of the jumps. Hence the name JAGS – Just Another Gibbs Sampler. You can only use this if it’s easy to calculate the conditional posterior distribution, which is often the case. But it also frees you from the Metropolis-Hasting need to have (and tune) a “proposal” distribution to choose potential jumps.

In the next post, I’ll cover pragmatics of running JAGS on a simple example, then look at the performance characteristics.

One of the nice things about scmutils is that it tracks units, so you can’t accidentally add 10 seconds to 5 metres.

(+ 
 (& 10 &second)
 (& 5 &meter))
=> Units do not match: + (*with-units* 10 (*unit* SI ... 1)) (*with-units* 5 (*unit* SI ... 1))

When dealing with angles, it initially seems to do the right thing too:

(+
 (& pi/2 &radian)
 (& 90 °ree))
=> 3.141... (ie. its converting everything to radians)

But this is less cool:

(/ (& pi &radian) (& 1 &second))
=> (& 3.141592653589793 &hertz)

Err, pi radians should be 0.5Hz. The trouble is, scmutils treats radians as a unit-less number.

To check whether this was a reasonable thing to do, I checked my old favourite Frink. In frink’s units.txt files, we have the following:

// Alan’s editorializing:
// Despite what other units programs might have you believe,
// radians ARE dimensionless units and making them their own
// unit leads to all sorts of arbitrary convolutions in
// calculations (at the possible expense of some inclarity if
// you don’t know what you’re doing.)
// If you really want radians to be a fundamental unit,
// replace the above with “angle =!= radian”
// (This will give you a bit of artificiality in calculations.)
//
// The radian was actually a fundamental base unit in the SI
// up until 1974, when they changed it, making it no longer
// be a special unit, but just a dimensionless number (which
// it is.) See the definition of the “Hz” below for a
// discussion of how this broke the SI’s definitions of
// basic circular / sinusoidal measures, though.

And down a bit, on the section about hertz, we have:

//
// Alan’s Editorializing: Here is YET ANOTHER place where the SI made a
// really stupid definition. Let’s follow their chain of definitions, shall
// we, and see how it leads to absolutely ridiculous results.

// The Hz is currently defined simply as inverse seconds. (1/s).
// See: http://physics.nist.gov/cuu/Units/units.html
//
// The base unit of frequency in the SI *used* to be “cycles per second”.
// This was fine and good. However, in 1960, the BIPM made the
// change to make the fundamental unit of frequency to
// be “Hz” which they defined as inverse seconds (without qualification.)
//
// Then, in 1974, they changed the radian from its own base unit in the SI
// to be a dimensionless number, which it indeed is (it’s a length divided by
// a length.) That change was correct and good in itself.
//
// However, the definition of the Hz was *not* corrected at the same
// time that the radian was changed. Thus, we have the conflicting SI
// definition of the radian as the dimensionless number 1 (without
// qualification) and Hz as 1/s. (Without qualification.)
//
// This means that, if you follow the rules of the SI,
// 1 Hz = 1/s = 1 radian/s which is simply inconsistent and violates basic
// ideas of sinusoidal motion, and is simply a stupid definition.
// The entire rest of the world, up until that point, knew that 1 Hz needs to
// be equal to *2 pi* radians/s or be changed to mean *cycles/second* for
// these to be reconcilable. If you use “Hz” to mean cycles/second, say,
// in sinusoidal motion, as the world has done for a century, know that the SI
// made all your calculations wrong. A couple of times, in different ways.
//
// This gives the wonderful situation that the SI’s Hz-vs-radian/s definitions
// have meant completely different things in the timeperiods:
//
// * pre-1960
// * 1960 to 1974
// * post-1974
//
//
// Thus, anyone trying to mix the SI definitions for Hz and angular
// frequencies (e.g. radians/s) will get utterly wrong answers that don’t
// match basic mathematical reality, nor match any way that Hz was ever used
// for describing, say, sinusoidal motion.
//
// Beware the SI’s broken definition
// of Hz. You should treat the radian as being correct, as a fundamental
// dimensionless property of the universe that falls out of pure math like
// the Taylor series for sin[x], and you should treat the Hz as being a
// fundamental property of incompetence by committee.
//
// One could consider the CGPM in 1960 to have made the original mistake,
// re-defining Hz in a way that did not reflect its meaning up to that point,
// or the CGPM in 1974 to have made the absolutely huge mistake that made
// the whole system inconsistent and wrong, and clearly broke the definition
// of Hz-vs-radian/s used everywhere in the world, turning it into a broken,
// self-contradictory mess that it is now.
//
// Either way, if I ever develop a time machine, I’m going to go back and
// knock both groups’ heads together. At a frequency of about 1 Hz. Or
// better yet, strap them to a wheel and tell them I’m going to spin one group
// at a frequency of 1 Hz, and the other at 1 radian/s and let them try to
// figure out which one of those stupid inconsistent definitions means what.
// Hint: It’ll depend on which time period I do it in, I guess, thanks to
// their useless inconsistent definition changes.
//
// It’s as if this bunch of geniuses took a well-understood term like “day”
// and redefined it to mean “60 minutes”. It simply breaks every historical
// use, and present use, and just causes confusion and a blatant source of
// error.
//
// In summary: Frink grudgingly follows the SI’s ridiculous, broken definition
// of “Hz”. You should not use “Hz”. The SI’s definition of Hz should be
// considered harmful and broken. Instead, if you’re talking about circular
// or sinusoidal motion, use terms like “cycles/sec” “revolutions/s”,
// “rpm”, “circle/min”, etc. and Frink will do the right thing because it
// doesn’t involve the stupid SI definition that doesn’t match what any
// human knows about sinusoidal motion. Use of “Hz” will cause communication
// problems, errors, and make one party or another look insane in the eyes
// of the other.

I’m doing the Stanford “Quantum Physics for Engineers” online course just now. Separately, a few months ago I was reading the Sussman “Structure And Interpretation of Classical Mechanics” book which is notable for using scheme as a mathematical notation, thereby avoiding a lot of the ambiguities of ‘normal’ maths notation (a big win in Lagrangian mechanics, which makes heavy use of partial derivatives).

Anyhow, the Stanford Quantum course requires you to do various exercises, such as the following:

An electron has a 1nm wavelength. Is it reasonable to treat this electron as an approximately non-relativistic particle (i.e. traveling much slower than the speed of light)?

As usual, this requires plugging the supplied numbers and a bunch of physics constants into the right equation. At school, I would’ve done this by hand – hopefully remembering constants like ‘c’ (3e8 m/s) and h (6.62e-34).

But I can also do this using scheme, as per the SICM book. The ‘scmutils’ library comes with a bunch of built-in constants, with the correct units:

:c
=> (& 299792458. (* &meter (expt &second -1)))

:h
=> (& 6.62606896e-34 (* (expt &meter 2) &kilogram (expt &second -1)))

In scmutils, the ampersand function attaches units to a number.

So now I can use de Broglie’s wavelength relation to find velocity as a function of mass and wavelength:

(define (velocity mass wavelength) (/ :h (* mass wavelength)))

then plug in the appropriate values to find the velocity:

(velocity :m_e (& 1e-9 &meter))
=> (& 727389.4676462485 (* &meter (expt &second -1)))

The question actually asked “can you treat it as non-relativistic” so we want to know if it’s close to the speed of light or not:

(/ (velocity :m_e (& 1e-9 &meter)) :c)
=> 2.43e-3

So it’s much slower than the speed of light, and the answer is “yes, it’s reasonable to treat this as a non-relativistic particle). But thanks to scheme/scmutils, I’m also pretty confident I haven’t made errors with units (because scheme tracked them for me) or constants (because I didn’t have to enter them).

Although not required for this exercise, the scmutils package also handles symbolic differentiation which is pretty nifty! For example:

(define (foo x) (log x))

(foo 'a)
 => (log a)

((D foo) 'x)
 => (/ 1 x)

The scmutils library is very elegant once you realise how it works. The definition of the scheme ‘foo’ function is just that – a scheme function. You can use it in one of two ways. You can pass a number to it – eg. (foo 5) – and it’ll evaluate it numerically – eg. 1.609. Or you can pass that same function a symbol, such as ‘a, and it’ll give you back a symbolic expression – eg. “log a”. It has a built-in simplifier too, as seen here:

(define (addaddadd x) (+ x x x))
=> #| addaddadd |#

(addaddadd 'a)
=> #| (* 3 a) |#

A while ago, I wrote an emacs ‘reading mode’. It highlights a single sentence at a time, fading the rest of the text into a gentle grey, and a keypress moves onto the next sentence. It retains the familiarity and consistency of normal text layout, but provides additional cues about the extent of the current sentence.

Tonight, I played with the idea of including smarter parsing into this reading mode. The Stanford Parser parses english sentences. It tells you about the grammatical structure (noun phrases, verb phrases, etc) and dependencies between words. This is just about enough to do what I had in mind – a “superfluous word” highlighter. The whole world is absolutely packed full of so many documents with wholly unnecessary words. Ideally, I’d like to just delete the pointless words. But it’s rare for a word to be completely devoid of semantic meaning. So, my compromise is just to highlight those decorative words – adjectival and adverbial modifiers – which are commonly guilty.

Here’s some examples, not completely perfect, but useful nonetheless:

I REALLY want some SUPER TASTY chocolate.
The system has been VERY CAREFULLY designed, and will cope admirably with all 
  CONCEIVABLE combinations of circumstances.
I wanted to leave my SMALL pond and see HOW I'd fare in a BIG one, with some 
  of the BEST developers in the world.
You define HOW you want your data to be structured ONCE, THEN you can 
  use SPECIAL GENERATED source code to EASILY write and read your STRUCTURED data.