11institutetext: Ludwig–Maximilians Universtät München, Fakultät für Physik,
Schellingstr. 4, D-80799 München
11email<EMAIL_ADDRESS>
# Improving the accuracy of estimators for the two-point correlation function

Martin Kerscher

(March 24, 2022)

###### Abstract

Aims. We show how to increase the accuracy of estimates of the two-point
correlation function without sacrificing efficiency.

Methods. We quantify the error of the pair-counts and of the Landy & Szalay
estimator by comparing with exact reference values. The standard method, using
random point sets, is compared to geometrically motivated estimators and
estimators using quasi Monte-Carlo integration.

Results. In the standard method the error scales proportional to
$1/\sqrt{{N_{\mathrm{r}}}}$, with ${N_{\mathrm{r}}}$ the number of random
points. In our improved methods the error is scaling almost proportional to
$1/{N_{\mathrm{q}}}$, where ${N_{\mathrm{q}}}$ is the number of points from a
low discrepancy sequence. In an example we achieve a speedup by a factor of
$10^{4}$ over the standard method, still keeping the same level of accuracy.
We also discuss how to apply these improved estimators to incompletely sampled
galaxy catalogues.

###### Key Words.:

Methods: statistical, data analysis, Cosmology: large-scale structure of
Universe

##  1 Introduction

Often a statistical approach is used to characterise the large-scale structure
of the galaxy distribution. Then we assume that the distribution of galaxies
is a realisation of a point process (see e.g. Neyman & Scott 1958; Peebles
1980). Observations give us the positions of galaxies in space. From this set
of points we estimate the moments of the point process, specifically the two-
point correlation function $\xi(r)$ or the two-point density

$\displaystyle\varrho_{2}(\bm{x},\bm{y})=\varrho^{2}\
\left(1+\xi(|\bm{x}-\bm{y}|)\right),$ (1)

the probability of finding two galaxies at $\bm{x}$ and $\bm{y}$, where
$\varrho$ is the number density. In a homogeneous and isotropic point process
$\xi(r)$ only depends on the separation $r=|\bm{x}-\bm{y}|$. To determine
$\xi(r)$ from a galaxy catalogue within a finite domain
$W\subset\mathbb{R}^{3}$ we use estimators. In cosmology mostly estimators
based on random point sets are used. They rely on the data-data DD, data-
random DR, and random-random RR pair-counts (see below for the definition).

The two-point correlation function of the galaxy distribution is often used to
constrain models of structure and galaxy formation and to estimate parameters
of cosmological models. In current and upcoming galaxy samples the positions
of millions up to billions of galaxies will be observed (Dawson et al. 2013,
BOSS, Ross et al. 2020, eBOSS, Abbott et al. 2022, DES, Aghamousa et al. 2016,
DESI, Alarcon et al. 2021, PAUS, Ivezić et al. 2019, LSST, Amendola et al.
2018, Euclid). Fast and reliable methods for calculating the two-point
correlation function are needed. These large samples allow us to reduce the
statistical error from cosmic variance. But for the error budget of the two-
point correlation function we also have to control the systematic errors. One
systematic contribution is the error from random sets used in the pair-counts
DR and RR. As an example consider the baryon acoustic oscillations (BAO) which
lead to a peak in the two-point correlation function of galaxies at a scale of
about 100$\text{Mpc}/h$ (Eisenstein et al., 2005; Bautista et al., 2021). This
BAO peak has a height of approximately 0.01 above zero (compare Fig. 2 in
Eisenstein et al. 2005). For a percent level accuracy we need to calculate the
two-point correlation function with an absolute accuracy of less than
$10^{-4}$. We will show how to reduce the systematic error to this level of
accuracy without sacrificing efficiency.

Several estimators for the two-point correlation function have been developed
(Peebles & Hauser, 1974; Hewett, 1982; Davis & Peebles, 1983; Rivolo, 1986;
Landy & Szalay, 1993; Hamilton, 1993). By comparing these estimators to a
reference result from a cosmological simulation Kerscher et al. (2000) found
that the Landy & Szalay (1993)-estimator is the preferred estimator with the
smallest deviation from the reference and also that its bias is negligible
compared to its variance.

We focus on methods for increasing the numerical accuracy of the pair-counts
as used in these estimators. The random point set, shared in the pair-counts
RR and DR is used to correct for boundary (finite-size) and inhomogeneous
sampling effects. We will see that RR and DR are Monte-Carlo volume
integration schemes. As expected from standard Monte-Carlo integration, the
error of these pair-counts RR and DR scales as $1/\sqrt{{N_{\mathrm{r}}}}$,
where ${N_{\mathrm{r}}}$ is the number of random points used. This slow
convergence rate makes increasing the accuracy costly, sometimes unfeasible.
To improve the accuracy of the pair-counts without sacrificing efficiency we
follow two directions:

  * •

The pair-counts can be expressed as averages of specific volume fractions. We
use this to propose special adapted volume integration schemes which, in turn,
can be calculated more efficiently than the standard approach.

  * •

We replace the standard Monte-Carlo scheme with a quasi Monte-Carlo
integration, which leads to an improved scaling of the error almost
proportional to $1/{N_{\mathrm{q}}}$, where ${N_{\mathrm{q}}}$ is the number
of points from a low discrepancy sequence.

We compare the standard and the new methods to exactly known reference values.
This allows us to empirically validate the asymptotic scaling of the errors.

A variety of approaches have been suggested to improve the speed and accuracy
of estimators for the two-point correlation function. Keihänen et al. (2019)
show that at fixed computational cost, a split random catalogue improves the
accuracy of estimators for the two-point correlation function. For galaxy
catalogues Demina et al. (2018) achieve a speedup by factorising the
calculations in radial (redshift) and angular coordinates (see also Breton &
de la Torre 2021). Perhaps closest to our work are the investigations by
Dávila-Kurbán et al. (2021). They use glass-like point sets instead of the
random catalogues, where we use low discrepancy sequences. As reference values
in our comparisons we use the exact results from Baddeley et al. (1993) and
Kerscher (1999) for a rectangular box as summarised in Appendix B (for
periodic boxes see Appendix A.1). Additional to these exact results He (2021)
also discusses approximations.

Other approaches focus on the computational problem of calculating the pair-
counts. Tree based methods can be significantly faster than a direct
implementation of the pair-counts, specifically for small radii (Moore et al.,
2001). The double loop in the pair-count calculations can be parallelised.
Alonso (2012) show how to obtain a speedup by a factor of 100 over the direct
implementation by using multi-threading on multi-core CPU’s or utilizing many
cores in GPU’s. It is well known from matrix computations, that the memory
layout of the data can have dramatic consequences on the running time of
algorithms (see e.g. Anderson et al. 1999). Also for the pair-counts a clever
layout of the coordinates in the memory can lead to a significant speedup
(Donoso, 2019). A similar approach can be combined with multi-threading and
vectorization resulting in a blazingly fast code (Sinha & Garrison, 2020). Our
conceptual improvements can be combined with these computational speed-ups.

In Sec. 2 we give the definition of the pair-counts and discuss the geometry
of the expected pair-counts. In Appendix A we give the details and in Appendix
B we recap some geometric results for simple sample geometries. Together with
Sec. 2 this enables us to calculate exact reference values for the pair
counts. At the end of Sec. 2 we give a short introduction to the quasi Monte
Carlo method, as used in our improved estimators. In Sec. 3 and 4 we compare
the standard and improved versions of the pair-counts with the exactly known
reference values. We discuss the scaling of the error with the number of
points. We put this together in Sec. 5 and show how an improved version of the
Landy & Szalay (1993)-estimator can be constructed and test its performance on
a simulated cluster catalogue. In Sec. 5 we discuss how these improvements
have to be adapted to estimate the two-point correlation function from an
inhomogeneously sampled galaxy distribution. In Appendix C we discuss details
of the implementation and give a link to the code.

##  2 Paircounts, geometry, and quasi Monte-Carlo

The set of the $N$ data points (e.g. galaxies) is $\\{\bm{x}_{i}\\}_{i=1}^{N}$
with all points $\bm{x}_{i}\in W$ inside the observation window
$W\subset\mathbb{R}^{3}$, i.e. inside the unmasked area. The number density is
estimated with $\widehat{\rho}=\tfrac{N}{|W|}$, where $|W|$ is the volume of
$W$. Then we define

$\textsf{DD}(r)=\frac{1}{N^{2}}\sum_{i=1}^{N}\sum_{j=1,j\neq i}^{N}\ \delta
k_{r}^{\delta}(|\bm{x}_{i}-\bm{x}_{j}|),$ (2)

the normalised number of data–data pairs with a distance of
$r=|\bm{x}_{i}-\bm{x}_{j}|$ in the interval $[r,r+\delta]$. We use a
rectangular kernel

$k_{r}^{\delta}(s)=\tfrac{1}{\delta}\mathbb{1}_{[r,r+\delta]}(s),$ (3)

with the indicator function of the set $A$ defined as

$\mathbb{1}_{A}(q)=\begin{cases}1&\text{if }q\in A,\\\
0&\text{else}.\end{cases}$ (4)

Also other kernels with $\int k_{r}^{\delta}(s)\mathrm{d}s=1$ are possible
(e.g. triangular, truncated Gaussian, or Epanechnikov). We consider
${N_{\mathrm{r}}}$ randomly distributed points
$\\{\bm{y}_{j}\\}_{j=1}^{{N_{\mathrm{r}}}}$, all inside the sample geometry
$\bm{y}_{j}\in W$. The normalised number of data–random pairs with a distance
in $[r,r+\delta]$ is denoted by

$\textsf{DR}(r)=\frac{1}{N{N_{\mathrm{r}}}}\sum_{i=1}^{N}\
\sum_{j=1}^{{N_{\mathrm{r}}}}\ \delta k_{r}^{\delta}(|\bm{x}_{i}-\bm{y}_{j}|)$
(5)

Similarly,

$\textsf{RR}(r)=\frac{1}{{N_{\mathrm{r}}}^{2}}\sum_{i=1}^{{N_{\mathrm{r}}}}\sum_{j=1,j\neq
i}^{{N_{\mathrm{r}}}}\ \delta k_{r}^{\delta}(|\bm{y}_{i}-\bm{y}_{j}|),$ (6)

is the normalised number of random–random pairs. The Landy & Szalay
(1993)-estimator is defined by

$\displaystyle\xi_{\text{LS}}(r)$
$\displaystyle=\frac{\textsf{DD}(r)-2\textsf{DR}(r)+\textsf{RR}(r)}{\textsf{RR}(r)}.$
(7)

Also the estimators due to Peebles & Hauser (1974); Hewett (1982); Davis &
Peebles (1983); Hamilton (1993) can be defined in terms of the pair-counts and
our results apply accordingly.

###  2.1 Geometry of pair-counts

The expectation of the pair-counts $\textsf{DR},\textsf{RR}$ can be expressed
in terms of geometric quantities depending on the sample window $W$ and on the
point set (for DR, Kerscher 1999). First consider the set–covariance

$\displaystyle\gamma_{W}(\bm{x})$ $\displaystyle=|W\cap
W_{\bm{x}}|=\int_{\mathbb{R}^{3}}\mathbbold{1}_{W}(\bm{y})\mathbbold{1}_{W}(\bm{y}+\bm{x})\mathrm{d}\bm{y}.$
(8)

