Sigma clip by axis (addressing #1037)

I liked suggestion #1037 by @lisafogarty of allowing one to `sigma_clip` along a given axis. The implementation was not completely straightforward since as written the procedure assumed 1-D arrays throughout. My alternative uses masked arrays; it works well, but has the slight disadvantage that user functions for determining a center and variance now have to be able deal with masked arrays. Nevertheless, it is neat to be able to do:
```
from astropy.stats import sigma_clip
from numpy.random import normal
from numpy import arange, diag, ones
data = arange(5)+normal(0.,0.05,(5,5))+diag(ones(5))
masked = sigma_clip(data, axis=0, sig=2.3)
```
And note that `sigma_clip(data, axis=1, sig=2.3)` will, correctly, not mask anything.

Note also that for clipping along an axis I made `maout=True` the default (i.e., return a masked array, since for that case a 1-d flattened array makes no sense; indeed, I raise an exception on `maout=False`). I wonder more generally about the use case of returning just the good points; for most astronomical applications, it would seem much more useful to keep the shape and thus get the masked array. Personally, I'd remove the `maout` option and replace it with a `copy`.