Stupid Google Calculator Tricks

How many fortnights would it take light to travel the circumference of the earth?

What about the speed of light in furlongs per fortnight?

They've got the mass of the proton and the mass of the electron so you can easily calculate their ratio.

I guess we've no need for that gigantic CRC Handbook any more...but what if you need to look up a constant and you're on a desert island with no internet connection, better still pack that CRC just to be safe.

SCons for Latex

Latex is really great for publishing papers with nicely typeset equations and graphics. It automates all of the little formatting stuff and lets you spend time developing the content of your paper rather than wasting effort fiddling with format. I've been using GNU Make for years to process my Latex sources into dvi's and pdf's, but it's always been pretty painful and clumsy (mostly because I'm not a Make guru).

The makefiles for a Latex document can get pretty ugly and complicated especially if your document depends on various data or graphics files. SCons is a new and improved build-tool that implements its configuration files as Python scripts, the file is called SConstruct now instead of Makefile. The great thing from my point of view is that it has built-in support for Latex. Making a pdf from my Latex source file is simple:

env = Environment()
env.PDF('main.tex')

The PDF builder is smart enough to know it needs to run pdflatex and then based on the output possibly run bibtex and then pdflatex again to track down undefined references (or not). All in two lines of Python.

The really neat thing is how easy it is to extend the built-in build environments.

env = Environment()
plot_bld = Builder(action = 'gnuplot < $SOURCE > $TARGET,
suffix = '.png'
src_suffix = '.gp')
env.Append(BUILDERS = {'Plot' : plot_bld})
env.PDF('main.tex')

Now my SCons can automatically update my png graphics files based on gnuplot commands in any files ending with .gp.

Practical Development Environments has a nice little write-up in Chapter 5 about SCons.

EM Flash Cards

I just discovered a great new plug-in (in time for mid-terms) for OpenOffice.org Impress called Open Cards. It is uses the titles of the slides as the question and the body of the slide as the answer to the flashcards. It keeps track of your learning state, has a couple of different learning modes, and seems to use some pretty sophisticated training algorithms.

I put together a set of flashcards for electrodynamics if you're interested in trying it out (of course you need Open Office installed to use the file).

Complex Step

Finite differences are cool, but you are limited by subtractive cancellation (section 3.4: subtractive cancellation exercises). What if you want to make a numerical estimate of the derivative of a function, like a 2D Gaussian:


A finite difference approximation would be (in Octave):

dfdx = ( f(x+h) - f(x) )/ h;

Complex step is even cooler, because you don't have any subtraction (there's no difference), so you can choose a very small step size without loosing accuracy due to subtractive cancelation:

complex_step = complex( 0, 1e-320 );

The derivative is approximated by just the imaginary part:

x = x + complex_step;
dfdx = imag( f(x) )/imag( complex_step );

So the derivative with respect to x looks like this:

And the estimated derivative with respect to y:


This approach is really useful for design sensitivity analysis, and since modern Fortran supports complex types we can even use this method for serious number crunching!

2008 Election Outcome

This is not a political blog, but the Economist has a pretty good article out on John McCain that I felt like commenting on. The last line sums it up: "...pity, there are few better men in American politics".

While he is a very different man than his fellow Naval Academy grad Jimmy Carter, both are men of high integrity, who took that academy training in "duty, honor, country" seriously, and lived it in their public lives after their military service. Though I wasn't alive for the Carter presidency, so I only know it was ineffective from reading about it, I have been alive for President Carter's active work with Habitat for Humanity and his attempts at peace brokering around the world. I think it's safe to say that he embodies Christian ideals of service and peace-seeking. While John McCain is not a paragon of Christian values, he is an exemplar of honesty and sacrifice in public service.

It's too bad that such good men are such unsuccessful politicians in our country. Hopefully Senator McCain's public service after his unsuccessful presidential bid will be as good for the world as President Carter's continues to be following his unsuccessful presidency. Either way, they are a couple of Academy grads I'm proud of (even if they did go to the wrong academy).

Go Air Force, Sink Navy, Beat Army!