where $W_{\bm{x}}$ is the shifted window $W$, i.e. the set of all points from
$W$ shifted by the vector $\bm{x}$. $|W\cap W_{\bm{x}}|$ is the volume of the
set $W\cap W_{\bm{x}}$. The isotropised set–covariance
$\overline{\gamma_{W}}(r)$ can be calculated from $\gamma_{W}(\bm{x})$:

$\overline{\gamma_{W}}(r)=\frac{1}{4\pi}\int_{0}^{\pi}\int_{0}^{2\pi}\,\gamma_{W}(\bm{x}(r,\theta,\phi))\
\sin(\theta)\mathrm{d}\theta\mathrm{d}\phi.$ (9)

Here
$\bm{x}(r,\theta,\phi)=(r\cos(\phi)\sin(\theta),r\sin(\phi)\sin(\theta),r\cos(\theta))$.
For a large number of random points ${N_{\mathrm{r}}}$ one obtains (see
Appendix A):

$\displaystyle\textsf{RR}(r)\rightarrow\ \mathcal{RR}(r)$
$\displaystyle=\frac{4\pi}{|W|^{2}}\int_{r}^{r+\delta}\\!\\!\overline{\gamma_{W}}(s)\
s^{2}\mathrm{d}s$ (10) $\displaystyle\approx\frac{4\pi r^{2}\delta}{|W|^{2}}\
\overline{\gamma_{W}}(r)\quad\text{for }\delta\text{ small}.$

$\mathcal{RR}(r)$ is the expectation value of the pair-count $\textsf{RR}(r)$
illustrating its geometric nature.

Ripley (1988) used a local area weight in an estimator for his $K$-function
(the normalised integrated two-point density) and Rivolo (1986) considered a
similar weight in his estimator for the two-point correlation function. This
weight is inversely proportional to the fraction of the surface area of a
sphere $B_{r}(\bm{y})$ with radius $r$ centred on the point $\bm{y}$ inside
$W$:

$\text{area}(\partial B_{r}(\bm{y})\cap
W)=\int_{0}^{\pi}\int_{0}^{2\pi}\,\mathbbold{1}_{W}(\bm{y}+\bm{x}(r,\theta,\phi))\
\sin(\theta)\mathrm{d}\theta\mathrm{d}\phi.$ (11)

For a large number of random points ${N_{\mathrm{r}}}$ one obtains (see
Appendix A):

$\displaystyle\textsf{DR}(r)\rightarrow\ \mathcal{DR}(r)$
$\displaystyle=\frac{1}{|W|N}\sum_{i=1}^{N}\int_{r}^{r+\delta}\\!\\!\\!\\!\\!\text{area}(\partial
B_{s}(\bm{x}_{i})\cap W)\,\mathrm{d}s$ (12)
$\displaystyle\approx\frac{1}{|W|N}\sum_{i=1}^{N}\text{area}(\partial
B_{r}(\bm{x}_{i})\cap W)\,\delta,\ \text{for }\delta\text{ small}.$

As before, $\mathcal{DR}(r)$ is the expectation value of the pair-count
$\textsf{DR}(r)$ illustrating its geometric nature. But now both
$\textsf{DR}(r)$ and $\mathcal{DR}(r)$ are depending on the points
$\\{\bm{x}_{i}\\}_{i=1}^{N}$ under consideration.

In Appendix B we give expressions for $\overline{\gamma_{W}}(r)$ and
$\text{area}(\partial B_{r}(\bm{y})\cap W)$ if $W$ is a rectangular box or a
sphere. From these expressions we calculate the reference values
$\mathcal{DR}$ and $\mathcal{RR}$. This allows us to compare different
integration schemes for the pair-counts DR and RR and we can investigate the
scaling of the accuracy with the number of points used in these methods.

###  2.2 Quasi Monte-Carlo

In our improved methods for estimating the pair-counts we use quasi Monte-
Carlo integration. In a standard Monte-Carlo integration scheme one uses
${N_{\mathrm{r}}}$ random points
$\\{\bm{y}_{1},\ldots,\bm{y}_{N_{\mathrm{r}}}\\}$ to estimate the integral
$\int_{[0,1]^{d}}f(\bm{x})\mathrm{d}\bm{x}$ by
$\frac{1}{{N_{\mathrm{r}}}}\sum_{i=1}^{N_{\mathrm{r}}}f(\bm{y}_{i})$. The
accuracy can be estimated using the Chebyshev inequality, which tells us that
the probability of an error exceeding a given threshold decreases with
$1/\\!\sqrt{{N_{\mathrm{r}}}}$. In other words, the standard error of a Monte-
Carlo integration scales as $1/\\!\sqrt{{N_{\mathrm{r}}}}$

With quasi Monte-Carlo methods we use
$\frac{1}{{N_{\mathrm{q}}}}\sum_{i=1}^{N_{\mathrm{q}}}f(\bm{q}_{i})$ to
numerically integrate $\int_{[0,1]^{d}}f(\bm{x})\mathrm{d}\bm{x}$ (see e.g.
Niederreiter 1992). This estimate of the integral almost looks identical to
the Monte-Carlo integration above. But for a quasi Monte-Carlo integration the
points in $Q=\\{\bm{q}_{1},\ldots,\bm{q}_{N_{\mathrm{q}}}\\}$ are not random.
It is essential for the application of quasi Monte-Carlo integration, that for
a given point set $Q$ the error bound

$\left|\int_{[0,1]^{d}}f(\bm{x})\mathrm{d}\bm{x}-\frac{1}{{N_{\mathrm{q}}}}\sum_{i=1}^{N_{\mathrm{q}}}f(\bm{q}_{i})\right|\leq
V(f)\,D(Q).$ (13)

factorises into a measure of variation $V(f)$ — only depending on properties
of $f$, and a measure of discrepancy $D(Q)$ — only depending on the properties
of the point set $Q$. If we consider functions of bounded variation the
equation (13) is called the Koksama-Hlawka bound and the measure of
discrepancy is the star-discrepancy (see e.g. L’Ecuyer & Lemieux 2002). To
control the error bound (13), we have to control $D(Q)$ ($V(f)$ is not
depending on the point set $Q$). Low discrepancy sequences, like the Halton
sequence, have been constructed with that in mind. For such sequences Halton
(1960) showed that

$D(Q)\propto\frac{(\log{N_{\mathrm{q}}})^{d}}{{N_{\mathrm{q}}}}.$ (14)

For small dimensions $d$ this compares favourable to a Monte-Carlo integration
where the standard error only scales as $1/\\!\sqrt{{N_{\mathrm{r}}}}$. Halton
sequences can be used to estimate integrals over indicator functions which in
turn define volumes like the set-covariance. Since indicator function have
bounded total variation Eq. (13) applies.

Upper bounds like Eq. (13) are worst case bounds. Owen & Rudolf (2021) derive
an analogue to a strong law of large numbers for randomised low discrepancy
sequences. This further justifies the procedure for scrambling the Halton
sequence developed by Owen (2017), where the scaling from Eq. (14) in Eq. (13)
still gives an upper bound, but on average smaller errors are expected. We
will use these randomised Halton sequences in our calculations (see Appendix
C).

##  3 RR

First we investigate how the accuracy of the standard RR, as given in Eq. (6),
scales with the number of random points used. The expectation value
$\mathcal{RR}$ of RR can be expressed in terms of the isotropised set-
covariance, see Eq. (10). Using the Eqs. (53),(55) from Appendix B.1 for the
isotropised set covariance we calculate $\mathcal{RR}(r)$ as a reference value
for rectangular boxes $W$.

Figure 1: The relative error $\varepsilon_{\textsf{RR}}$ against the number of
random points ${N_{\mathrm{r}}}$ used in the standard procedure for
calculating $\textsf{RR}(r)$ in rectangular windows $W$. The black dotted line
is proportional to $1/\\!\sqrt{{N_{\mathrm{r}}}}$.

In Fig. 1 we show the scaling of the relative error
$\varepsilon_{\textsf{RR}}=|\mathcal{RR}(r)-\textsf{RR}(r)|/\mathcal{RR}(r)$
with the number of random points used. The error $\varepsilon_{\textsf{RR}}$
is the mean value calculated from 100 samples of the random point sets. Only
for the largest ${N_{\mathrm{r}}}$ we use 10 samples. This applies accordingly
to the errors calculated in the following sections. We show the results for a
cuboid $W=[0,1]^{3}$ and for a thin rectangular box
$W=[0,1]\times[0,0.1]^{2}$. With the thin box we investigate how the error is
affected in a sample $W$ where the boundary effects are more dominant. As
expected, for a large number of random points ${N_{\mathrm{r}}}$, the error
scales as $1/\\!\sqrt{{N_{\mathrm{r}}}}$. For small ${N_{\mathrm{r}}}$ and
small radii we see larger errors, especially in $W=[0,1]^{3}$. This increased
error is due to the bias from ”zero-counts” in some of the shells. Simply
using a three dimensional low discrepancy sequence instead of the random
points in Eq. (6) is not feasible as we will see below.

From Eqs. (A) and (38) we get

$\mathcal{RR}(r)=\frac{1}{|W|^{2}}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\mathbb{1}_{W}(\bm{y})\mathbb{1}_{W}(\bm{z})\,\delta
k_{r}^{\delta}(|\bm{y}-\bm{z}|)\,\mathrm{d}\bm{y}\,\mathrm{d}\bm{z}.$ (15)

This allows a more flexible approach. Consider two (random) point sets
$P_{1}=\\{\bm{y}_{i}\\}_{i=1}^{{N_{\mathrm{r}}}}$ and
$P_{2}=\\{\bm{z}_{j}\\}_{j=1}^{{N_{\mathrm{r}}}}$, whith the points
$\bm{y}_{i}\in W$, $\bm{z}_{j}\in W$ all inside the sample geometry. Then we
define

$\displaystyle\textsf{RR}_{2x3d}(r)$
$\displaystyle=\frac{1}{{N_{\mathrm{r}}}^{2}}\sum_{i=1}^{{N_{\mathrm{r}}}}\sum_{j=1}^{{N_{\mathrm{r}}}}\
\delta k_{r}^{\delta}(|\bm{y}_{i}-\bm{z}_{j}|),$ (16)

which is also an estimate of $\mathcal{RR}(r)$. If the points in $P_{1}$ and
$P_{2}$ are drawn from a Poisson process, the estimate from Eq. (16) is
essentially the same as the estimate from Eq. (6), since points in a Poisson
process are independent. See also Dávila-Kurbán et al. (2021), who observe
that one cannot use $P_{1}\equiv P_{2}$ for $\textsf{RR}_{2x3d}$ when they
construct estimator for the two-point correlation function using ‘glass like’
point sets. The representation of $\mathcal{RR}$ in Eq. (15) as a six-
dimensional integral suggest a six-dimensional (quasi) Monte Carlo approach.
Consequently we use six-dimensional random points or a six-dimensional
randomised Halton sequence $\\{\bm{q}_{i}\\}_{i=1}^{{N_{\mathrm{q}}}}$, which
we split as $\bm{q}_{i}=(\bm{y}_{i},\bm{z}_{i})$ into two three-dimensional
sequences. We scale the points in the sequences
$\\{\bm{y}_{i}\\}_{i=1}^{{N_{\mathrm{q}}}}$ and
$\\{\bm{z}_{i}\\}_{i=1}^{{N_{\mathrm{q}}}}$ such that each $\bm{y}_{i}\in W$
and $\bm{z}_{j}\in W$ are uniformly distributed inside the sample geometry.

