## Have my old job!

My old job at the Health & Safety Laboratory is being advertised, and at a higher pay grade to boot. (Though it is still civil service pay, and thus not going to make you rich.)

You’ll need to have solid mathematical modelling skills, particularly solving systems of ODEs, and be proficient at writing scientific code, preferably R or MATLAB or acslX. From chats with a few people at the lab, management are especially keen to get someone who can bring in money so grant writing and blagging skills are important too.

It’s a smashing place to work and the people are lovely. Also, you get flexitime and loads of holiday. If you are looking for a maths job in North West* England then I can heartily recommend applying.

*Buxton is sometimes North West England (when we get BBC local news) and sometimes in the East Midlands (like when we vote in European elections).

## Viewing the internals of MATLAB Matrices

A cool undocumented trick I just learnt from The MathWorks’ Bob Gilmore. If you type

format debug

Then printing any vector reveals information about its internal representation. For example:

x = magic(3) x = Structure address = 6bc1ab0 m = 3 n = 3 pr = d8dccf0 pi = 0 8 1 6 3 5 7 4 9 2

The structure address is the address in memory where the matrix is stored, `m`

and `n`

are the number of rows and columns respectively of the matrix, and `pr`

and `pi`

are pointers to the addresses of the matrices storing the real and imaginary components of the matrix.

One interesting thing to look at is the representation of scalar numbers.

y = 1 y = Structure address = 6bc31e0 m = 1 n = 1 pr = d790b90 pi = 0 1

Yep: they are stored in exactly the same way as matrices: in the same way the “everything in R is a vector”, everything in MATLAB is a matrix. To finish up, here are some more examples for you to explore:

% higher dimensional arrays rand(2, 3, 4) % cell arrays (unfortunately not that revealing) {1, magic(3)} % sparse matrices (very interesting) sparse(ones(3))

## MATLAB’s stand out new feature

It’s been a while since my last MATLAB post, not because I don’t love the language, but more because I do most of my blogging from home, where I have no license, and because (mostly thanks to R-bloggers) I get ten times as many page views for the R posts. (TODO: Create MATLAB-bloggers service.)

Having returned from holiday (it was lovely, thanks for asking) I’ve been trying out the latest release of MATLAB – R2011b. So far, the standout new feature is the automatic variable renaming. If you change the name of a variable at the point where it was declared, then pressing Shift+Enter lets MATLAB rename all other instances. IDEs for statically-typed languages have had this feature for years, but to see it in a dynamically-typed language is very impressive.

## Friday Function: nclass

When you draw a histogram, an important question is “how many bar should I draw?”. This should inspire an indignant response. You didn’t become a programmer to answer questions, did you? No. The whole point of programming is to let your computer do your thinking for you, giving you more time to watch videos of fluffy kittens.

Fortunately, R contains three functions to automate the answer, namely `nclass.Sturges`

, `nclass.scott`

and `nclass.FD`

. (FD is short for Freedman-Diaconis; watch out for the fact that `scott`

isn’t capitalised.)

The differences depend upon length and spread of data. For longer vectors, Scott and Freedman-Diaconis tend to give bigger answers.

short_normal <- rnorm(1e2) nclass.Sturges(short_normal) #8 nclass.scott(short_normal) #8 nclass.FD(short_normal) #12

long_normal <- rnorm(1e5) nclass.Sturges(long_normal) #18 nclass.scott(long_normal) #111 nclass.FD(long_normal) #144

For strongly skewed data, you are best to use some sort of transformation before you draw a histogram, but for the record, Freedman-Diaconis again gives bigger answers for highly skewed (and thus wider) vectors.

short_lognormal <- rlnorm(1e2) nclass.Sturges(short_lognormal) #8 nclass.scott(short_lognormal) #9 nclass.FD(short_lognormal) #20

long_lognormal <- rlnorm(1e5) nclass.Sturges(long_lognormal) #18 nclass.scott(long_lognormal) #443 nclass.FD(long_lognormal) #1134

My feeling is that since each of the three algorithms is rather dumb, it is safest to calculate all three, then pick the middle one.

nclass.all <- function(x, fun = median) { fun(c( nclass.Sturges(x), nclass.scott(x), nclass.FD(x) )) } log_islands hist(log_islands, breaks = nclass.all(log_islands))

I also wrote a MATLAB implementation of this a couple of years ago.

It is worth noting that ggplot2 doesn’t accept a number-of-bins argument to `geom_histogram`

, because

In practice, you will need to use multiple bin widths to

discover all the signal in the data, and having bins with

meaningful widths (rather than some arbitrary fraction of the

range of the data) is more interpretable.

That’s fine if you are interactively exploring the data, but if you want a purely automated solution, then you need to make up a number of bins.

calc_bin_width <- function(x, ...) { rangex <- range(x, na.rm = TRUE) (rangex[2] - rangex[1]) / nclass.all(x, ...) } p <- ggplot(movies, aes(x = votes)) + geom_histogram(binwidth = calc_bin_width(log10(movies$votes))) + scale_x_log10() p

## supercalifragilisticexpialidocious = 1

I notice that the latest version of R has upped the maximum length of variable names from 256 characters to a whopping 10 000! (See `?name`

.) It makes the 63 character limit in MATLAB look rather pitiful by comparison. Come on MathWorks! Let’s have the ability to be stupidly verbose in our variable naming!

## Presentation on testing frameworks

It’s testing week in the Software Carpentry course. To commemorate the occasion, I’ve uploaded a presentation I gave a couple of years ago to evangelise the use of testing frameworks to make your life easier. Just for you, I’ve recorded an all new audio track to accompany the slides. The example uses MATLAB but the principle is applicable to any language.

## When 1 * x != x

Trying to dimly recall things from my maths degree, it seems that in most contexts the whole point of the number one is that it is a multiplicative identity. That is, for any x in your set, 1 * x is equal to x. It turns out that when you move to floating point numbers, *in some programming lanugages*, this isn’t always true.

In R, try the following

x <- Inf + 0i 1 * x [1] Inf + NaNi

Something odd is happening, and my first reaction was that this is a bug. In fact, I reported it as such (in a slightly different form) but the terrifyingly wise Brian Ripley came up with the cause.

When you multiply two complex numbers, this is what happens

In this instance

So, there is a reasonable argument that R is doing the right thing.

MATLAB works a little differently since all it’s numbers are implicitly complex. The inner workings of MATLAB are somewhat opaque for commercial reasons, but in order for compiled MATLAB code to work, I think it’s reasonable to assume that MATLAB stores its variables in a class that is similar to the MWArray class used by compiled code. As far as I can tell, each `double`

matrix contains two matrices for storing real and imaginary components, and has an `IsComplex`

property that is true whenever the complex part has any nonzero values. If `IsComplex`

is false, only the real matrix is used for arithmetic.

Thus in MATLAB, when you define `x = Inf + 0i`

it immediately resolves itself to simply `Inf`

, and we don’t get the same weirdness.

x = Inf + 0i 1 * x Inf

Other languages are divided on the subject: python and ruby agree with R but Mathematica (or at least Wolfram Alpha) agrees with MATLAB.

Thinking about this from a purely mathematical sense,

This concurs with the MATLAB answer and I’m inclined to agree that it makes more sense, but the issue isn’t entirely clear cut. Changing the way complex arithmetic works for a single edge case is likely to introduce more confusion than clarity. It is perhaps sufficient for you to remember that complex infinity is a bit pathological with arithemtic in some languages, and add checks in your code if you think that that will cause a problem.