Gaussian Processes for Machine Learning

What a great resource for learning about Gaussian Processes: The Gaussian Processes Web Site.

• The Gaussian Processes for Machine Learning book.
  
• Software around the web, and to go with the book
  
• Data Sets
• Tutorials
  

Scientific Computing Going Mobile/Embedded

I follow both the Octave and Maxima user mailing lists (I learn more from the Maxima list, but the Octave one still throws up a gem every now and then). I am also a founding member of the Dayton Diode hackerspace. So it was interesting to see some things of interest to makers/hackers pop-up on those two lists.

The first was a thread on the Octave list about installing Octave and the Octave Forge packages on a Raspberry Pi. In theory it should be as easy as any other Debian system using the package manager for the distribution. No word back from the OP yet on whether it works in practice or not.

The second thread to pop up was on the Maxima list about installing Maxima on Android systems. Maxima-on-Android is apparently a bundle on Google Play of Maxima 5.28.0, Gnuplot 4.6, Mathjax 2.1 and Qepcad 1.69. I have not played with it yet on any of my Android devices.

Mathworks Fighting the GPL

From the Octave mailing list:

from Jordi Gutiérrez Hermoso to Octave-help date Wed, Feb 10, 2010 at 5:55 PM subject [OT]: Mathworks fighting the GPL mailing list Filter messages from this mailing list unsubscribe Unsubscribe from this mailing-list

Apologies for the offtopic chatter, but I just noticed this as I was
trying to look for Emacs's Matlab mode to see if I could get some
ideas for Octave mode. It appears that the Mathworks only allows the
BSD license on their servers now, which affects much free software
under other licenses (in particular GPL) that they were hosting.

It's difficult to find exactly what happened and how because they seem
to have spread the word about this relatively secretly by email to
people with code on their servers, but they do have a FAQ:

    http://www.mathworks.com/matlabcentral/FX_transition_faq.html

This happened July last year.

I don't know about anyone else, but I'm deeply disturbed that
Mathworks is now actively fighting against copyleft.

- Jordi G. H.

This seems to be so they can transition to a  more 'cloud-like' model.  From the FAQ:

Why is the File Exchange adding licensing?

Licensing clarifies the rights you have as an author and as a user of the code available on the File Exchange. Licensing details how the file can be used and addresses common questions around rights to modification, distribution, and commercial use.

What happens if I don't do anything? Do I have to put a license on my code?

We have no plans to remove unlicensed submissions, but they will be prominently marked as unlicensed. In addition, unlicensed contributions will not be available for use with any future tools that access the File Exchange from within MathWorks products.

Well being able to download and use code on the Exchange 'automagically' without worrying about license restrictions is probably a good thing.  Of course this also means MathWorks can then ship binaries based on your code without sharing the source.  BSD vs. GPL: which is freedomier, flame on!

JASSM Bootstrap Reliability II

I previously did an analysis of JASSM's reliability using Octave. It's a very simple exercise using the built-in Octave function empirical_rand() and public statements by various folks concerning the flight tests.

As reported by Reuters, JASSM has had success in recent tests. To be quantitative, 15 out of 16 test flights were a success. Certainly good news for a program that's been bedevilled by flight test failures and Nunn-McCurdy breaches.

The previous analysis can be extended to include this new information. This time we'll use Python rather than Octave. There is not (that I know of) a built-in function like empirical_rand() in any of Python's libraries, but it's relatively straightforward to use random integers to index into arrays and accomplish the same thing.

def bootstrap(data, nboot):
"""Draw nboot samples from the data array randomly with replacement, ie
a bootstrap sample."""
bootsample = sp.random.random_integers
return(data[bootsample(0, len(data)-1, (len(data), nboot))])

Applying this function to our new data vectors to get bootstrap samples is easy.

# reliability data for July 2009, based on Reuters report
d09_jul_tests = 19
d09_jul = sp.ones(d09_jul_tests, dtype=float)
d09_jul[0:3] = 0.0

# generate the bootstrap samples:
d09_jul_boot = sp.sum(bootstrap(d09_jul, nboot), axis=0) / float(d09_jul_tests)
# find the number of unique reliabilities to set the number of bins for the histogram:
d09_jul_nbins = sp.unique(d09_jul_boot).size