In Fig. 2 we compare the scaling of the relative error
$\varepsilon_{\textsf{RR},2x3d}=|\mathcal{RR}(r)-\textsf{RR}_{2x3d}(r)|/\mathcal{RR}(r)$
with the number of (quasi) random points. For random points we see the
expected scaling $1/\\!\sqrt{N}$ of the the error. Using the quasi Monte Carlo
approach we see that the error scales proportional to $1/N$, which is even
faster than the theoretical expectation according to Eq. (14). Using a low
discrepancy sequence we gain more than two orders of magnitude in the accuracy
compared to the random point sets.

Figure 2: The relative error $\varepsilon_{\textsf{RR},2x3d}$ calculated for
rectangular windows $W$ with a pure Monte Carlo integration (solid lines) and
with a randomised Halton sequence (dashed lines). The black dotted lines are
proportional to $1/\\!\sqrt{{N_{\mathrm{r}}}}$ (upper) and
$1/{N_{\mathrm{q}}}$ (lower).

##  4 DR

In Sect. 2.1 we show that the expectation of $\textsf{DR}(r)$ is
$\mathcal{DR}(r)$ and that it can be calculated from the area fraction
$\text{area}(\partial B_{r}(\bm{x}_{i})\cap W)$ . For a rectangular box
Baddeley et al. (1993) gave explicit expressions for $\text{area}(\partial
B_{r}(\bm{x}_{i})\cap W)$ (see Appendix B.1). We will use these expressions to
calculate $\mathcal{DR}(r)$ according to Eqs. (42), (43). In the numerical
integration of Eq. (42) we make sure to achieve a relative error of at least
$10^{-10}$ for our reference value $\mathcal{DR}(r)$. Both $\textsf{DR}(r)$
and $\mathcal{DR}(r)$ are depending on the point set under consideration. We
need a realistic data set $\textsf{D}=\\{\bm{x}_{i}\\}_{i=1}^{N}$ to calculate
both the theoretical reference $\mathcal{DR}(r)$ and the different estimates
for $\textsf{DR}(r)$. For this purpose we use a sample of simulated galaxy
clusters111Specifically we use the simulated galaxy clusters from the snapshot
Box2/hr, snap_136, z=0.066340191 downloaded from
http://www.magneticum.org/data.html\\#FULL_CATALOUGES. from the Magneticum
simulation (Hirschmann et al., 2014; Ragagnin et al., 2017). The side length
of the simulation box is $325\,\text{Mpc}/h$ and we use the real space
positions of the 10429 simulated clusters. We rescale the coordinates by the
side-length of the box such that all points are inside $W=[0,1]^{3}$. In the
smaller window $W=[0,1]\times[0,0.1]^{2}$ only 86 clusters are left.

Now we are set to determine the scaling of the relative error
$\varepsilon_{\textsf{DR}}=|\mathcal{DR}(r)-\textsf{DR}(r)|/\mathcal{DR}(r)$
with the number of (quasi) random points used. In Fig. 3 we compare the
standard approach with ordinary random numbers in $\textsf{DR}(r)$ to a pair-
count determined using a low discrepancy sequence instead of the random
points. Using a randomised Halton sequence we only observe a minor gain for
small $N$, but for large $N$ the error is reduced by an order of magnitude. In
the small window $W=[0,1]\times[0,0.1]^{2}$ this is more pronounced. The
scaling follows the expected behaviour from Eq. (14).

Figure 3:  The relative error $\varepsilon_{\textsf{DR}}$ calculated for
rectangular windows $W$ with a standard Monte Carlo integration (solid lines)
and with a quasi Monte Carlo scheme using a randomised Halton sequence (dashed
lines). The black dotted line are proportional to $1/\\!\sqrt{N}$ (upper) and
$(\log N)^{3}/N$ (lower).

We can do better if we consider the geometry of $\textsf{DR}(r)$. From Eq.
(42) and Eq. (43) we know that expectation value
$\mathcal{DR}(r)=\frac{1}{|W|N}\sum_{i=1}^{N}\mathcal{V}_{r}^{\delta}(\bm{x}_{i})$,
where

$\displaystyle\mathcal{V}_{r}^{\delta}(\bm{x}_{i})$
$\displaystyle=\int_{r}^{r+\delta}\\!\\!\\!\\!\\!\text{area}(\partial
B_{s}(\bm{x}_{i})\cap W)\,\mathrm{d}s=|S_{r}^{\delta}(\bm{x}_{i})\cap W|$ (17)

is the volume of the spherical shell with a radial range in $[r,r+\delta]$
around $\bm{x}_{i}$ inside the sample geometry $W$. As already suggested by
Rivolo (1986) this directly leads to a Monte Carlo scheme. With $M$ points
$\\{\bm{y}_{i}\\}_{i=1}^{M}$ (quasi) randomly distributed in the shell
$S_{r}^{\delta}(\bm{x}_{i})=\\{\bm{y}\in\mathbb{R}^{3}\,|\,s<|\bm{y}-\bm{x}_{i}|\leq
s+\delta\\}$ around $\bm{x}_{i}$ we define a (quasi) Mote Carlo estimate of
$\mathcal{V}_{r}^{\delta}(\bm{x}_{i})$

$\displaystyle\textsf{V}_{s}^{\delta}(\bm{x}_{i})$
$\displaystyle=\frac{|S_{r}^{\delta}|}{M}\sum_{j=1}^{M}\mathbb{1}_{W}(\bm{y}_{j})$
(18)

with the volume
$|S_{r}^{\delta}|=\frac{4\pi}{3}\left((r+\delta)^{3}-r^{3}\right)$ of the
shell. Then we get
$\textsf{V}_{s}^{\delta}(\bm{x}_{i})\longrightarrow\mathcal{V}_{s}^{\delta}(\bm{x}_{i})$
for a large number $M$ of (quasi) random points. Consequently we compare

$\displaystyle\textsf{DR}_{\text{shell}}(r)=\frac{1}{|W|N}\sum_{i=1}^{N}\textsf{V}_{r}^{\delta}(\bm{x}_{i})$
(19)

with $\mathcal{DR}(r)$. This is not a pair-count, but we still have a double
sum over $N$ data points and now $M$ points in the shell.

In Fig. 4 we show
$\varepsilon_{\textsf{DR},\text{shell}}=|\mathcal{DR}(r)-\textsf{DR}_{\text{shell}}(r)|/\mathcal{DR}(r)$.
A comparison with DR in Fig. 3 shows that using $\textsf{DR}_{\text{shell}}$
leads to a reduction of the error by almost two orders of magnitude even for
ordinary random points. Again this can be improved by using a low discrepancy
sequences. Then we additionally gain at least another order of magnitude.

Figure 4:  The relative error $\varepsilon_{\textsf{DR},\text{shell}}$
calculated for rectangular windows $W$ with a pure Monte Carlo integration
(solid lines) and with a randomised Halton sequence (dashed lines). The black
dotted line are proportional to $1/\sqrt{M}$ (upper) and $(\log M)^{3}/M$
(lower).

##  5 Estimating $\xi$

Now we join the improved pair-count estimates together and compare results
from the standard Landy & Szalay (1993)-estimator to our new estimator for the
pair correlation function. We use the simulated galaxy clusters from the
Magneticum simulation to illustrate the behaviour of the different estimators
(Hirschmann et al. 2014; Ragagnin et al. 2017, compare also Sect. 4). The
10429 clusters are in a box with side-length $325\text{Mpc}/h$ and we use
their real space positions. The standard Landy & Szalay (1993)-estimator was
defined in Eq. (7):

$\xi_{\text{LS}}(r)=\frac{\textsf{DD}(r)-2\textsf{DR}(r)+\textsf{RR}(r)}{\textsf{RR}(r)}.$

We will use the same ${N_{\mathrm{r}}}$ random points to calculate RR and DR
(see Eqs. (5) and (6)). As an exact reference we have

$\displaystyle\Xi(r)=\frac{\textsf{DD}(r)-2\mathcal{DR}(r)+\mathcal{RR}(r)}{\mathcal{RR}(r)}.$
(20)

with the $\mathcal{RR}$ and $\mathcal{DR}$ given in Eqs. (10), (12), using the
results from Appendix B.1 for a rectangular box. In Sect. 3 and 4 we discuss
alternative possibilities to calculate the pair-counts $DR$ and $RR$. We focus
on the following combination

$\displaystyle\widetilde{\xi}_{\text{LS}}(r)$
$\displaystyle=\frac{\textsf{DD}(r)-2\textsf{DR}_{\text{shell}}(r)+\textsf{RR}_{\text{2x3d}}(r)}{\textsf{RR}_{\text{2x3d}}(r)},$
(21)

resembling the Landy & Szalay (1993)-estimator, but now with improved pair-
count estimates. We use a low discrepancy sequence with $N_{\text{shell}}$
3d-points to calculate $\textsf{DR}_{\text{shell}}$ (see Eq. 19), and another
6d low discrepancy sequence with $N_{\text{2x3d}}$ points to calculate
$\textsf{RR}_{\text{2x3d}}$ (see Eq. 16).

As an illustrative example we compare these estimates in Fig. 5 using an
insufficient number of points. Not unexpected we observe that a standard Landy
& Szalay (1993)-estimator with only $10^{4}$ random points shows deviations
from the exact result. But with $10^{5}$ random points the LS-estimator starts
to follow the exact result. Visually one can see that the estimator using the
low discrepancy sequences achieves a higher accuracy already with
$N_{\text{2x3d}}=10^{4}=N_{\text{shell}}$ points.

Figure 5:  The two-point correlation function calculated for a simulated
galaxy cluster sample. The exact $\Xi$ is compared to standard
$\xi_{\text{LS}}$ with ${N_{\mathrm{r}}}=10^{4}$ and ${N_{\mathrm{r}}}=10^{5}$
random points, and to the new $\widetilde{\xi}_{\text{LS}}$ using
$N_{\text{2x3d}}=10^{4}=N_{\text{shell}}$ points from a randomised Halton
sequence.

Now we make this quantitative and investigate the absolute error
$\Delta_{\xi}(r)=|\Xi(r)-\xi_{\text{LS}}(r)|$. In Fig. 6 we see that the error
of the standard Landy & Szalay (1993)-estimator from Eq. (7) scales as
$1/\sqrt{{N_{\mathrm{r}}}}$. We compare this to the error of the new estimator
Eq. (21) with fixed $N_{\text{shell}}=10^{5}$. In Fig. 6 we see observe that
the error scales as $1/N_{\text{2x3d}}$. Depending on the number of points
used, the absolute error from the new estimator is smaller than the standard
estimator by one up to three orders of magnitude. For
$r\in[100,105]\text{Mpc}/h$ we see that the error of
$\widetilde{\xi}_{\text{LS}}$ shows a turnover towards a constant error for
$N_{\text{2x3d}}>10^{6}$. This is an effect from $\textsf{DR}_{\text{shell}}$.
We would have to increase $N_{\text{shell}}$ to obtain a smaller error. We
checked that an $N_{\text{shell}}=10^{5}$ is sufficient to obtain the required
accuracy for the rest of the results as shown in Fig. 6. For a reduced
accuracy goal and a small $N_{\text{2x3d}}$ also a smaller $N_{\text{shell}}$
can be used, but we did not optimise this.

