Drawing random numbers from a multivariate Gaussian distribution

I like software libraries. A software library means that someone has already done some of the heavy lifting for me. A software library means that I won’t have to reinvent the wheel. It means that I might not have to spend my afternoon fixing bugs and stringently validating my own code*.

It also means that I don’t necessarily have a good idea of what the computer is doing “behind the scenes”; that I’m potentially unaware of the clever (or not) algorithm that is being implemented, and what it’s limitations (or bugs) are.

For a long time now, I’ve been blithely accepting that numpy.random.multivariate_normal applies some sensible, but mysterious, procedure in order to generate random variables from a multivariate Gaussian (normal) distribution. And it does – I just never tried to unravel said mystery. So, prompted by Dag Sverre, and courtesy of Wikipedia,  I present a Python snippet that draws multivariate Gaussians given (a) some way of generating a set of standard univariate Gaussian random numbers, and (b) a Cholesky decomposition routine. In all, it’s about 10% faster than the NumPy routine, which applies a bunch of validity tests that I don’t.

import numpy as np

# Define mean, covariance matrix, and no. of samples required
mean = ... # some array of length N
cov = ... # some positive-definite NxN matrix
Ndraws = 1000

# Do factorisation (can store this and use it again later)
L = np.linalg.cholesky(cov)

# Get 3*Ndraws Gaussian random variables (mean=0, variance=1)
norm = np.random.normal(size=Ndraws*3).reshape(3, Ndraws)

# Construct final set of random numbers (with correct mean)
rand = mean + np.dot(L, norm)


Of course, this now introduces two new mysteries: How do you efficiently implement Cholesky decomposition, and how do you generate univariate Gaussian random numbers? Those are exercises for another day, I think.

* It depends on how sloppy I’m feeling; of course, you should validate everything that matters, regardless of who wrote it.

A good ol’ fashioned seminar tour

I’m in Oslo at the moment, on the first leg of my World Tour Of Physics Departments. My supervisor, Pedro, is keen on sending his students out to give talks as they near the final year of their PhD, and it’s increasingly the “done thing” in the US when postdoc application season comes around too, apparently. As far as I’m concerned, it’s an opportunity to meet lots of interesting people and promote/develop my research, so the “job tour” concept gets two thumbs up from me. I’ll be in Helsinki and Geneva next week, followed by Lyon and Heidelberg the week after. I’ll also be visiting Johns Hopkins, Pittsburgh, Princeton, Stanford, Berkeley, and finally Bielefeld in late October/early November. Following that, I’ll probably have to wander off to plant a few trees and apologise to Mother Earth for my carbon footprint.

As for Oslo, I’ve spent the last few days working with my (spectacularly generous) host, Hans Kristian Eriksen, and others on a soon-to-be-published paper on the Sunyaev-Zel’dovich Effect. Watch this space. I gave my seminar talk today (it was about the acceleration stuff I did with Tim Clifton), and I’ve also had a number of interesting discussions with Cython magician Dag Sverre Seljebotn about a possible rewrite of Hans Kristian’s Commander CMB analysis code. But most importantly: I saw a moose at the weekend. Mission accomplished, that’s what I say.

Book review: The Artist’s Guide to GIMP

I just reviewed The Artist’s Guide to GIMP over at my open source blog. Overall, it’s a nice introduction to the free image editing powerhouse that is the GIMP – the book is packed full of interesting projects that seem pretty easy to follow, written in a refreshingly waffle-free manner. It would make a good read for web developers, computer artists, and Photoshop whizzes interested in moving over to the right side of the proprietary software divide alike (although the bits on photography have some rough edges).

See here for the full review.

Healpix map2gif options

I’m doing a lot of work with CMB maps at the moment, most of which are pixelised using Healpix. It’s good to be able to visualise the maps, and the Healpix distribution comes with a utility to do just that: map2gif. Unfortunately, map2gif has somewhat sparse documentation – there are a few options that aren’t documented and that sort of thing. I ended up poking through the source code to find some options, so I figured I should share them here to make other peoples’ lives a little easier.

Rough full-sky map of the Sunyaev-Zel’dovich effect, extracted from WMAP data.

(You can access much of this information by typing map2gif -hlp to print usage instructions, but they’re quite basic.)

Basic usage

At its simplest, you need to supply filenames for an input FITS file (in the Healpix format) and an output GIF:

$ map2gif -inp input.fits -out map.gif

If you have a CMB map, you’re probably most interested in seeing the CMB anisotropies. Unfortunately, residual foregrounds can mess up the scaling of the map, so you have to add a scale yourself. For a map in units of micro-Kelvin, the following scale should produce good results for the anisotropies:

$ map2gif -inp input.fits -out map.gif -min -300 -max 300

For most uses, you probably want a colour scale to be displayed on the output image too. Use the -bar option for this (it’s a boolean value, so specify True):

$ map2gif -inp input.fits -out map.gif -bar True

A rough (Gnomonic-projected) map of the Sunyaev-Zel’dovich effect centred around the Virgo cluster, extracted from WMAP data.

Projection and centering

For a lot of applications, you’ll be wanting to zoom in on some region of the map rather than seeing the full sky. Healpix supports the Mollweide (default) and Gnomonic projections, which you can switch between using the -proj switch. You only need to type the first three letters of each projection. If you want to zoom in on some specific area, use Gnomonic and specify a different map origin using the (FITS file coordinate system-dependent) -lat and -lon switches:

$ map2gif -inp input.fits -out map.gif -proj gno -lat 45 -lon 87

By default this will look quite zoomed-in, so you might want to use the resolution (-res) switch to zoom out a little (higher numbers zoom out more):

$ map2gif -inp inp.fits -out map.gif -proj gno -lat 45 -lon 87 -res 4

Presentation (colours, labels, image size)

There are a number of switches for controlling the presentation of the output image, most of which are explained sufficiently well in the Healpix documentation. For example, you can add a title to the plot (-ttl), increase its width in pixels (-xsz; the aspect ratio is kept the same) or use a different colormap (-col; try different numbers to get different colour maps).

Other options

There are a few other options that are a bit more specialised; for example, you can scale the signal by arcsinh (-ash) for some reason. A particularly interesting one for CMB applications is the map selection option (-sig); it’s typical to include maps for several Stokes parameters (I, Q, and U) inside a single FITS file for experiments that do polarisation mapping, so you can use this switch to select between them. For example, to pick the Q map out of a WMAP IQU map, use:

$ map2gif -inp input.fits -out map.gif -sig 2

There are also some basic rescaling operations that you can apply using the addition (-add) and multiplication (-mul) switches. Read more about those in the Healpix docs.

Listing of command line options

Here’s a full listing of command-line options to map2gif, taken directly from the output of the -hlp switch:

usage: MAP2GIF
[-inp input_file]
[-out output_file]
[-max usermax] [-min usermin]
[-add offset] [-mul factor]
[-pro projection] [-res gridresn]
[-lon lon0] [-lat lat0]
[-col color_table] [-sig signal]
[-xsz xsize] [-log logflg]
[-ash ashflg]
[-bar add color bar] [-ttl title]
[-hlp]

Weird gfortran compiler error of the week

In response to the following compiler error from gfortran, the correct response is apparently to run make clean and then make again:

Error: More actual than formal arguments in procedure call at (1)

God knows what was actually going wrong.