So, is it a 0.9 missile or a 0.8 missile? We should probably be a little more modest in the inferences we wish to draw from such small samples. As reported by Bloomberg the Air Force publicly contemplated cancellation for a missile with reliability distributions like the blue histogram shown below (six of ten failed), while stating that if 13 of 16 were successful (the green histogram) that would be acceptable. The figure below shows these two reliability distributions along with the actual recent performance.

This seems to be a reasonably supported inference from these sample sizes, the fixes that resulted in the recent 15 out of 16 successes had a measurable effect when compared to the 4 out of 10 performance.

Not that binary outcomes for reliability are a great measure, they just happen to be easily scrape-able from press releases. The figure below illustrates the problem, there just isn't much info in each 1 or 0, so really cranking up the number of samples only slowly improves the power (reduces the area of overlap between the two distributions).

JASSM bootstrap reliability

The Joint Air-to-Surface Stand-off Missile has been in the news lately. It's the same old story that's been going on ever since it was "fielded". A few missiles fail to function properly in testing, the Air Force says Lockheed-Martin needs to bring up the reliability or the program will be terminated. Lockheed-Martin agrees that more work needs to be done ("keep the money flowing") to improve the reliability, and they'll work diligently with their Air Force team members to serve the warfighter.

That's military-industrial complex business as usual, not that interesting. The neat part is that there's enough information in most of those press releases and news articles to do some resampling statistics on the reliability of the missile. So, in the interest of supporting "an alert and knowledgeable citizenry" (see the Eisenhower video at the bottom of this post) here are some bootstraps on JASSM reliability.

It is really easy to do in Octave using the empirical_rnd() function.

n_09 = 19; % number of tests, based on stated 79% success rate
f_09 = 4; % number of failures
nbins_09 = 14;
t_09 = ones(n_09, 1);
t_09(1:f_09) = 0;
boot_09 = empirical_rnd(t_09, n_09, nboot); % really easy to do
% bootstraps
reliability_09 = sum(boot_09, 1) / n_09;

The sample size n_09 and number of failures f_09 are based on the recent Reuters story:

Four JASSM missiles tested in November, January and February did not detonate on impact or had other problems, raising fresh questions about the program. But the missile still had a reliability rate of 79 percent, and was on track to reach the 90 percent goal, the Air Force said. --Reuters, 6 Jul 2009

So, according to the recent reports JASSM should work about four times out of five (the AF wants it to work nine times out of ten). Here's the "4 out of 5" reliability distribution (based on the bootstrap shown above) with 19 samples:

Back in 2004 JASSM had a claimed reliability of 76%:

We have had 29 launches of JASSM and we have a 76 percent success rate. --Judy Stokely, deputy of acquisition at the Air Armament Center

Which means it should work three times out of four. Here's a "3 out of 4" missile's reliability distribution with 29 samples:

The sampling distributions are too large to measure a change as small as the difference between 0.76 reliability and 0.79 reliability, so it's statistically the same missile now that it was back in 2004. In both distributions the desired 0.9 reliability is out on the tail of the distributions, i.e. you can't claim it's a "9 out of 10" missile with much credence.

A more interesting question is what sort of sampling distribution would you get from nine successes and one failure (the desired nine out of ten missile)?

With a sample size of only ten it would be pretty hard to tell the difference between a "9 out of 10" missile and a "4 out of 5" missile. Based on the press releases apparently a "4 out of 5" missile is unacceptable, but what level of confidence will the Air Force place on knowing that they have 0.9 reliability?

Eisenhower on the military industrial complex:


Here's the Octave file with the code for doing the bootstraps.

Also, thanks to Michael J.T. O'Kelly's bootstrap.py for showing how easy it is to resample with replacement from an array in Python using SciPy.

Open Math Tools Custom Search Engine

I found a really simple how-to for creating custom Google search engines and used it to make my own little Open Source Math Tools search engine. I find it really annoying when I'm searching for documentation or examples in Maxima that the results are full of Nissan Maxima pages. Things are even worse when trying to search for stuff relevant to R.

Using the custom search engine I don't have to specify that I'm searching for Maxima, Octave, R or SciPy stuff, because I already told Google to favor results from those pages. Right now the custom search engine (CSE) only has about twenty sites included here's the big ones (list updated 10 Nov 2009):