Figure 6:  The absolute error $\Delta_{\xi}$ of the standard
$\xi_{\texttt{LS}}$ is shown against ${N_{\mathrm{r}}}$ (solid lines). We
compare with the absolute error $\Delta_{\xi}$ from
$\widetilde{\xi}_{\texttt{LS}}$ plotted against $N_{\text{2x3d}}$ (dashed
lines). The black dotted line are proportional to $1/\sqrt{{N_{\mathrm{r}}}}$
(upper) and $(\log N_{\text{2x3d}})^{3}/N_{\text{2x3d}}$ (lower).

Ultimately we are not interested in the scaling of the error with the number
of (quasi) random points. We want to obtain reliable results _fast_. The speed
of the calculation depends on the hardware platform, on the algorithms, on
details of the implementation and further effects. Nevertheless we prepare an
example to investigate what one can expect as a speedup. Our goal is to obtain
an absolute error of $\Delta_{\xi}\leq 10^{-4}$ for a radius in
$[30,35]\,\text{Mpc}/h$. From Fig. 6 we reed off that we need at least
${N_{\mathrm{r}}}=10^{7}$ random points for the standard $\xi_{\text{LS}}$,
whereas $N_{\text{2x3d}}=10^{5}$ points from a low discrepancy sequences are
sufficient with $\widetilde{\xi}_{\text{LS}}$. As above we checked that
$N_{\text{shell}}=10^{5}$ is sufficient to obtain the required accuracy. We
ran our mildly optimised code (see Appendix C) on a current small
workstation222This workstation was equipped with an Intel Xeon W1350 processor
with six cores and a clock frequency of at least 3,3 Ghz. Memory was not an
issue in these calculations. With the above parameters the standard
$\xi_{\text{LS}}$ calculation took about 19 hours. The improved
$\widetilde{\xi}_{\text{LS}}$ finished in less than 0.01 hours (less then 36
sec). We observe a speed improvement by a factor larger than $10^{4}$. With a
similar approach for the radius in $[100,105]\,\text{Mpc}/h$ we observe a
speedup by a factor of $10^{3}$.

##  6 Incomplete sampling

The random points in RR and DR are not only used to correct for finite-size
effects but also to correct for incomplete sampling. We will follow the ideas
presented in Baddeley et al. (2000) and Shaw et al. (2021) and adapt them for
the pair-counts. In this way we still can use the new improved estimators and
also correct for incomplete sampling.

Typically the galaxy distribution is observed incompletely. For example
unobserved regions around bright stars are masked. This leads to holes in the
observational window $W$ and can be dealt with the methods already described.
But a partial sampling of the galaxy distribution in crowded fields leaves us
with an inhomogeneous selection of the galaxy sample. Also in magnitude
limited samples we have a systematic selection of the galaxies, depending on
the distance from us. We model this by $p(\bm{x})\in(0,1]$, the probability of
including/observing a galaxy at position $\bm{x}$. We assume that this
$p(\bm{x})$ is statistically independent from any other point in galaxy
distribution. For example, in a crowded field only a fraction of the galaxies
is targeted for spectroscopic follow up. For a galaxy in this field we have a
$p(\bm{x})$ equal to the fraction of targeted galaxies in the field. In a
magnitude limited sample we observe all galaxies down to a limiting brightness
$l_{\text{lim}}$. The luminostiy of a galaxy $L(d(\bm{x}),l)$ at position
$\bm{x}$ can be determined from the brightness $l>l_{\text{lim}}$ and its
luminosity distance $d(\bm{x})$ from our galaxy (for simplicity we neglect
absorbtion). Hence at a distance of $d$ we only include galaxies with
luminosity $L>L(d,l_{\text{lim}})$ in our catalogue. Now consider $F[L]$ the
distribution function of the absolute luminosities of all the galaxies (i.e.
the normalised cumulative luminosity function). The fraction of galaxies
included at a distance $d(\bm{x})$ is then
$p(\bm{x})=1-F[L(d(\bm{x}),l_{\text{lim}})]$. Hence we need a good model for
the luminosity function. The exclusion of a galaxy due to fibre collision
cannot be modelled with the independent thinning (see below), and the
following approach may only serve as an approximation. Further selection and
sampling effects are discussed in Ross et al. (2012).

To investigate the effects of incomplete sampling we assume an unobserved
homogeneous and isotropic galaxy distribution $G$ with number density
$\varrho_{G}$ and two-point correlation function $\xi_{G}(r)$. Our goal is to
estimate this $\xi_{G}(r)$. The inhomogeneous sampling, described by
$p(\bm{x})$, leads to an observed inhomogeneous galaxy catalogue and is
modelled in the following way:

$\textsf{D}^{\prime}=\\{\bm{x}_{i}\in G|\bm{x}_{i}\in W\text{ and }u_{i}\leq
p(\bm{x}_{i})\\},$ (22)

where $u_{i}$ are independent random variables, specifically independent from
the points, uniformly distributed on $[0,1]$. This closely follows the
construction of an inhomogeneous Markov point process by independent thinning
as discussed by Baddeley et al. (2000). The observed inhomogeneous point set
$\\{\bm{x}_{i}^{\prime}\\}_{i=1}^{N}$ is considered a realisation of
$\textsf{D}^{\prime}$. With the number density $\varrho_{G}$ of the
homogeneous galaxy distribution, the inhomogeneous number density of this
point process is $\varrho^{\prime}(\bm{x})=p(\bm{x})\,\varrho_{G}$. The two-
point density, the probability of observing a point at $\bm{x}$ and $\bm{y}$,
is

$\displaystyle\varrho_{2}^{\prime}(\bm{x},\bm{y})=\varrho^{\prime}(\bm{x})\varrho^{\prime}(\bm{y})\left(1+\xi^{\prime}(|\bm{x}-\bm{y}|)\right).$
(23)

This is well defined for a certain class of point process as defined in
Baddeley et al. (2000). Since the thinning is assumed to be independent from
the points we have

$\xi_{G}(|\bm{x}-\bm{y}|)=\xi^{\prime}(|\bm{x}-\bm{y}|)=\frac{\varrho_{2}^{\prime}(\bm{x},\bm{y})}{\varrho_{G}^{2}\
p(\bm{x})p(\bm{y})}-1,$

for the two-point correlation function.

In full analogy to Sec. 2 we define the ‘inhomogeneous’ pair counts
$\textsf{DD}^{\prime}$. The points $\\{\bm{x}_{i}^{\prime}\\}_{i=1}^{N}$ are
the incompletely sampled galaxies with positions $\bm{x}_{i}^{\prime}$ inside
the observation window $W$. Then

$\textsf{DD}^{\prime}(r)=\frac{1}{N^{2}}\sum_{i=1}^{N}\sum_{j=1,j\neq i}^{N}\
\delta k_{r}^{\delta}(|\bm{x}_{i}^{\prime}-\bm{x}_{j}^{\prime}|).$ (24)

To account for the effect of incomplete sampling one often applies the
thinning as described in Eq. (22) to the random points too. In this way our
random point set $\\{\bm{y}_{i}^{\prime}\\}_{i=1}^{N_{\mathrm{r}}}$ is a
realisation of an inhomogeneous Poisson process with number density
$\tfrac{{N_{\mathrm{r}}}}{N}\varrho^{\prime}(\bm{x})$. The pair-counts
involving these inhomogeneous random points are

$\textsf{RR}^{\prime}(r)=\frac{1}{{N_{\mathrm{r}}}^{2}}\sum_{i=1}^{{N_{\mathrm{r}}}}\sum_{j=1,j\neq
i}^{{N_{\mathrm{r}}}}\ \delta
k_{r}^{\delta}(|\bm{y}_{i}^{\prime}-\bm{y}_{j}^{\prime}|),$ (25)

and

$\textsf{DR}^{\prime}(r)=\frac{1}{N{N_{\mathrm{r}}}}\sum_{i=1}^{N}\
\sum_{j=1}^{{N_{\mathrm{r}}}}\ \delta
k_{r}^{\delta}(|\bm{x}_{i}^{\prime}-\bm{y}_{j}^{\prime}|).$ (26)

In Appendix A.2 we show, that the Peebles & Hauser (1974)-estimator using
$\textsf{DD}^{\prime}$ and $\textsf{RR}^{\prime}$ is ‘ratio unbiased’. For the
Landy & Szalay (1993)-estimator we combine these pair-counts
$\textsf{DD}^{\prime}$, $\textsf{DR}^{\prime}$ and $\textsf{RR}^{\prime}$ to

$\displaystyle\xi_{\text{LS}}^{\prime}(r)$
$\displaystyle=\frac{\textsf{DD}^{\prime}(r)-2\textsf{DR}^{\prime}(r)+\textsf{RR}^{\prime}(r)}{\textsf{RR}^{\prime}(r)},$
(27)

in full analogy to Eq. (7). It remains to show how we can calculate these
inhomogeneous pair-counts using our new approach based on quasi Monte Carlo
methods. For this we make a detour and consider the expectation value of
$\textsf{RR}^{\prime}$ and then $\textsf{DR}^{\prime}$.

For a large number of random points ${N_{\mathrm{r}}}$ from an inhomogeneous
Poisson process we obtain (see Appendix A.2):

$\displaystyle\textsf{RR}^{\prime}(r)$ $\displaystyle\rightarrow\
\mathcal{RR}^{\prime}(r)=\frac{4\pi}{N^{2}}\int_{r}^{r+\delta}\overline{\Gamma_{W}}(s)\
s^{2}\mathrm{d}s$ (28)

where $\overline{\Gamma_{W}}(s)$ is given in Eq. (49).
$\mathcal{RR}^{\prime}(r)$ is the expectation value of the pair-count
$\textsf{RR}^{\prime}(r)$. There is no longer a direct correspondence with
geometric objects like the set-covariance, however we can write
$\mathcal{RR}^{\prime}$ as the following integral (compare to Eq. (15)):

$\mathcal{RR}^{\prime}(r)=\frac{1}{N^{2}}\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\mathbb{1}_{W}(\bm{y}_{1})\mathbb{1}_{W}(\bm{y}_{2})\,\varrho^{\prime}(\bm{y}_{1})\varrho^{\prime}(\bm{y}_{2})\\\
\delta
k_{r}^{\delta}(|\bm{y}_{1}-\bm{y}_{2}|)\,\mathrm{d}\bm{y}_{1}\,\mathrm{d}\bm{y}_{2}.$
(29)

This suggests a (quasi) Monte Carlo approach similar to Eq. (16). We use
homogeneously sampled point sets $P_{1},P_{2}$ as described in Sec. 3. Then

