WIP/ENH: add repeated measures twoway anova function

This addresses #226 and some aspects of #535

Adapted from my gist: 
https://gist.github.com/dengemann/5427106

which is a translation of MATLAB code by Rik Henson:
http://www.mrc-cbu.cam.ac.uk/people/rik.henson/personal/repanova.m

and to some lesser extend by Python code from pvttbl by Roger Lew.
http://code.google.com/p/pyvttbl/

While there is a new WIP PR in statsmodels related to this (https://github.com/statsmodels/statsmodels/pull/786), this minimal version aims at supporting our (mass-univariate) use case.

Some features:
- supports joblib (on my 2011 macbook air it took me 7 minutes to
compute 1.000.000 repeated measures anovas for 18 subjects and
all three effects from a 2 x 2 design using 2 jobs). I might find a more efficient way to do this, but I think this is definitely start.
- supports sphericity correction for factor levels > 2

Current limitations are:
- to keep things simple, I constrained this function to only estimate 2 factors. This should make sense in an MEG context.
- both factors are expected  to repeated, so no 'between-subject' effects.

Test and examples are coming, so please wait with strict reviewing. 
For the API, however I could already need some feedback. I was wondering whether to either add some iterator support for retrieving data matrices for the effects desired sequentially (assuming you don't want all three 20484 * e.g. 150 matrices at once when applied in a source space analysis). Another option would be to add some mini formula language allowing to request one effect e.g. using some R conventions: 'A'  or 'B'  for the main effects,  'A:B'  for just the interaction or  'A*B' for the full model. I think both iteration and formula will be too much.