• ntrs.nasa.gov
  
• mit.edu
  
• umich.edu
  
• projects.scipy.org/pipermail/*
  
• nabble.com
  
• nist.gov
  
• http://www.csulb.edu/~woollett/*
  
• http://beshenov.ru/maxima/*
  
• axiom-developer.org
  
• clusterbuilder.org
  
• linuxhpc.org
  
• clustermonkey.net
  
• matplotlib.sourceforge.net
  
• f2py.org
  
• http://cens.ioc.ee/projects/f2py2e/*
  
• scipy.org
  
• http://maxima.sourceforge.net/email-archives/*
  
• http://www.math.utexas.edu/pipermail/maxima/*
  
• maxima.sourceforge.net
  
• www.gnu.org/software/octave
  
• octave.org
  
• r-project.org
  

And of course the mailing list archives for those projects as well. The idea of this CSE is to help the user find articles, tutorials and documentation about using open source math tools to do specific tasks.

The results are pretty good. I choose the option of including broader web search results, so high ranking results from Wikipedia or MathWorld still show up, but the first page is full of relevant links to the favoured sites. For instance the results for 'bayes model averaging' have a mix of about half and half papers from academic websites and links to various R project pages, or the results for the Maxima function augcoefmatrix() which links to the Maxima manual as well as several web sites with usage examples.

It really improves the signal to noise ratio. Leave a comment if there's a site you think I'm missing.

Performance optimization (allocation inside a for loop)

An interesting discussion about an old topic just popped up on the Octave lists. It's not interesting because someone noticed that failing to pre-allocate a vector before entering a loop is slow (doesn't anyone RTFM?), but because of the follow-up discussion on indexing and range objects.

Here's the interesting bit:

octave:1> tic(); n=1e5; retval=1:n; toc()
Elapsed time is 0.000528962 seconds.
octave:2> tic(); n=1e5; retval = (1:n)(1:n); toc
Elapsed time is 0.00593709 seconds.
octave:3> tic();n=1e5;retval=[1:n]; toc
Elapsed time is 0.010952 seconds.

Why the significant difference in performance? According to jwe:

In Octave, an expression like 1:n creates a range object, which contains only the base, limit, and increment as double precision values, so no matter how many elements are in the range, it only takes a few bytes of storage (24 for the data plus some overhead for the internal octave_value object itself).

If you write [1:n], you force a Matrix object with N elements to be created. It will require 8*N bytes of storage, plus the overhead for the internal octave_value object itself.

Similar to the difference between range() and xrange() in Python.

Vectorize Matrix Slice Extraction

There was an interesting thread on the Octave mailing lists recently about vectorizing matrix slice extraction. That is avoiding the dreaded for-loop in something like:

%% create a random sparse matrix
n = 10000; %% matrix dimension
nz = 300; %% number of non-zero elements

ia = floor(rand(nz,1)*n); %% dimension 1 indices
ja = floor(rand(nz,1)*n); %% dimension 2 indices
A(:,:) = sparse(ia, ja, rand(nz,1), n, n, n);

%% taking slices with a for loop, slow, poor scaling to large matrices:
for jj = nz:-1:1
slow_columns(:,jj) = A(:,ja(jj));
endfor

%% example using built-in sub2ind, fast, good scaling to large matrices:
jj = sub2ind(size(A)(1),ja);
fast_columns = A(:,jj);

This is a simplified example, the one in the thread shows how you can access a 5-dimensional matrix with 4 indexes. Knowing the built-in functions that come with your high-level language is well worth it, even if it means you have to do weird things to stay in the high-level language and still get reasonable speed. This allows you to avoid slow looping and use the quick compiled functions that someone else has already spent time and effort developing.

Unblurring Images

Tikhonov regularization is a pretty cool concept when applied to image processing. It lets you 'de-blur' an image, which is an ill-posed problem. The results won't be as spectacular as CSI, but still cool nonetheless.

Original

The original image, CSI fans pretend it is a license plate or a criminal's face in a low quality surveillance video, was made by taking a screenshot of some words in a word processor. The image doesn't have to be character based, but characters make good practice problems.

Blurring

We then blur the image with a Toeplitz matrix. In Octave that's pretty easy:

im = imread("foo-fiddle.png");
[m, n] = size(im);

sigma_true = 10;
p = exp(-(0:n-1).^2 / (2 * sigma_true^2))';
p = p + [0; p(n:-1:2)];
p = p / sum(p);

G_true = toeplitz(p);

im = G_true*double(im)*G_true';

Generally after doing some floating point operations on an image you won't end up with something still scaled between 0 and 255, so the little function shown below rescales a vector's values to the proper range for a grayscale image.

function im = image_rescale(double_vector)
%% scale values linearly to lie between 0 and 255
m = max(double_vector);
n = min(double_vector);
scale = max(double_vector) - min(double_vector);
im = double_vector - (min(double_vector) + scale/2);
im = im ./ scale;
im = (255 .* im) + 127.5;
endfunction

De-Blurring

This is where the regularization comes in. Since the problem is poorly conditioned we can't just recover the image by straightforward inversion. We'll use the singular value decomposition of our blurring matrix to create a regularized inverse, this is similar to the pseudoinverse. If the singular value decomposition of the blurring matrix is


Then the regularized inverse is


Where the entries in the diagonal matrix are given by


So as alpha approaches zero we approach the un-regularized problem. Chapter 4 of Watkins has good coverage of doing SVDs in Matlab. The Octave to perform this process is shown below in the function tik_reg().

function im_tik = tik_reg(im,sigma,alpha)
[m,n] = size(im);
p = exp(-(0:n-1).^2 / (2 * sigma^2))';
p = p + [0; p(n:-1:2)];
p = p / sum(p);

G = toeplitz(p);

[U, S, V] = svd(G,0);

Sinv = S / (S^2 + alpha*eye(length(S),length(S))^2);

im_tik = V*Sinv*U'*double(im)*V*Sinv*U';
endfunction

This gives us two parameters (sigma and alpha) we have to guess in a real image processing problem (because we won't be the ones blurring the original). Another weakness of the method is that we assumed we knew the model that produced the blur, if the real physical process that caused the blur is not similar to Gaussian noise, then we won't be able to extract the original image.

Here's a couple for-loops in Octave that try several combinations of alpha and sigma.

sigma = [9,10,11];
alpha = [1e-6,1e-7,1e-8];

for j = 1:length(sigma)
for i = 1:length(alpha)
unblur = tik_reg(im, sigma(j), alpha(i));
unblur = reshape(image_rescale(reshape(unblur,m*n,1)),m,n);
fname = sprintf("unblur_%d%d.png",i,j);
imwrite(fname, uint8(unblur));
endfor
endfor

Optical character recognition probably wouldn't work on the output, but a person can probably decipher what the original image said.

Kronecker Product of Sparse Matrices

Gilbert Strang discusses the use of a neat built-in function in Octave, kron(), that's pretty useful for doing multidimensional finite difference methods on PDEs:


kron() takes the Kroenecker product of two matrices. This is especially useful when constructing the large, sparse matrices needed for finite difference approximations.

K is the tridiagonal matrix created for a central difference approximation to the second derivative:

n = 5;
K = spdiag(2*ones(n,1),0) + spdiag(-ones(n-1,1),-1) + spdiag(-ones(n-1,1),1);
spy(K)

K is the differences for a 1D row of points, now we can create the 2-dimensional version for a plane of points with kron():

K2D = kron(speye(n),K) + kron(K,speye(n));
spy(K2D);

The extension to 3 dimensions is exactly as before, now we have a matrix describing the differences between points on the same line in the first sub and superdiagonal, differences across lines in the next bands out from those, and differences across planes in the furthest out bands:

K3D = kron(speye(n),K2D) + kron(K2D,speye(n));
spy(K3D);

Pretty simple! Now you're ready to solve multidimensional partial differential equations.

Fibonacci Numbers

Gilbert Strang has an interesting example in one of his linear algebra lectures on eigenvalues that is relevant to the articles I've written on speeding up calculations in Octave.

You might be familiar with the Fibonacci numbers, they are calculated by a two term recurrence:


We've already seen that recasting a calculation in terms of a matrix vector multiply or matrix vector solve is an important skill for writing good fast code in Octave. Strang demonstrates this at about halfway through the above lecture. The first order system that describes the generation of Fibonacci numbers is:


The nth term of the sequence is calculated by taking the nth power of A and applying it to the initial conditions:


When you diagonalize a matrix the powers of the matrix become very easy to calculate:

If you aren't familiar with the concept of eigenvalues and eigenvectors and diagonalizing a matrix I would highly recommend the linear algebra lectures on MIT's OCW website.

Here's the Octave code to accomplish the above. If you want, you can use tic() and toc() for timing. On my laptop the matrix method is about three times faster than the for loop method. What kind of results do you get? This advantage will grow as you go further into the sequence because the amount of work stays constant for the matrix method, but grows linearly for the recurrence method.

%% Our 'fibonacci' operator:
A = [1,1;1,0]
%% our initial conditions:
u_0 = [0;1];

%% get the matrix of eignevectors, V, and the matrix of eigenvalues, lambda
[V,lambda] = eig(A)

%% A = V lambda V^{-1}

Vinv = inverse(V)

%% This gives us an idea of the error introduced by calculating the
%% inverse with finite precision arithmetic:
Error = A - V*lambda*Vinv

n = 1400;
F = zeros(n,2);
%% Find the 100th Fibonaci number:
F(100,:) = V * (lambda^100) * Vinv * u_0;

%% check our answer with a simple loop:
fibo = zeros(n,1);
fibo(1) = 0;
fibo(2) = 1;
for i=3:n
fibo(i) = fibo(i-1) + fibo(i-2);
endfor

for i=2:n
F(i,:) = V * (lambda^i) * Vinv *u_0;
endfor

relative_error = abs( ( F(3:n,2) - fibo(3:n) ) ./ fibo(3:n) );

One thing to note is that one of the eigenvalues of our matrix is greater than one. This has implications for the numerical stability of the above method for calculating terms in the Fibonacci sequence. It means that the initial tiny round-off error in our eigenvector matrix inverse will be amplified as we continue to calculate higher and higher terms in the sequence. This behaviour is shown in the plot below, luckily the error starts out very small, and doesn't grow too large before we run out of range in the 32bit floating point numbers on the machine.

What's the answer? The 100th Fibonacci number is approximately 2.1892e20.

How to Speed Up Octave

Speed seems like a continuing theme on the Octave mailing lists, a recent thread on the list of "comparing the executive speed with Matlab" demonstrates that there really is very little speed overhead when you use Octave smartly. If you read down far enough, it also demonstrates that the Intel Fortran compiler is awesome!

There are four basic methods of speeding up a chunk of Octave code (in order of increasing speed):

1. Allocate memory smartly (using for loops)


2. Vectorize your code (avoiding for loops)
   


3. Use sparse matrix operators in place of for loops you can't vectorize


4. Write a compiled function in C/C++/Fortran to call from Octave

The first three options keep you in pure Octave. The fourth option, while requiring you to use a compiled language is really not too bad because you can use all of the types (vectors, matrices, etc) defined in the Octave libraries to write your function. That's almost as easy as doing it in Octave.

Speeding Up For Loops

This one really amounts to not being foolish. The only way to speed up for loops is to make sure you allocate all of your memory (using zeros() for instance) before you enter the for loop. The Octave Manual covers this topic here. The better thing to do is vectorize your code and avoid for loops all together.

Vectorizing

What you really want to do is cast your calculation as a series of vector and matrix operations, or use the built-in (usually compiled) functions that operate directly on vectors. The post Is Octave Slow? demonstrates the sorts of speed up you can get when you get rid of for loops.

Sparse Matrices for Data Dependancy

What about for loops you can't vectorize? The ones where the result at iteration n depends on what you calculated in iteration n-1, n-2, etc. You can define a sparse matrix that represents your for loop and then do a sparse matrix solve. This is much faster because you are now using compiled functions. The drawback is that you are trading some increase (sometimes substantial) in memory usage to buy this speed up.

• Is Octave Slow? covers this approach as well.


• More on Speeding up Octave covers the details on sparse matrix allocation.

Compile It

The manual has a good intro on using mkoctfile to compile your function and make it callable from Octave. Probably the most useful thing is that you can use the Octave types, so you don't have to code up your own matrix or vector type. At this point the speed really has little to do with Octave and everything to do with the quality of the compiler you use and how efficiently you code up your calculation.

Conclusion

The real power of Octave isn't that you can get your calculations to be nearly as fast as if you wrote them all in a compiled language. It's that you can prototype a calculation quickly and make sure that it is correct, and then there is a fairly well defined path to follow to speed up your correct calculation.

More on Speeding Up Octave

In a previous post we found that speeding up Octave amounted mostly to avoiding for loops. The simplest way to do that is to operate directly on vectors using the built-in operators and functions (which are fast because they are compiled) or cast your problem as a sparse matrix solve. This second option is especially helpful when subsequent calculations depend on previous iterations (the code can't really be vectorized). This option is a very general way of avoiding loops in those cases.

The simple example given previously was a numerical integration of x-squared:

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

The way to recast the problem as a sparse matrix solve is to think of p1 as the vector of unknowns, and each iteration of the loop as an equation in the system we want to solve. Writing down the system gives us:


The important detail to remember is to use the functions in Octave to allocate the sparse matrix, or you could easily find yourself writing more really slow for-loops just to create the sparse matrix which you hoped would save lots of time by avoiding for-loops. Talk about the long, slow way around!

Two very useful functions are speye() and spdiag(). The first returns a sparse identity matrix, which is often a good initial building block for many useful operators. The second allows you to place vectors (allocated quickly with the usual vector suspects such as ones() and zeros() ) along the diagonals of the sparse matrix.

# create the main diagonal
A = speye( N );
# alternatively could use spdiag:
# A = spdiag( ones(N,1), 0 );
A = A + spdiag( -ones(N-1,1), -1 ); # add the first sub-diagonal

If your operator isn't banded then you'll need to use spconvert(), which takes as its argument a three (or four if you need a complex result) column matrix. Each row of the argument defines a non-zero entry in the sparse matrix. The first column is the row index, the second column is the column index and the third column is the entry value (fourth column being the imaginary part of the value, if necessary).

Using those three functions should allow you to allocate a sparse matrix in Octave without resorting to for loops (which was why we embarked on this journey to begin with).

This method is counter-intuitive to folks who come from a Fortran (or other compiled language) background, because writing down the loop is the simple, efficient way to solve the problem. It also seems like 'wasting' memory to store all of those redundant coefficients. The timing results speak for themselves though, if you want to stay completely in Octave (or Matlab or Python) sometimes you have to do weird things to get reasonable performance. Of course, if speed really becomes a problem then the inner loops of your calculations need to move to a compiled language that can then be called from Octave, this is a bit more complicated than our little sparse matrix method.

LED Light Board

The first part of my aeroponics project (overview) is turning those blue and red LED Christmas lights into a grow lamp.

The basic design is a 2 dimensional array of alternating red and blue lights mounted to a piece of white foam board.

Here's a plot of the layout (generated in Octave of course):

And here is the light board as built:

You'll notice two empty holes on the right hand side of the board, the wire between the lights was a really tight fit to meet the spacing I chose (should have checked that a little more closely before drilling all those holes).

Those mini-LEDs look pretty bright, but I'm not sure if it is producing enough light to actually grow a plant, might need to add a string or two. The packaging claims that you can safely plug 43 (!!) of the 60 light strings into a single outlet.

Poor Man's ADC

Can you use a commodity webcam as an el-cheapo analog to digital converter (ADC)? The folks at slashdot want to know so they can reduce their energy consumption.
Here's the example image:
Octave has plenty of image processing functions in the image package from Octave Forge. Here's the edges detected with the built-in edge() function:

Now we need to perform a smart correlation with the gauge needles so we can map the angle of the needle to a decimal value in [0,10). One way to do it would be to zoom in on the center of the gauge:

Then we can perform a simple least squares regression on the non-zero pixel locations to get a rough estimate of the "slope" of the needle.

needleI = pngread( "needle.png" );
[i,j] = find(needleI(1:35,1:35));
y = 40-j;
x = i;
X = [ ones(length(x),1), x ];
A = X'*X;
b = (X')*y;
beta = A\b;

For the first needle pictured above we get a slope of -0.19643. This seems a tad low, so we probably need to correct for the slant in the original image.

Octave script with all the calculations demonstrated above.

Also check out this book or this page for lots of image processing examples in MATLAB/Octave.

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!

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)