$\displaystyle\textsf{RR}_{2x3d}^{\prime}(r)$
$\displaystyle=\frac{1}{N^{2}}\sum_{i=1}^{{N_{\mathrm{r}}}}\sum_{j=1}^{{N_{\mathrm{r}}}}\
\varrho^{\prime}(\bm{y}_{i})\varrho^{\prime}(\bm{z}_{j})\ \delta
k_{r}^{\delta}(|\bm{y}_{i}-\bm{z}_{j}|),$ (30)

is an estimate of $\mathcal{RR}^{\prime}(r)$ (remember, the number density of
the inhomogeneous Poisson process is
$\tfrac{{N_{\mathrm{r}}}}{N}\varrho^{\prime}(\bm{x})$). As described in Sec. 3
the point sets $P_{1},P_{2}$ can be generated from a six-dimensional random
sequence or six-dimensional low discrepancy sequence.

In Eq. (26) we use a thinned random point set to calculate
$\textsf{DR}^{\prime}$. Similarly, we could use a thinned low discrepancy
sequence in Eq. (26), but in Sec. 4 we saw that this first approach was only
mildly successful. Hence, we follow the second approach by considering Eq.
(A.2) and find

$\displaystyle\textsf{DR}^{\prime}(r)$ $\displaystyle\rightarrow\
\mathcal{DR}^{\prime}(r)=\frac{1}{N^{2}}\sum_{i=1}^{N}\mathcal{U}_{r}^{\delta}(\bm{x}_{i}),$
(31)

where $\mathcal{U}_{r}^{\delta}(\bm{x}_{i})$ is given in Eq. (52). Similar to
Eq. (18) we use $M$ points $\\{\bm{y}_{i}\\}_{i=1}^{M}$ (quasi) randomly
distributed in the shell
$S_{r}^{\delta}(\bm{x}_{i})=\\{\bm{y}\in\mathbb{R}^{3}\,|\,s<|\bm{y}-\bm{x}_{i}|\leq
s+\delta\\}$ around $\bm{x}_{i}$ to estimate
$\mathcal{U}_{r}^{\delta}(\bm{x}_{i})$ by

$\displaystyle\textsf{U}_{s}^{\delta}(\bm{x}_{i})$
$\displaystyle=\frac{|S_{r}^{\delta}|}{M}\sum_{j=1}^{M}\mathbb{1}_{W}(\bm{y}_{j})\
\varrho^{\prime}(\bm{y}_{j}).$ (32)

This results in a new estimate for $\mathcal{DR}^{\prime}(r)$:

$\displaystyle\textsf{DR}_{\text{shell}}^{\prime}(r)=\frac{1}{N^{2}}\sum_{i=1}^{N}\textsf{U}_{r}^{\delta}(\bm{x}_{i}).$
(33)

Similar to Eq. (21) we can use $\textsf{DD}^{\prime}$,
$\textsf{RR}_{2x3d}^{\prime}$ and $\textsf{DR}_{\text{shell}}^{\prime}$ to
construct a Landy & Szalay (1993)-type estimator
$\widetilde{\xi}_{\text{LS}}^{\prime}$ for inhomogeneously sampled gaxay
catalogues.

Shaw et al. (2021) discuss how to estimate the $\xi^{\prime}(r)$ _and_ a non-
parametric model for $\varrho^{\prime}(\bm{x})$ directly from the data. In a
typical application to galaxy catalogues we often have a good model for
$\varrho^{\prime}(\bm{x})$ (equivalently for $p(\bm{x})$). Parameters of this
model are fixed by the sampling strategy, and some are determined from the
galaxy distribution (e.g. from an estimate of the luminosity function).

##  7 Summary and outlook

First, we focussed on the scaling of the error in estimates of the pair-counts
RR and DR with the number of (quasi) random points. For the standard approach,
with ordinary random numbers, we confirm the expected slow shrinking of the
error proportional to $1/\sqrt{N}$. A reformulation of the pair-counts made a
quasi Monte Carlo integration possible. There we found that the error shrinks
almost proportional to $1/{N_{\mathrm{q}}}$, where ${N_{\mathrm{q}}}$ is the
number of points from a low discrepancy sequence. This scaling of the error
not only holds for bulky samples but is even more pronounced in a thin sample
with prevalent boundary effects. Hence we are confident that our improved
methods are also applicable in more complicated sample geometries. We combine
these improved pair-counts into a new Landy & Szalay (1993)-type estimator and
compare with the standard one. The new estimator inherits the favorable
scaling $1/{N_{\mathrm{q}}}$. We can turn this observation around. For a fixed
maximum error we read off how many points are necessary in the new and in the
standard estimator. In an example we compare the running times of the
estimators and obtaine a speedup by a factor of $10^{4}$ for our improved
estimator.

The random points sets are not only used for boundary corrections but also to
correct for incomplete sampling. Typically the selection and sampling effects,
as present in the galaxy distribution, are modelled onto the random point set
used in the pair-counts. We discuss how to adapt this for the new improved
estimator. Essentially we are using the probability of observing the galaxies
as a weight in the calculation of the pair-counts. This is only a first step.
Using wheighting schemes one could envision the construction of a minimum
variance estimate for inhomogeneously sampled galaxies in the line of Feldman
et al. (1994) and Colombi et al. (1998), but now for the improved estimators
using low discrepancy sequences.

Our usage of the randomised Halton sequence is similar to how Dávila-Kurbán et
al. (2021) use glass-like point sets. Generating glass-like point sets can be
computationally challenging. In contrast a randomised Halton sequence is
generated easily (see Appendix C). Halton sequences are probably the simplest
choice, but also other low discrepancy sequences could be used (L’Ecuyer &
Lemieux, 2002). Interesting alternatives could be so called blue noise random
sets, which are used in computer graphics for efficiently sampling from
surfaces (see e.g. Heck et al. 2013). The computational improvements, as
referenced in the introduction, are complimentary to our approach and can be
used similarly for the low discrepancy sequences. This should lead to a
further speed-up.

###### Acknowledgements.

Many thanks to Adrian Baddeley for sharing code and the comments on the
expressions for the area fraction. I would like to thank Klaus Dolag and
Antonio Ragagnin for providing public access to the simulated galaxy clusters
from the Magneticum simulation.

## References

  * Abbott et al. (2022) Abbott, T. M. C., Aguena, M., Alarcon, A., et al. 2022, Phys. Rev. D, 105, 023520
  * Aghamousa et al. (2016) Aghamousa, A., Aguilar, J., Ahlen, S., et al. 2016, arXiv e-prints, arXiv:1611.00036
  * Alarcon et al. (2021) Alarcon, A., Gaztanaga, E., Eriksen, M., et al. 2021, Mon. Not. Roy. Astron. Soc., 501, 6103
  * Alonso (2012) Alonso, D. 2012, arXiv e-prints, arXiv:1210.1833
  * Amendola et al. (2018) Amendola, L., Appleby, S., Avgoustidis, A., et al. 2018, Living Reviews in Relativity, 21, 2
  * Anderson et al. (1999) Anderson, E., Bai, Z., Bschof, C., et al. 1999, LAPACK Users Guide, 3rd edn. (Philadelphia: SIAM)
  * Baddeley & Turner (2005) Baddeley, A. & Turner, R. 2005, Journal of Statistical Software, 12, 1
  * Baddeley et al. (1993) Baddeley, A. J., Moyeed, R. A., Howard, C. V., & Boyde, A. 1993, Appl. Statist., 42, 641
  * Baddeley et al. (2000) Baddeley, A. J., Møller, J., & Waagepetersen, R. 2000, Statistica Neerlandica, 54, 329
  * Bautista et al. (2021) Bautista, J. E., Paviot, R., Vargas Magaña, M., et al. 2021, Mon. Not. Roy. Astron. Soc., 500, 736
  * Breton & de la Torre (2021) Breton, M.-A. & de la Torre, S. 2021, A&A, 646, A40
  * Colombi et al. (1998) Colombi, S., Szapudi, I., & Szalay, A. S. 1998, Mon. Not. Roy. Astron. Soc., 296, 253
  * Dagum & Menon (1998) Dagum, L. & Menon, R. 1998, Computational Science & Engineering, IEEE, 5, 46
  * Daley & Vere-Jones (2003) Daley, D. J. & Vere-Jones, D. 2003, An Introduction to the Theory of Point Processes (Berlin: Springer Verlag)
  * Dávila-Kurbán et al. (2021) Dávila-Kurbán, F., Sánchez, A. G., Lares, M., & Ruiz, A. N. 2021, Mon. Not. Roy. Astron. Soc., 506, 4667
  * Davis & Peebles (1983) Davis, M. & Peebles, P. J. E. 1983, ApJ, 267, 465
  * Dawson et al. (2013) Dawson, K. S., Schlegel, D. J., Ahn, C. P., et al. 2013, AJ, 145, 10
  * Demina et al. (2018) Demina, R., Cheong, S., BenZvi, S., & Hindrichs, O. 2018, Mon. Not. Roy. Astron. Soc., 480, 49
  * Donoso (2019) Donoso, E. 2019, Mon. Not. Roy. Astron. Soc., 487, 2824
  * Eisenstein et al. (2005) Eisenstein, D. J., Zehavi, I., Hogg, D. W., et al. 2005, ApJ, 633, 560
  * Feldman et al. (1994) Feldman, H. A., Kaiser, N., & Peacock, J. A. 1994, ApJ, 426, 23
  * Fiksel (1988) Fiksel, T. 1988, Statistics, 19, 67
  * Halton (1960) Halton, J. 1960, Numer. Math., 2, 84
  * Hamilton (1993) Hamilton, A. 1993, ApJ, 417, 19
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357
  * He (2021) He, C.-C. 2021, ApJ, 921, 59
  * Heck et al. (2013) Heck, D., Schlömer, T., & Deussen, O. 2013, ACM Transactions on Graphics, 32 (3), 1
  * Hewett (1982) Hewett, P. C. 1982, Mon. Not. Roy. Astron. Soc., 201, 867
  * Hirschmann et al. (2014) Hirschmann, M., Dolag, K., Saro, A., et al. 2014, Mon. Not. Roy. Astron. Soc., 442, 2304
  * Ivezić et al. (2019) Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111
  * Jakob et al. (2017) Jakob, W., Rhinelander, J., & Moldovan, D. 2017, pybind11 – Seamless operability between C++11 and Python, https://github.com/pybind/pybind11
  * Keihänen et al. (2019) Keihänen, E., Kurki-Suonio, H., Lindholm, V., et al. 2019, A&A, 631, A73
  * Kerscher (1999) Kerscher, M. 1999, Astron. Astrophys., 343, 333
  * Kerscher et al. (2000) Kerscher, M., Szapudi, I., & Szalay, A. 2000, Ap. J., 535, L13
  * Landy & Szalay (1993) Landy, S. D. & Szalay, A. S. 1993, ApJ, 412, 64
  * L’Ecuyer & Lemieux (2002) L’Ecuyer, P. & Lemieux, C. 2002, Recent Advances in Randomized Quasi-Monte Carlo Methods (New York, NY: Springer US), 419–474
  * Moore et al. (2001) Moore, A. W., Connolly, A. J., Genovese, C., et al. 2001, in Mining the Sky, ed. A. J. Banday, S. Zaroubi, & M. Bartelmann, 71
  * Neyman & Scott (1958) Neyman, J. & Scott, E. L. 1958, J. R. Stat. Soc., 20, 1
  * Niederreiter (1992) Niederreiter, H. 1992, Random Number Generation and Quasi-Monte Carlo Methods (Philadelphia, Pennsylvania: SIAM)
  * Ohser (1983) Ohser, J. 1983, Math. Operationsforsch. u. Statist., Ser. Statist., 14, 63
  * Owen (2017) Owen, A. B. 2017, arXiv e-prints, arXiv:1706.02808
  * Owen & Rudolf (2021) Owen, A. B. & Rudolf, D. 2021, SIAM Rev., 63(2), 360
  * Peebles (1980) Peebles, P. J. E. 1980, The Large Scale Structure of the Universe (Princeton, New Jersey: Princeton University Press)
  * Peebles & Hauser (1974) Peebles, P. J. E. & Hauser, M. G. 1974, Ap. J. Suppl., 28, 19
  * Ragagnin et al. (2017) Ragagnin, A., Dolag, K., Biffi, V., et al. 2017, Astronomy and Computing, 20, 52
  * Ripley (1976) Ripley, B. D. 1976, J. Appl. Prob., 13, 255
  * Ripley (1988) Ripley, B. D. 1988, Statistical Inference For Spatial Processes (Cambridge: Cambridge University Press)
  * Rivolo (1986) Rivolo, A. R. 1986, ApJ, 301, 70
  * Ross et al. (2020) Ross, A. J., Bautista, J., Tojeiro, R., et al. 2020, Mon. Not. Roy. Astron. Soc., 498, 2354
  * Ross et al. (2012) Ross, A. J., Percival, W. J., Sánchez, A. G., et al. 2012, Mon. Not. Roy. Astron. Soc., 424, 564
  * Shaw et al. (2021) Shaw, T., Møller, J., & Waagepetersen, R. 2021, Australian & New Zealand Journal of Statistics, 63, 93
  * Sinha & Garrison (2020) Sinha, M. & Garrison, L. H. 2020, MNRAS, 491, 3022
  * Stoyan et al. (1995) Stoyan, D., Kendall, W. S., & Mecke, J. 1995, Stochastic Geometry and its Applications, 2nd edn. (Chichester: John Wiley & Sons)
  * Stoyan & Stoyan (1994) Stoyan, D. & Stoyan, H. 1994, Fractals, Random Shapes and Point Fields: Methods of Geometrical Statistics (Chichester: John Wiley & Sons)
  * Stoyan & Stoyan (2000) Stoyan, D. & Stoyan, H. 2000, Scandinavian Journal of Statistics, 27, 641
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261

