math: anatomy of a multivariate Gaussian

these are really good lecture notes, explaining the linear algebra, covariance matrices and so on. Directly helps for coding.

stats: the covariance matrix

The covariance matrix is N-dimensional generalisation of the scalar variance in 1 dimension. Here's a simple explanation of it which I needed in order to code up a test case for fitting a 2D Gaussian with intrinsic scatter.

Python: bootstrapping with sklearn

sklearn.cross_validation.Bootstrap returns indices of random bootstrap samples from your data. I am not yet well familiar with the use of bootstrapping for 2 sample comparison, so I'm using means as a way to characterise the two distributions (that's better than 2D K-S test, which doesn't exist anyway!). This is of course an incomplete statistic (I have a reason to suspect that two of the samples I am working with have different skewnesses). I'm working with three sets of absolute r magnitudes (M_r) here.
Here's how it looks in practice:

from sklearn import cross_validation

#let's call our sample array data. It can be N-dimensional.

#len(data) -- total number of data points in my dataset, 
#nBoot -- number of bootstrap samples, 
#train_size = bootstrap sample size (proportion of the whole sample, or just number)
#Here we create an empty array we will store the means of bootstrapped 

means = np.empty((nBoot, ))

#then we get the class instance, drawing half of the sample with replacement every time

bs = cross_validation.Bootstrap(len(data), nBoot, train_size=0.5, random_state=0)

#Filling the means array while iterating over the bootstrap samples (indexed by train_index):

i = 0
for train_index, test_index in bsJ:
    means[i] = np.mean(data[train_index])
    i+=1

I've repeated it for all three distributions I'm interested in, and here is a histogram of all the bootstrap sample means. The difference is just what what I've expected: basically, some faint galaxies (with M_r < 20 or so) were manually rejected from the green and red distributions, so the mean of them is shifted towards brighter absolute magnitudes. It remains to be seen how important this is. If we think this is worth it, I'll try using other distributions' statistics some other time.

Python: scikit-learn MoG

Probably should have tried this, before reinventing the wheel.

MCMC -- how to sample a multi-modal Gaussian

see this.

lit: Some aspects of measurement error in linear regression of astronomical data, Kelly 2007

Kelly 2007:
I describe a Bayesian method to account for measurement errors in linear regression of astronomical data. The method allows for heteroscedastic and possibly correlated measurement errors and intrinsic scatter in the regression relationship. The method is based on deriving a likelihood function for the measured data, and I focus on the case when the intrinsic distribution of the independent variables can be approximated using a mixture of Gaussian functions. I generalize the method to incorporate multiple independent variables, non-detections, and selection effects (e.g., Malmquist bias).

Daily paper: The Tully-Fisher relations of early-type spiral and S0 galaxies

I gave a talk at our group meeting yesterday, this was one of several useful suggestions and references I got. bWilliams et al. 2010, The Tully-Fisher relations of early-type spiral and S0 galaxies
Offsets which are ascribed to fading or brightening could in fact be due to systematic biases that differ between the measures of rotation used. For typical slopes of the TFR, a relatively small systematic difference in velocity of ∼0.1 dex introduces an offset in the TFR that is indistinguishable from a ∼ 1mag difference in luminosity.
Our conclusions in this section have three consequences for previous works that measure the TFR of S0 galaxies using stellar kinematics corrected for asymmetric drift (e.g. Neistein et al. 1999; Hinz et al. 2001, 2003, BAM06): (i) the asymmetric drift correction they use yields similar results to detailed Jeans modelling; (ii) the drift correction does not seem to introduce a systematic bias between spirals and S0s and, if applied to samples of both spirals and S0s, can indeed be used to compare the TFRs of the two classes; (iii) however, a TFR derived from asymmetric drift-corrected stellar kinematics cannot be directly compared to TFRs derived from global HI line widths or resolved emission line PVDs. ==== There's a really interesting section on the ADC, and a nice comparison of gas/HI/AD corrected stellar kinematics. Also I should check their section on TFR fitting, intrinsic scatter in it, etc.In addition, I ought to take a look at the references therein, esp. the introduction.

The difference between MCMC and simulations

from PSU astrostatistics lecture notes:
MCMC is a form of simulation, because MCMC generates random numbers and bases the inference on them. However, there are two important differences between MCMC and the simulation that we shall discuss here.
a) In MCMC we do not simulate from the target distribution directly. Instead, we simulate a Markov Chain that converges to the target distribution. Thus, MCMC performs approximate simulation.
b) The primary use of MCMC is in Bayesian statistics where we seek to simulate from the posterior distribution of the parameters. Classical simulation on the other hand simulates fresh data.

Uncertainties are distributions, not errorbar values