A winning combination: Intel and Heineken

I just got a little Atom-based Intel board to play around with and it currently resides in a rubermaid container as its case. I've also been enjoying the new Heineken draught kegs, so I have a few empties sitting in my kitchen floor. This led to the obvious question:

Can a little Atom motherboard

fit inside a Heineken Draught keg?

These new Atom based boards are little mini-itx form factor, that's ~6.93x6.93". The little micro ATX power supplies are ~5x4x3" and a hardrive is ~3.5x5x1". So it might be just a tad to close on the diameter to be able to slide the motherboard right in.

That would be a pretty sweet case though...

The Shuffle

The shuffle, also known as the Fisher's Exact Test, is a permutation test that can be used to estimate the sampling distribution of a statistic without relying on parametric assumptions. This is especially important when sample sizes are small. The other neat thing about permutation tests is that you don't have to know what the distribution of your statistic is. So if you have a really odd function of your data that you want to use rather than one of the classical statistics you can.

Octave has a great built-in function called nchoosek which makes shuffling a breeze. Called with scalar arguments it returns the value of the binomial coefficient, which is the number of ways you can choose k things from n things (k<=n). For fixed n, nchoosek is maximum when k=n/2. That is plotted below for n=2:24.

Notice that's a log-log scale, as n increases it quickly becomes intractable to perform a complete permutation test.

Suppose we have two data vectors, and we want to know if they are from populations with different means. We can use the normrnd() function to draw 8 samples from a normal distribution with 0 mean and 8 samples from a normal distribution with a mean of 1.

n=8;
x1 = normrnd(0,1,n,1);
x2 = normrnd(1,1,n,1);
perms = nchoosek( [x1;x2], n );

We stored all 12870 permutations of the vectors (16 choose 8) in the matrix perms. Now we use the built-in function complement to find the difference of the means under each of those labellings.

for(i=1:length(m1) )
comps(i,:) = complement( perms(i,:), [x1;x2] );
endfor
m1 = mean( perms, 2 ); % take row-wise averages, DIM=2
m2 = mean( comps, 2 ); % take row-wise averages, DIM=2
m_shuff = m1 - m2;

Now we can use the distribution of m_shuff to decide if the difference in the means of our two data vectors is significant.


Octave script with all of these calculations: shuffle.m.

Simple Bayes

I like Bayes theorem, it's really useful. The most intuitive and accessible explanation I've found of using Bayes theorem to solve a problem is in Russell and Norvig's classic, Chapter 20 (pdf) (I just own the first edition, the second edition looks even better).

The initial example they give is about pulling different flavoured candy out of a sack (remember the balls and urn from your basic stats?). They also provide a really good discussion showing how standard least-squares regression is a special case of maximum-likelihood for when the data are generated by a process with Gaussian noise of fixed variance.

Their first example is for estimating parameters in a discrete distribution of candy, but we can apply the same math to estimating the variance of a continuous distribution. Estimating variance is important, lots of times in industrial or business settings the variance of a thing matters as much or more than the average, just check-out all of the press those Six Sigma guys get. That's because it gives us insight into our risk. It helps us answer questions like, "What's our probability of success?" And maybe, if we're lucky, "What things is that probability sensitive to?"

Bayes theorem is a very simple equation:

Where P(h) is the prior probability of the hypothesis, P(d|h) is the likelihood of the data given the hypothesis, and P(h|d) is the posterior probability of the hypothesis given the data.

Octave has plenty of useful built-in functions that make it easy to play around with some Bayesian estimation. We'll set up a prior distribution for what we believe our variance to be with chi2pdf(x,4), which gives us a Chi-squared distribution with 4 degrees of freedom. We can draw a random sample from a normal distribution with the normrnd() function, and we'll use 5 as our "true" variance. That way we can see how our Bayesian and our standard frequentist estimates of the variance converge on the right answer. The standard estimate of variance is just var(d), where d is the data vector.

The likelihood part of Bayes theorem is:

% likelihood( d | M ) = PI_i likelihood(d_i, M_j)
for j=1:length(x)
lklhd(j) = prod( normpdf( d(1:i), 0, sqrt( x(j) ) ) );
endfor
lklhd = lklhd/trapz(x,lklhd); % normalize it

Then the posterior distribution is:

% posterior( M | d ) = prior( M ) * likelihood( d | M )
post_p = prior_var .* lklhd;
post_p = post_p/trapz(x,post_p); % normalize it

Both of the estimates of the variance converge on the true answer as n approaches infinity. If you have a good prior, the Bayesian estimate is especially useful when n is small.

It's interesting to watch how the posterior distribution changes as we add more samples from the true distribution.


The great thing about Bayes theorem is that it provides a continuous bridge from what we think we know to reality. It allows us to build up evidence and describe our knowledge in a consistent way. It's based on the fundamentals of basic probability and was all set down in a few pages by a nonconformist Presbyterian minister and published after his death in 1763.

Octave file for the above calculations: simple_bayes.m

Is Octave Slow?

The answer, like most things worth asking the question about, is "Well, it depends..."

I'm not going to talk about "computer time" vs. "programmer time"; though that's probably one of the most important considerations most of the time (check out Paul Graham's essays if you need convincing). I'll just talk about what takes a long time to compute in Octave, and what trade-offs can be made to improve the situation. We'll dwell mostly on the first tip from the Octave Manual (you'll see it's first for good reason). The other important thing to do is always allocate memory before looping over a vector (eg. x=zeros(n,1)).

Suppose we have a simple re-sampling problem. We would like to estimate the sampling distribution of a statistic of our data. The data could be a simple vector of 1 or 0 representing success or failure of a trial, and we want to estimate the probability of success.

x = [1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0];

One method to solve the problem would be to use a for loop:

X = empirical_rnd(N,x); % get a bootstrap sample from x
for(j=1:10) % timing loop, measure a couple times to get an idea of the
% variation
tic();
for( i=1:N-n ) % loop over the bootstrap samples
p(i) = sum( X(i:i+n-1) )/n;
endfor
fortime(j)=toc();
endfor

With N=20000, this takes 1.09 seconds on my old laptop, that seems kind of slow. The sampling distribution for the probability of success is shown in the figure below.

It is always a good idea to vectorize what you can when in Octave, so the method below does just that.

X = empirical_rnd(x,N,n);
for(j=1:10)
tic();
P = sum(X,2)/n;
vecttime(j) = toc();
endfor

With n=20000, this takes 0.009 seconds on the same old laptop, that's a speed-up of 2 orders of magnitude!

Well, those are trivial examples of vectorizing a calculation to avoid the dreaded for-loop. What happens when subsequent calculations depend on the results of previous iterations?

tic();
p1(1) = 0;
for(i=2:N)
t = t + dt;
p1(i) = p1(i-1)+dt*2*t;
endfor
toc()

This loop integrates x*x numerically, with N=20000 it takes ~0.75 seconds. We can't vectorize, so how do we speed it up? Octave has some good sparse matrix capabilities, maybe we could recast the problem as a sparse matrix solve.

Now we are trading memory for speed. In the for-loop implementation we just have to store the vector p. If we want to use Octave's sparse matrix facilities we need to store the two diagonals of our operator, so that roughly triples the memory requirements. Given the enormous size of modern computer memory, most toy problems should fit (if your problem doesn't fit, why are you still using Octave?).

tic();
A = spdiag( ones(N,1), 0 );
A = A + spdiag( -ones(N-1,1), -1 );
p2 = A\(dt*2*t);
toc()

The direct solution of the simple bi-diagonal system, with N=20000 takes ~0.019 seconds, better than 2 orders of magnitude speed-up over the for-loop implementation. For more complex sparse operators one of the iterative schemes might be appropriate.


The moral: if your problem is small enough to fit into memory, cram it all in and don't use for-loops!

Is Octave slow? Well, it depends, if you use for-loops it is.

Further Reading

A couple more posts on optimizing your Octave code.

• More on Speeding Up Octave (24 Dec 2008)
  
• How to Speed Up Octave (5 Jan 2009)