##  Appendix A Expectation of pair-counts

Estimators for the two-point density and the correlation function using
geometrical weights have been developed in spatial statistics (e.g. Stoyan &
Stoyan 1994). These geometrical weights can be derived with an application of
the Campbell-Mecke formula to the pair-counts. You can find these ideas in
several places, Ripley (1976); Ohser (1983); Fiksel (1988); Stoyan et al.
(1995) and quite recently for inhomogeneous point sets in Shaw et al. (2021).
We give analogous derivations to show the connection between the pair-counts
DR and RR and geometric quantities like the set-covariance and the area
fraction (Kerscher, 1999; Stoyan & Stoyan, 2000).

The Campbell-Mecke formula connects the expectation value over realisations of
a point process $\Phi$ to integrals over $n$-point densities (e.g. Stoyan et
al. 1995; Daley & Vere-Jones 2003). For suitable functions $f(\bm{x})$ and
$g(\bm{x},\bm{y})$ we have

$\displaystyle\mathbb{E}\left[\sum_{\bm{x}\in\Phi}f(\bm{x})\,\right]$
$\displaystyle=\int_{\mathbb{R}^{3}}f(\bm{x})\,\varrho\,\mathrm{d}\bm{x},$
(34)
$\displaystyle\mathbb{E}\left[\sum_{\begin{subarray}{c}\bm{x},\bm{y}\in\Phi\\\
\bm{x}\neq\bm{y}\end{subarray}}g(\bm{x},\bm{y})\,\right]$
$\displaystyle=\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}g(\bm{x},\bm{y})\,\varrho_{2}(\bm{x},\bm{y})\,\mathrm{d}\bm{x}\mathrm{d}\bm{y}.$
(35)

The Campbell-Mecke formula allows us to interchange the expectation of sums
with an integration over the number density $\varrho$ or the two-point density
$\varrho_{2}(\bm{x},\bm{y})$. For a simple point process $\Phi$ we define

$\displaystyle S$
$\displaystyle=\sum_{\begin{subarray}{c}\bm{x},\bm{y}\in\Phi\\\
\bm{x}\neq\bm{y}\end{subarray}}\mathbbold{1}_{W}(\bm{x})\mathbbold{1}_{W}(\bm{y})\delta
k_{r}^{\delta}(|\bm{x}-\bm{y}|).$ (36)

Considering the galaxy distribution as a realisation of a point process $\Phi$
we get $\textsf{DD}=\frac{S}{N^{2}}$. We apply the Campbell-Mecke formula to
$\mathbb{E}[S]$, where we assume homogeneity and isotropy
$\varrho_{2}(\bm{x},\bm{y})=\varrho_{2}(|\bm{x}-\bm{y}|)$.

$\displaystyle\mathbb{E}[S]$
$\displaystyle=\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}\mathbbold{1}_{W}(\bm{x})\mathbbold{1}_{W}(\bm{y})\delta
k_{r}^{\delta}(|\bm{x}-\bm{y}|)\varrho_{2}(|\bm{x}-\bm{y}|)\,\mathrm{d}\bm{x}\mathrm{d}\bm{y}$
$\displaystyle=\int_{\mathbb{R}^{3}}\underbrace{\int_{\mathbb{R}^{3}}\mathbbold{1}_{W}(\bm{z}+\bm{y})\mathbbold{1}_{W}(\bm{y})\mathrm{d}\bm{y}}_{=\gamma_{W}(\bm{z})}\
\delta k_{r}^{\delta}(|\bm{z}|)\varrho_{2}(|\bm{z}|)\,\mathrm{d}\bm{z}$
$\displaystyle=\int_{0}^{\infty}\underbrace{\int_{0}^{\pi}\\!\\!\\!\int_{0}^{2\pi}\\!\\!\\!\gamma_{W}(\bm{z}(s,\theta,\phi))\sin(\theta)\mathrm{d}\theta\mathrm{d}\phi}_{=4\pi\overline{\gamma_{W}}(s)}\
\delta k_{r}^{\delta}(s)\varrho_{2}(s)\,s^{2}\mathrm{d}s$
$\displaystyle=4\pi\int_{r}^{r+\delta}\\!\\!\overline{\gamma_{W}}(s)\
\varrho^{2}(1+\xi(s))\ s^{2}\mathrm{d}s,$ (37)

with $\varrho_{2}(s)=\varrho^{2}(1+\xi(s))$ and the number density $\varrho$
of the point process.

The random point set used in RR is a realisation of a Poisson process with
number density $\varrho_{\text{r}}=\frac{{N_{\mathrm{r}}}}{|W|}$ and by
definition a vanishing two-point correlation function $\xi(r)=0$. Applying the
Campbell-Mecke formula to
$\textsf{RR}=\frac{S_{\text{r}}}{{N_{\mathrm{r}}}^{2}}$ we obtain the
connection between the pair-count RR and the isotropised set-covariance
(Kerscher, 1999; Stoyan & Stoyan, 2000):

$\displaystyle\mathbb{E}[\textsf{RR}(r)]=\mathcal{RR}(r)$
$\displaystyle=\frac{4\pi}{{N_{\mathrm{r}}}^{2}}\int_{r}^{r+\delta}\\!\\!\overline{\gamma_{W}}(s)\
\varrho_{\text{r}}^{2}\ s^{2}\mathrm{d}s$ (38) $\displaystyle\approx\
\frac{4\pi r^{2}\delta}{|W|^{2}}\ \overline{\gamma_{W}}(r)\quad\text{for
}\delta\text{ small}.$

Using the isotropised set covariance (see Appendix B.1) we can calculate
$\mathcal{RR}(r)$ as a geometrical reference value for rectangular boxes.

As mentioned before, we consider the galaxy distribution as a realisation of a
point process with number density $\varrho=\frac{N}{|W|}$ and we seek to
estimate its two-point correlation function $\xi(r)$. From

$\displaystyle\frac{\mathbb{E}[\textsf{DD}(r)]}{\mathbb{E}[\textsf{RR}(r)]}$
$\displaystyle=\frac{{N_{\mathrm{r}}}^{2}\
4\pi\int_{r}^{r+\delta}\overline{\gamma_{W}}(s)\
\varrho^{2}(1+\xi(s))s^{2}\mathrm{d}s}{N^{2}\
4\pi\int_{r}^{r+\delta}\overline{\gamma_{W}}(s)\
\varrho_{\text{r}}^{2}s^{2}\mathrm{d}s}$
$\displaystyle=1+\frac{\int_{r}^{r+\delta}\overline{\gamma_{W}}(s)\ \xi(s)\
s^{2}\mathrm{d}s}{\int_{r}^{r+\delta}\overline{\gamma_{W}}(s)\
s^{2}\mathrm{d}s}$ (39) $\displaystyle\approx\ 1+\xi(r)\quad\text{for
}\delta\text{ small},$

we see that the estimator
$\xi_{\text{PH}}(r)=\frac{\textsf{DD}(r)}{\textsf{RR}(r)}-1$ of Peebles &
Hauser (1974) is a ‘ratio unbiased’ estimator for the two-point correlation
function $\xi(r)$.

To derive an analogous relation between DR and the average surface area we
consider

$\displaystyle T$
$\displaystyle=\sum_{i=1}^{N}\sum_{\bm{y}\in\Phi}\mathbbold{1}_{W}(\bm{y})\,\delta
k_{r}^{\delta}(|\bm{x}_{i}-\bm{y}|).$ (40)

whith a point process $\Phi$ and $\\{\bm{x}_{i}\\}_{i=1}^{N}$ a given set of
points inside the sample geometry $\bm{x}_{i}\in W$. With $\Phi$ a Poisson
process and using the Campbell-Mecke formula for
$\textsf{DR}(r)=\frac{T}{N{N_{\mathrm{r}}}}$ we obtain

$\displaystyle\mathbb{E}[T]$
$\displaystyle=\sum_{i=1}^{N}\int_{\mathbb{R}^{3}}\mathbbold{1}_{W}(\bm{y})\,\delta
k_{r}^{\delta}(|\bm{x}_{i}-\bm{y}|)\,\varrho_{\text{r}}\,\mathrm{d}\bm{y}=\frac{{N_{\mathrm{r}}}}{|W|}\sum_{i=1}^{N}\mathcal{V}_{r}^{\delta}(\bm{x}_{i})$
(41)

where $\varrho_{\text{r}}=\frac{{N_{\mathrm{r}}}}{|W|}$ is the number density
of the Poisson process and

