# Measuring error in the lag when cross-correlating data

Cross-correlation can be used to measure a 'lag' (or shift or offset) between two data sets (e.g., data streams). Is there a standard way to measure the errors in the 'lag'? Or better yet, is there a good way to measure those errors?

Some more details: I'm performing 2D cross-correlations between images in order to achieve sub-pixel image registration. I would like to know how accurate the measured offsets are (i.e., ARE they sub-pixel?). These are images of astrophysical sources, but not stars - all of the emission is extended. For images with stars, you can measure the centroids of the stars to determine offsets between images: that is not possible for this data.

The Gaussian approach suggested by Nestor is one that I've attempted, but I'm not entirely clear on some points, particularly how to account for individual pixel errors. I suspected that other fields might use this same technique but call it by another name, which is why I asked this question here.

-
What do you mean by error in lag. If you do crosscorrelation the highest lagged correlation might be taken as a measure of how one time series lags behind the other. But what are you looking for in this? Is it that you want to address the situation when two or more adjacent lags have high crosscorrelation values but the highest is only highest by chance? –  Michael Chernick Aug 15 '12 at 16:00
@MichaelChernick - yes, essentially. Given that there is measurement error in each element within a data stream, there is some uncertainty on the highest value in the cross-correlation function [using terminology from Nestor's post) –  keflavich Aug 15 '12 at 20:15

Edit:

Now that you edited your question, I should add that the procedure of Tonry & Davis that I describe below is pretty general. In fact, appart from reducing the squared difference between the convoluted, shifted and scaled template and the observed spectrum (which can be thought of as a Maximum Likelihood estimator of the parameters with equal errors, which can be accounted for [see my comments]) the only part where they make strong assumptions is on the shape of the convolution of the template and on the shape of the largest peak in the CCF. You described that your objects may have different shapes; have you tried actually plotting the shape of the CCF? Maybe an extended gaussian or a Moffat function are more suitable choices in your case in order to model the 2D CCF peak. Analytically tracking the shape of the CCF, I think, is the easier path; this would allow you to analytically track the errors on the CCFs.

At least in astronomy we "have" a standard way of doing this, and it's the method of Tonry & Davis (1979). The part that may interest you in that paper is Section III, but let me explain to you a little of the background of their method in order to have an easier reading :-).

The idea in astronomy is that you have a source that emits certain flux (energy per meter$^2$ per second) at different wavelengths, say, $f_{\text{source}}(\lambda)$: this is called the spectrum of the source. However, because the objects that Tonry & Davis measure (galaxies, by the way) are receding from our point of view in our galaxy, all the known features of this object (say, bumps at certain wavelengths because of atomic absorption or emission) are generally shifted towards the red (i.e., we see them at longer wavelengths. For example, if we expected to see a bump at $\lambda=4500$, we may actually see it at $\lambda_{\text{source}}=4510$). Their work, then, is focused on calculating this wavelength shift and, of course, measuring the error on this shift.

Note that in their work, Tonry & Davis make a conversion between wavelength and bins (because we bin this flux as a function in wavelength in pixels; both to have higher signal to noise ratio and because a CCD camera is the best instrument to date for measuring flux) and from there they measure this wavelength shift. Maybe in your case you don't actually need this conversion between the bins Tonry & Davis talk about and wavelength, so you may want to change equation (1) in their paper to suit your needs. Another important feature of this paper is that they weight any deviations from zero quadratically, because these "bumps" I was talking about in the first paragraph of this answer are clearly more important. However, you can modify their $\chi^2$ reduction scheme to suit your needs ;-).

Finally, note that "all they do" is to approximate the largest peak in the Cross-Correlation Function (CCF) by a gaussian, and from there do the error analysis. In practice, I've seen this approximation to work pretty well, but always check the shape of your CCF, just in case.

PS: I didn't post this answer on that post you cited because JBWhitmore appeared to work on astronomy (or be an astronomer), and almost every astronomer I know knows the paper of Tonry & Davis, so I thought he was searching for something else.

-
I'm a Galactic astronomer, so I've never had to do spectral cross-correlation =). Anyway, while my question is closely related to spectral cross-correlation, I've found that the data I'm working with have CCFs that don't look very gaussian... I'll explain more in an edit to the original question –  keflavich Aug 15 '12 at 19:34
@Nestor What do you mean by approximating the CCF by a Gaussian? Are you assuming zero cross correlation with the non zero outcome just due to white noise or do you mean a nonzero cross correlation function with gaussian error about the true crosscorrelation? –  Michael Chernick Aug 15 '12 at 20:27
@keflavich, the individual pixel errors can be accounted for directly on the $\chi^2$ reduction. You can think in Equation (7) in the paper of Tonry & Davis as a search of the maximum likelihood estimator with unit variance, so you can rephrase it as: $$\chi^2=\sum_{n}\frac{1}{\sigma_n^2}[\alpha t*b(n-\delta)-g(n)]^2,$$ where $\sigma_n^2$ is the variance of $g(n)$ at bin $n$ and where you are assuming gaussian uncorrelated errors (this is ok as long as you have enough photon flux. Of course, if the flux is low you should use the Poissonian Likelihood and from there calculate the MLE). –  Néstor Aug 16 '12 at 1:13
@MichaelChernick, I meant to say that the largest peak on the CCF is approximated by a gaussian. I'll edit that on my original post (sorry, just used to the astronomical jargon...). –  Néstor Aug 16 '12 at 1:16