Spatial Cross-correlation Function

I have code that calculates a bias-adjusted estimator xi(r) defined as a measure of the excess probability dP, above what is expected for an unclustered random Poisson distribution, of finding a galaxy in a volume element dV at a separation r from another galaxy.

The inputs are two spatial point sets:

1. The galaxy data catalogue (size $n$)
2. The random point catalogue (size $n_r$)

The Landy and Szalay estimator is

where $DD$ is the count of data-data pairs (within distance $r$), $RR$ is the count random-random pairs (within distance $r$), and $DR$ is the count of data-random pairs (within distance $r$). All counts are suitably normalised.

Q1. What are data-random pairs?

Q2. How does one calculate DR, the count of (data-random) pairs?

1 Answer

The estimator you are referring to comes from Bias and Variance of Angular Correlation Functions.

$D$ is an empirical sample of galaxies, typically captured as a CCD image.

$R$ is a simulated point distribution with the same mean density and sampling geometry as $D$.

$DD$ are the number of galaxy pairs (within radius r) from distribution $D$.

$RR$ are the number of point pairs (within radius r) from distribution $R$.

$DR$ are the number of point-galaxy pairs (within radius r) from the joint distribution of $R$ and $D$.

Using the notation of the paper, there are $n$ points in $D$, $n_r$ points in $R$, and $n$$n_r$ point-galaxy pairs.

So for $i \in \{1,...,n_r \}$ take $R_i$, for $j \in \{1,...,n\}$ count $D_{ij}$ if it lies within radius r, sum these counts $\forall i,j$, normalise using $nn_r$.

Example

Conceptually, using euclidean 2D space with a square geometry and using $n = n_r$:

set.seed(1)

n  <- 96
nr <- n
k  <- 8
r  <- 5

D <- data.frame(x = unlist(lapply(1:k, function(i) rnorm(n/k, runif(1)*i^2))),
y = unlist(lapply(1:k, function(i) rnorm(n/k, runif(1)*i^2))))

R <- data.frame(x = runif(nrow(D), min(D$x), max(D$x)),
y = runif(nrow(D), min(D$y), max(D$y)))

plot(D, col='red', main='Simulated Joint point-galaxy Distribution')
points(R, col='blue')


## normalised counts
DD <- sum(dist(D)<r) / (n  * (n  - 1) / 2 )
RR <- sum(dist(R)<r) / (nr * (nr - 1) / 2 )
DR <- 0
for (i in seq(nr))
for (j in seq(n))
DR <- DR + ifelse(sqrt((R$x[i]-D$x[j])^2 + (R$y[i]^2-D$y[j])^2) < r,1,0)
DR <- DR / (n * nr)

(xi <- (DD - 2 * DR + RR) / RR)

[1] 2.98162

• I am not a statistics student, so as far I can understand, you are saying. Take a radius r, position at a random point in the data and count the number of galaxy pairs within that radius which is DD. Do the same for the random data and the pairs counted are RR. Now I count the total number of pairs of both data and randoms within the same fixed radius, which will be my DR??? – ThePredator Nov 29 '14 at 8:59
• You need to be careful about how you phrase it but you're on the right track. For DR, you want to take each point in R, one at a time, and calculate distances to every point in D. You will produce $nn_r$ distance pairs. The subset of these that are shorter than r is the set DR. – goangit Nov 29 '14 at 10:39
• So to visualise, can I say this. I overplot data points (D) and random points (R), then as you said, I take one random point Ri and pair them with all the data points and find out which pairs nnr are within this distance r. So all these pairs nnr that are within this distance is DR. Am I correct? Then I do so for the remaining randoms points – ThePredator Nov 29 '14 at 11:16
• So DR is the number of subsets of these pairs within r? – ThePredator Nov 29 '14 at 11:19
• Frankly, I'm running out of ways to say the same thing. I'm pretty sure that if you try and program this you will discover what is meant by the above. Find all $n$ x $n_r$ distance pairs from $D$ and $R$. Count the ones shorter than r. That's DR. – goangit Nov 29 '14 at 11:44