$\displaystyle\mathcal{V}_{r}^{\delta}(\bm{x}_{i})$
$\displaystyle=|S_{r}^{\delta}(\bm{x}_{i})\cap
W|=\int_{r}^{r+\delta}\\!\\!\\!\\!\\!\text{area}(\partial
B_{s}(\bm{x}_{i})\cap W)\,\mathrm{d}s$ (42)

is the volume of the spherical shell
$S_{r}^{\delta}(\bm{x}_{i})=\\{\bm{y}\in\mathbb{R}^{3}\,|\,s<|\bm{y}-\bm{x}_{i}|\leq
s+\delta\\}$ with a radial range in $[r,r+\delta]$ around $\bm{x}_{i}$ inside
the sample geometry $W$. For simple sample geometries the
$\text{area}(\partial B_{s}(\bm{x}_{i})\cap W)$ can be calculated explicitly
(see Appendix B.1). The integral in Eq. (42) can be evaluated using standard
numerical methods to obtain the expectation value of value $\textsf{DR}(r)$:

$\displaystyle\mathcal{DR}(r)$
$\displaystyle=\frac{1}{|W|N}\sum_{i=1}^{N}\mathcal{V}_{r}^{\delta}(\bm{x}_{i})$
(43) $\displaystyle\approx\
\frac{1}{|W|}\,\frac{1}{N}\sum_{i=1}^{N}\text{area}(\partial
B_{r}(\bm{x}_{i})\cap W)\,\delta\quad\text{for }\delta\text{ small}.$

Clearly, the expectation value of the estimators from Davis & Peebles (1983);
Hewett (1982), and Hamilton (1993) can be expressed in terms of the
istotropised set-covariance and the average area fraction (Kerscher, 1999).

###  A.1 Periodic boundaries

Lets assume that our window is a rectangular box
$W=[0,L_{1}]\times[0,L_{2}]\times[0,L_{3}]$ with periodic boundaries, i.e. $W$
has the topology of a three-torus. The majority of cosmological simulations
enforce these boundary conditions. In such a situation no boundary corrections
are needed in the calculation of the two-point correlation function, but we
have to respect the periodicity in each coordinate direction (Stoyan et al.,
1995). The distance between two points $\bm{x},\bm{y}\in W$ is
$d(\bm{x},\bm{y})=\sqrt{d_{1}^{2}+d_{2}^{2}+d_{3}^{2}}$, with
$d_{i}=\min\\{|x_{i}-y_{i}|,L_{i}-|x_{i}-y_{i}|\\}$. This is applicable if $r$
is smaller than any $L_{i}/2$. For the shifted periodic box we have
$W_{\bm{x}}=\\{\bm{y}\,|\,\bm{y}-\bm{x}\in W\\}=W$ and consequently the
(isotropised) set-covariance is constant
$\gamma_{W}(\bm{x})=\overline{\gamma_{W}}(|\bm{x}|)=|W|$. Similar to the
derivation in Eq. (A) the expectation value of $\textsf{DD}(r)$ can be
calculated using the Campbell-Mecke formula:

$\displaystyle\mathbb{E}[\textsf{DD}(r)]$
$\displaystyle=\frac{4\pi}{N^{2}}\int_{r}^{r+\delta}\\!\\!|W|\,\varrho^{2}(1+\xi(s))\,s^{2}\mathrm{d}s$
$\displaystyle=\frac{|S_{r}^{\delta}|}{|W|}+\frac{4\pi}{|W|}\int_{r}^{r+\delta}\\!\\!\xi(s)\,s^{2}\mathrm{d}s$
$\displaystyle\approx\ \frac{4\pi r^{2}\delta}{|W|}(1+\xi(s))\quad\text{for
}\delta\text{ small}.$ (44)

with the volume of the shell
$|S_{r}^{\delta}|=\frac{4\pi}{3}\left((r+\delta)^{3}-r^{3}\right)$. Hence, for
small $\delta$, $\frac{|W|}{|S_{r}^{\delta}|}\textsf{DD}(r)-1$ is an unbiased
estimate of the two-point correlation function $\xi(r)$ in a periodic box –
neither a random points set, nor a boundary correction with geometric factors
is needed.

###  A.2 Expectations of pair-counts for inhomogeneous point sets.

As described in Sect. 6 we are naturally confronted with an inhomogeneously
sampled galaxy distribution. Here we calculate the expectation values of
$\textsf{DD}^{\prime},\textsf{DR}^{\prime}$ and $\textsf{RR}^{\prime}$ for
such an inhomogeneous situation.

One also has Campbell-Mecke formulas for an ingomogeneous point process
$\Phi^{\prime}$ (Baddeley et al., 2000; Shaw et al., 2021):

$\displaystyle\mathbb{E}\left[\sum_{\bm{x}\in\Phi^{\prime}}f(\bm{x})\,\right]$
$\displaystyle=\int_{\mathbb{R}^{3}}f(\bm{x})\,\varrho^{\prime}(\bm{x})\,\mathrm{d}\bm{x},$
(45)
$\displaystyle\mathbb{E}\left[\sum_{\begin{subarray}{c}\bm{x},\bm{y}\in\Phi^{\prime}\\\
\bm{x}\neq\bm{y}\end{subarray}}g(\bm{x},\bm{y})\,\right]$
$\displaystyle=\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}g(\bm{x},\bm{y})\,\varrho_{2}^{\prime}(\bm{x},\bm{y})\,\mathrm{d}\bm{x}\mathrm{d}\bm{y}$
(46)
$\displaystyle=\int_{\mathbb{R}^{3}}\int_{\mathbb{R}^{3}}g(\bm{x},\bm{y})\,\varrho^{\prime}(\bm{x})\varrho^{\prime}(\bm{y})\,(1+\xi^{\prime}(|\bm{x}-\bm{y}|))\,\mathrm{d}\bm{x}\mathrm{d}\bm{y}.$

Again we consider

$\displaystyle S^{\prime}$
$\displaystyle=\sum_{\begin{subarray}{c}\bm{x},\bm{y}\in\Phi^{\prime}\\\
\bm{x}\neq\bm{y}\end{subarray}}\mathbbold{1}_{W}(\bm{x})\mathbbold{1}_{W}(\bm{y})\delta
k_{r}^{\delta}(|\bm{x}-\bm{y}|).$ (47)

and calculate the expectation (see Baddeley et al. 2000 and Shaw et al. 2021
for similar derivations)