This clicked on me while preparing for a talk, though it's an obvious thing: too often we are content with just plotting errorbars and not specifying what do they mean (for instance: 1 sigma-uncertainties, Gaussian distribution). This makes propagating uncertainties (and using the data for model selection) actually impossible.
To put it another way: we don't have datapoints, we operate with slippery, diffuse clouds of probabilities of data. My current task is to take full 2D uncertainties distribution for my datapoints and re-use them in the following steps of hierarchical analysis.

lit: Probability distributions of ellipticity

My pick for the weekly meeting: On the Probability Distributions of Ellipticity, by Viola et al. A review from the weak lensing perspective.

lit: Precise TFRs without galaxy inclinations

Precise TFRs without galaxy inclinations, by Obreschkow et al. A really interesting read about why photometric (axis ratio-based) inclination values are not optimal, what TFR is good for, fitting the TFR, etc. The main point of the article is recovering the TFR for samples with limited or no inclination data, but I haven't read it in detail yet.

Fitting the Tully-Fisher relation

There's one thing I can't quite wrap my head around while reading literature on the Tully-Fisher relation: why is the fitting always performed in the magnitude-log(velocity) space? It may be easier to fit an expression y = mx + b to data than y = c*x**z, but one would have to deal with asymmetric errors in two dimensions.
I'm currently trying to figure out a way to do the fit in semi-log space, namely, magnitude-velocity space, using errors in both directions and MCMC.

Math: Directional statistics, calculating average of angles in Python

How does one calculate the average of a distribution of angles? The mean of 360 deg and 0 deg is not 180 degrees.
I was quite certain that while it is possible to get around the problem programatically, there should be a host of mathematical methods to deal with functions that have circular domains. That's directional statistics and it deals with statistics in certain non-Euclidean spaces.
In a nutshell, that's how I calculated the average of angles in Python (and in radians!):

def getCircularMean(angles):
  n = len(angles)
  sineMean = np.divide(np.sum(np.sin(angles)), n)
  cosineMean = np.divide(np.sum(np.cos(angles)), n)
  vectorMean = math.atan2(sineMean, cosineMean)
  return vectorMean




Python: An Algorithm for Precise Aperture Photometry of Critically Sampled Images

An Algorithm for Precise Aperture Photometry of Critically Sampled Images, Bickerton and Lupton 2013.

Python: splitting a list into a number of roughly equal chunks

I used this solution, however, since I wanted to have a set number of chunks as opposed to chunks of a set length, the function ended up being slightly different:

def splitList(longList, chunkNumber):
    #Get the number of elements in a chunk we want
    n = int(math.floor(len(longList)/chunkNumber)) + 1
    print n, 'length of a chunk'
    return [longList[i:i+n] for i in range(0, len(longList), n)]

Daily Paper #10: Reconstructed density and velocity fields from the 2MASS Redshift Survey

Today's article was referenced to in the WISE TF paper I summarised yesterday. I wanted to find out how the peculiar velocities and density fields were reconstructed.
Title: Reconstructed density and velocity fields from the 2MASS Redshift Survey
Authors: Erdoğdu, Lahav et al.
Year: 2006, http://mnras.oxfordjournals.org/content/373/1/45.full
I have to admit I did not understand the mathematics of this article, and will read Peebles 1973 and Fisher 1995 to understand what are the spherical harmonics and Fourier-Bessel functions in detail. What they do in order to transform the the redshift space to the real space is decompose the redshift space distribution into radial and angular components. Peculiar velocities affect only the radial component, so it is easier to deconvolve them.
The relationship between mass and light distribution is usually assumed to be linear with some proportionality constant (there is an article by Vivienne arguing that this is not the case).
They provide density and velocity maps up to cz = 16 0000 km/s, which covers the CALIFA redshift range. I could ask them for the velocity data and obtain the peculiar velocities/infall corrected distances for the sample. I think it would be useful for me later, when working with the density fields and spin alignments. We would also have these distances and their errors for all the galaxies, not only the ones identified as belonging to structures.
Previously I used the Virgo-GA-Shapley infall corrected redshifts from NED (obtained using this model: http://iopscience.iop.org/0004-637X/529/2/786/fulltext/). It's a simpler model, taking only these three attractors into account.

lit: Asymmetry of galaxies

by Conselice, Bershady et al.

lit: 2dSWML algorithm, HI mass function

The 1000 Brightest HIPASS Galaxies by Zwaan et al.: has a nice section on Vmax corrections/fitting distribution functions, with a reference to the original Schmidt 1968 article.

lit: Spins and orientations of haloes in cosmic filaments

Zhang et al., many useful references on methods therein.

Fitting a model to data with errors in x, y directions

Here's a quick implementation to try.