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

enter image description here

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?


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$.


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


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')

enter image description here

## 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
| cite | improve this answer | |
  • $\begingroup$ 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??? $\endgroup$ – Srivatsan Nov 29 '14 at 8:59
  • $\begingroup$ 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. $\endgroup$ – goangit Nov 29 '14 at 10:39
  • $\begingroup$ 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 $\endgroup$ – Srivatsan Nov 29 '14 at 11:16
  • $\begingroup$ So DR is the number of subsets of these pairs within r? $\endgroup$ – Srivatsan Nov 29 '14 at 11:19
  • $\begingroup$ 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. $\endgroup$ – goangit Nov 29 '14 at 11:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.