$\displaystyle\mathbb{E}[S^{\prime}]=\int_{\mathbb{R}^{3}}\\!\int_{\mathbb{R}^{3}}\\!\\!\mathbbold{1}_{W}(\bm{x})\mathbbold{1}_{W}(\bm{y})\delta
k_{r}^{\delta}(|\bm{x}-\bm{y}|)\varrho_{2}^{\prime}(|\bm{x}-\bm{y}|)\,\mathrm{d}\bm{x}\mathrm{d}\bm{y}$
$\displaystyle=\int_{\mathbb{R}^{3}}\\!\int_{\mathbb{R}^{3}}\\!\\!\mathbbold{1}_{W}(\bm{z}+\bm{y})\mathbbold{1}_{W}(\bm{y})\varrho^{\prime}(\bm{z}+\bm{y})\varrho^{\prime}(\bm{y})\mathrm{d}\bm{y}\
\delta k_{r}^{\delta}(|\bm{z}|)(1+\xi^{\prime}(|\bm{z}|)\mathrm{d}\bm{z}$
$\displaystyle=4\pi\int_{r}^{r+\delta}\\!\\!\overline{\Gamma}_{W}(s)\
(1+\xi^{\prime}(s))\ s^{2}\mathrm{d}s,$ (48)

where $\overline{\Gamma}_{W}(s)$ is the density weighted isotropised set-
covariance

$\overline{\Gamma}_{W}(s)=\frac{1}{4\pi}\int_{0}^{\pi}\\!\\!\\!\int_{0}^{2\pi}\\!\\!\\!\int_{\mathbb{R}^{3}}\mathbbold{1}_{W}(\bm{z}(s,\theta,\phi)+\bm{y})\mathbbold{1}_{W}(\bm{y})\\\
\varrho^{\prime}(\bm{z}(s,\theta,\phi)+\bm{y})\varrho^{\prime}(\bm{y})\
\mathrm{d}\bm{y}\sin(\theta)\mathrm{d}\theta\mathrm{d}\phi$ (49)

For a homogeneous point distribution with $\varrho^{\prime}(\bm{x})=\varrho$
we consistently get
$\overline{\Gamma}_{W}(s)=\varrho^{2}\,\overline{\gamma_{W}}(s)$.

The random point set used in $\textsf{RR}^{\prime}$ is a realisation of an
inhomogeneous Poisson process with number density
$\tfrac{{N_{\mathrm{r}}}}{N}\varrho^{\prime}(\bm{x})$, and a vanishing two-
point correlation function $\xi^{\prime}(r)=0$. $N$ is the number of galaxies
and ${N_{\mathrm{r}}}$ the number of random points in $W$. From Eq. (A.2) we
get

$\displaystyle\mathbb{E}[\textsf{RR}^{\prime}(r)]=\mathcal{RR}^{\prime}(r)$
$\displaystyle=\frac{4\pi}{N^{2}}\int_{r}^{r+\delta}\overline{\Gamma_{W}}(s)\
s^{2}\mathrm{d}s$ (50) $\displaystyle\approx\ \frac{4\pi r^{2}\delta}{N^{2}}\
\overline{\Gamma_{W}}(r)\quad\text{for }\delta\text{ small}.$

Using Eq. (A.2) and Eq. (50) we can determine

$\displaystyle\frac{\mathbb{E}[\textsf{DD}^{\prime}(r)]}{\mathbb{E}[\textsf{RR}^{\prime}(r)]}$
$\displaystyle=\frac{\int_{r}^{r+\delta}\overline{\Gamma_{W}}(s)\
(1+\xi^{\prime}(s))s^{2}\mathrm{d}s}{\int_{r}^{r+\delta}\overline{\Gamma_{W}}(s)\
s^{2}\mathrm{d}s}$ $\displaystyle\approx\ 1+\xi^{\prime}(r)\quad\text{for
}\delta\text{ small},$

Hence we arrive at the well known result that the estimator
$\xi_{\text{PH}}^{\prime}(r)=\frac{\textsf{DD}^{\prime}(r)}{\textsf{RR}^{\prime}(r)}-1$
of Peebles & Hauser (1974) is ‘ratio unbiased’. $\xi_{\text{PH}}^{\prime}(r)$
is an estimate the two-point correlation function $\xi^{\prime}(r)=\xi_{G}(r)$
if we apply the same selection effects to the random point set as we find in
the incompletely sampled galaxy distribution.

For a given point set $\\{\bm{x}_{i}\\}_{i=1}^{N}$ inside the sample geometry
$\bm{x}_{i}\in W$ we calculate the expectation of $\textsf{DR}^{\prime}$.
Using the Campbell-Mecke formula and with $\Phi^{\prime}$ an inhomogeneous
Poisson process with number density
$\tfrac{{N_{\mathrm{r}}}}{N}\varrho^{\prime}(\bm{x})$ we obtain

$\displaystyle\mathbb{E}[\textsf{DR}^{\prime}]$
$\displaystyle=\mathcal{DR}^{\prime}(r)=\frac{1}{N{N_{\mathrm{r}}}}\sum_{i=1}^{N}\int_{\mathbb{R}^{3}}\mathbbold{1}_{W}(\bm{y})\,\delta
k_{r}^{\delta}(|\bm{x}_{i}-\bm{y}|)\,\tfrac{{N_{\mathrm{r}}}}{N}\varrho^{\prime}(\bm{y})\,\mathrm{d}\bm{y}$
$\displaystyle=\frac{1}{N^{2}}\sum_{i=1}^{N}\mathcal{U}_{r}^{\delta}(\bm{x}_{i}),$
(51)

where

$\displaystyle\mathcal{U}_{r}^{\delta}(\bm{x}_{i})=\int_{\mathbb{R}^{3}}\mathbbold{1}_{W}(\bm{y})\,\delta
k_{r}^{\delta}(|\bm{x}_{i}-\bm{y}|)\,\varrho^{\prime}(\bm{y})\,\mathrm{d}\bm{y}$
(52)

is the integral of the density $\varrho^{\prime}(\bm{y})$ in the volume of the
spherical shell
$S_{r}^{\delta}(\bm{x}_{i})=\\{\bm{y}\in\mathbb{R}^{3}\,|\,s<|\bm{y}-\bm{x}_{i}|\leq
s+\delta\\}$ inside $W$.

##  Appendix B Simple windows

The area fraction $\text{area}(\partial B_{r}(\bm{x})\cap W)$ and the
(isotropised) set covariance $\overline{\gamma_{W}}(r)$ can be calculated for
a rectangular box $W$. We use these expressions to determine the accuracy of
the (quasi) Monte-Carlo integration schemes. We also recap the result for a
sphere and for two-dimensions.

###  B.1 Rectangular box

For a point $\bm{x}=(x_{1},x_{2},x_{3})^{T}$ and a rectangular box
$W=[0,L_{1}]\times[0,L_{2}]\times[0,L_{3}]$ with side lengths
$L_{1}>|x_{1}|,L_{2}>|x_{2}|,L_{3}>|x_{3}|$ the set–covariance is

$\gamma_{W}(\bm{x})=(L_{1}-|x_{1}|)(L_{2}-|x_{2}|)(L_{3}-|x_{3}|),$ (53)

and the isotropised set-covariance is for $r<\min\\{L_{1},L_{2},L_{3}\\}$

$\overline{\gamma_{W}}(r)=L_{1}L_{2}L_{3}-\\\
-\frac{r}{2}(L_{1}L_{2}+L_{1}L_{3}+L_{2}L_{3})+\frac{2r^{2}}{3\pi}(L_{1}+L_{2}+L_{3})-\frac{r^{3}}{4\pi}$
(54)

(see e.g. Stoyan & Stoyan 1994; Kerscher 1999). A simple integration gives us
the integrated isotropised set-covariance:

$\int_{r}^{R}\overline{\gamma_{W}}(s)s^{2}\mathrm{d}s=\tfrac{1}{3}L_{1}L_{2}L_{3}(R^{3}-r^{3})-\\\
-\tfrac{1}{8}(L_{1}L_{2}+L_{1}L_{3}+L_{2}L_{3})(R^{4}-r^{4})+\\\
+\tfrac{2}{15\pi}(L_{1}+L_{2}+L_{3})(R^{5}-r^{5})-\tfrac{1}{24\pi}(R^{6}-r^{6}),$
(55)

which we need to calculate $\mathcal{RR}(r)$ according to Eq. (38).

It is more involved to calculate the surface area of a sphere inside a
rectangular box $W$. Using the inclusion-exclusion formula Baddeley et al.
(1993) derive an explicit expression for the area fraction
$\text{area}(\partial B_{r}(\bm{x})\cap W)$:

$\displaystyle\text{area}(\partial B_{r}(\bm{x})\cap W)=4\pi
r^{2}-\sum_{i=1}^{3}\left\\{A_{1}(x_{i},r)+A_{1}(L_{i}-x_{i},r)\right\\}$
$\displaystyle+\sum_{i=1}^{3}\sum_{j=i+1}^{3}\Big{\\{}\,A_{2}(x_{i},x_{j},r)+A_{2}(x_{i},L_{j}-x_{j},r)$
$\displaystyle\qquad\qquad\qquad+A_{2}(L_{i}-x_{i},x_{j},r)+A_{2}(L_{i}-x_{i},L_{j}-x_{j},r)\,\Big{\\}}$
$\displaystyle-
A_{3}(x_{1},x_{2},x_{3},r)-A_{3}(L_{1}-x_{1},x_{2},x_{3},r)-A_{3}(x_{1},L_{2}-x_{2},x_{3},r)$
$\displaystyle-
A_{3}(x_{1},x_{2},L_{3}-x_{3},r)-A_{3}(L_{1}-x_{1},L_{2}-x_{2},x_{3},r)$
$\displaystyle-
A_{3}(L_{1}-x_{1},x_{2},L_{3}-x_{3},r)-A_{3}(x_{1},L_{2}-x_{2},L_{3}-x_{3},r)$
$\displaystyle-A_{3}(L_{1}-x_{1},L_{2}-x_{2},L_{3}-x_{3},r),$ (56)

with

$\displaystyle A_{1}(t_{1},r)$ $\displaystyle=4r^{2}C(t_{1}/r,0,0)$
$\displaystyle A_{2}(t_{1},t_{2},r)$
$\displaystyle=2r^{2}C(t_{1}/r,t_{2}/r,0)$ (57) $\displaystyle
A_{3}(t_{1},t_{2},t_{3},r)$ $\displaystyle=r^{2}C(t_{1}/r,t_{2}/r,t_{3}/r)$

The following expression for $C(a,b,c)$ is almost Eq. (37) from Baddeley et
al. (1993), but with two typos corrected. These correct results have been
already used in the spatstat package for R (Baddeley & Turner 2005, see
sphefrac.c for the implementation). For $a,b,c>0$ and $a^{2}+b^{2}+c^{2}<1$ we
get

$\displaystyle
C(a,b,c)=\and\left(\frac{(1-a^{2}-c^{2})^{1/2}}{ac}\right)+\and\left(\frac{(1-b^{2}-c^{2})^{1/2}}{bc}\right)$
$\displaystyle+\and\left(\frac{(1-a^{2}-b^{2})^{1/2}}{ab}\right)-a\and\left(\frac{(1-a^{2}-c^{2})^{1/2}}{c}\right)$
$\displaystyle+a\and\left(\frac{b}{(1-a^{2}-b^{2})^{1/2}}\right)-b\and\left(\frac{(1-b^{2}-c^{2})^{1/2}}{c}\right)$
$\displaystyle+b\and\left(\frac{a}{(1-a^{2}-b^{2})^{1/2}}\right)-c\and\left(\frac{(1-a^{2}-c^{2})^{1/2}}{a}\right)$
$\displaystyle+c\and\left(\frac{b}{(1-b^{2}-c^{2})^{1/2}}\right)-\pi.$ (58)

For $a^{2}+b^{2}+c^{2}\geq 1$ we have $C(a,b,c)=0$. As special cases

$C(a,b,0)=\and\left(\frac{(1-a^{2}-b^{2})^{1/2}}{ab}\right)-a\and\left(\frac{(1-a^{2}-b^{2})^{1/2}}{b}\right)\\\
-b\and\left(\frac{(1-a^{2}-b^{2})^{1/2}}{a}\right),$ (59)

and $C(a,0,0)=\tfrac{\pi}{2}(1-a)$.

###  B.2 Sphere

For a spherical sample $W=B_{R}$ centred at the origin with radius $R$ and
with $r=|\bm{x}|<R$ the set-covariance reads

$\gamma_{W}(\bm{x})=\frac{4\pi}{3}\left(R^{3}-\frac{3}{4}rR^{2}+\frac{1}{16}r^{3}\right).$
(60)

Due to symmetry, the isotropised set-covariance is simply
$\overline{\gamma_{W}}(|\bm{x}|)=\gamma_{W}(\bm{x})$. The area fraction for a
point $\bm{q}\in W$ at a distance $|\bm{q}|=s>R-r$ from the origin, and with
$r<R$ is $\text{area}(\partial B_{r}(\bm{q})\cap W)=4\,\alpha\,r^{2}$, with
$\cos(\alpha)=\frac{r^{2}+s^{2}-R^{2}}{2rs}$.

###  B.3 Two-dimensions

Here we give the expressions for the set-covariance in two-dimensions. For a
point $\bm{x}=(x_{1},x_{2})^{T}\in\mathbb{R}^{2}$ and a rectangle
$W=[0,L_{1}]\times[0,L_{2}]$ with side lengths $L_{1}>|x_{1}|,L_{2}>|x_{2}|$
the set–covariance is

$\gamma_{W}(\bm{x})=(L_{1}-|x_{1}|)(L_{2}-|x_{2}|),$ (61)

and for $r<\min\\{L_{1},L_{2}\\}$ the isotropised set-covariance is (Ripley,
1988)

$\displaystyle\overline{\gamma_{W}}(r)$
$\displaystyle=\frac{1}{2\pi}\int_{0}^{2\pi}\gamma_{W}(\bm{x}(r,\phi))\,\mathrm{d}\phi$
$\displaystyle=L_{1}L_{2}-\frac{2r}{\pi}(L_{1}+L_{2})+\frac{r^{2}}{\pi}.$ (62)

The analogue to $\text{area}(\partial B_{r}(\bm{x})\cap W)$ is the arc-length
$\text{arc}(C_{r}(\bm{x})\cap W)$ of the circumference of a circle
$C_{r}(\bm{x})$ with radius $r$ centred on a point $\bm{x}\in W$ inside the
sample window $W$. It can be calculated from the intersection points of the
circle with the sides of $W$ and the angles between them and the centre
(Stoyan et al., 1995).

##  Appendix C Details on the implementation

We use Python with NumPy (Harris et al., 2020) and SciPy (Virtanen et al.,
2020) to implement the data processing, to calculate geometric quantities, to
generate random points, and to generate randomised Halton sequences (provided
in SciPy starting with version $\geq$1.7.0). The performance critical parts
are the pair-counting and a function for checking whether points are inside
the observational window. We implemented these routines as C++ functions
parallelised using OpenMP (Dagum & Menon, 1998). They are made callable from
python via pybind11 (Jakob et al., 2017). In the computations of DD and RR we
obviously avoid the double counting by using $\sum_{i}\sum_{j\neq
i}\ldots=2\sum_{i}\sum_{j>i}\ldots$. We implement the double sums without
using a tree, resulting in a code with a runtime quadratic in the number of
points. Using OpenMP we could achieve a reduction of the computing time almost
linear with the number of computing cores used (at least up to the 20 cores we
could easily muster, see also Alonso 2012). We tried the kd-tree provided with
SciPy but in our application with rather large radii our parallelised code
outperformed the single core cKDTree. This can be different for other
applications, especially if one is considering small radii only.

Only these basic optimisation are included in our code but perhaps it is
versatile enough to help with future implementations. You find a basic version
of the code at https://homepages.physik.uni-
muenchen.de/\~Martin.Kerscher/software/accuratexi/