Spatial Cross-correlation Function

Question

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?

Answer 1

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