{'fraction_non_alphanumeric': 0.030614379084967322, 'fraction_numerical': 0.011058823529411765, 'mean_word_length': 4.658431952662722, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A valid and fast spatial bootstrap for correlation functions

In this paper, we examine the validity of non-parametric spatial bootstrap as a procedure to quantify errors in estimates of N-point correlation functions. We do this by means of a small simulation study with simple point process models and estimating the two-point correlation functions and their errors. The coverage of confidence intervals obtained using bootstrap is compared with those obtained from assuming Poisson errors. The bootstrap procedure considered here is adapted for use with spatial (i.e. dependent) data. In particular, we describe a marked point bootstrap where, instead of resampling points or blocks of points, we resample marks assigned to the data points. These marks are numerical values that are based on the statistic of interest. We describe how the marks are defined for the two- and three-point correlation functions. By resampling marks, the bootstrap samples retain more of the dependence structure present in the data. Furthermore, this method of bootstrap can be performed much quicker than some other bootstrap methods for spatial data, making it a more practical method with large datasets. We find that with clustered point datasets, confidence intervals obtained using the marked point bootstrap has empirical coverage closer to the nominal level than the confidence intervals obtained using Poisson errors. The bootstrap errors were also found to be closer to the true errors for the clustered point datasets.


Introduction
In analyses of survey data, such as those of galaxies or quasars, N-point correlation functions are often estimated (e.g. Kulkarni et al. 2007;McCracken et al. 2007;Shen et al. 2007). These help to describe the structure of the observed objects, such as the filamentary structure of galaxies that has been observed, and to constrain the parameters of cosmological models. Such estimates are only as important as their associated errors, since it is the errors that indicate the amount of agreement between two sets of data or between data and simulations from a model.
In spatial point processes, the expressions for the standard errors of correlation (and similar) functions have been worked out only for the simplest of models. For example, Ripley (1988) found approximations for the variance of the K function, an integral of the two-point correlation function, for the Poisson process. These depend on such factors as the shape of the observation region and the type of correction method for boundary effects. Landy & Szalay (1993), Hamilton (1993) and others worked out approximations for standard errors of various estimators of the two-point correlation function under the Poisson or weakly clustered models. Thus, very often, approximations such as Poisson errors are used instead. However, point data arising in astronomy are typically clustered and non-Poisson. So while Poisson errors are useful and easy to compute, they only serve as rough indications of the size of errors.
Besides using Poisson errors, errors can also be estimated by using mock catalogs generated from a cosmological model. This method was employed in Eisenstein et al. (2005), where the initial conditions for the cosmological model were selected independently. In the statistics literature, this is referred to as parametric bootstrap, although the term more commonly refers to mock datasets generated from the model using parameter values fixed at the estimates from the data, instead of being independently chosen.
An alternative is to use non-parametric bootstrap. This involves generating new samples, called bootstrap samples, by resampling from the actual data, and computing estimates for these new samples. The distribution of these bootstrap estimates then serves as a proxy for the actual distribution of the data estimates, so that statistical inference, such as the construction of confidence intervals, can be performed. Note that the procedure does not make any specific model assumptions, thus the errors obtained by this method can serve as a check of model assumptions. Due to the simplicity and flexibility of the non-parametric bootstrap, the method is attractive. What is desirable then is to make the non-parametric bootstrap procedure work as well as possible for data that is correlated, and check that it performs satisfactorily, so that it can be useful as a tool in analysis. This paper thus examines the non-parametric bootstrap, specifically bootstrap of spatial data, where the dependence present in the data is of interest. There were some early misconceptions about how bootstrap should be applied to spatial data. The naive method of resampling individual points does not work in the spatial context. In order for spatial bootstrap to be valid, the underlying dependence structure has to be preserved as much as possible when generating bootstrap samples. Two common methods for doing this are the block bootstrap and subsampling where blocks of data, instead of individual data points, are resampled. We introduce these methods in Section 2 and describe their shortcomings. In Section 3 we describe the marked point bootstrap (Loh & Stein 2004) as a way to address these shortcomings. We describe how the marked point bootstrap can be used with the two-and three-point correlation function estimators, and by extension to estimators of Npoint correlation functions. In Section 4 we present results of a simulation study using simple point process models comparing the empirical coverage of confidence intervals obtained using non-parametric bootstrap and using normal approximations with Poisson errors.
In this paper, we restrict ourselves to constructing nominal 95% confidence intervals, i.e. these confidence intervals are supposed to contain the true value 95% of the time. The empirical coverage of the confidence intervals is the actual confidence level achieved by the confidence intervals. In a simulation study with a known model, the empirical coverage can be obtained by finding the number of confidence intervals that contain the true value and then compared with the nominal level. It is desirable, of course, for the empirical coverage to be close to the nominal level. Furthermore, it is often better for the empirical coverage to be higher instead of lower than the nominal level, so that the procedure is conservative.
Bootstrap is a computationally intensive procedure. With the large datasets now common in astronomy, even computing the N-point correlation functions pose computational challenges. For example, Eisenstein et al. (2005) avoided using the jackknife procedure for error estimation because of the size of the data they used. The way the marked point bootstrap is formulated, however, makes it much faster than subsampling (a generalization of the jackknife) and the block bootstrap, so that applying the procedure to large datasets is feasible as long as computing the actual estimates is feasible. In Section 4 we provide some time measurements of the procedure used in our simulation study.

Non-parametric bootstrap for spatial data
The non-parametric bootstrap was originally developed for independent data (Efron & Tibshirani 1994). The main idea is to draw new samples from the actual data by sampling with replacement a data point at a time. Bootstrap estimates of the same statistic are computed from the bootstrap samples. With these bootstrap estimates, confidence intervals, for example, can then be constructed. This can be done in a variety of ways. Suppose K,K and K * i , i = 1, . . . B are respectively the quantity of interest, the estimate of K computed from the data and the bootstrap estimates, with B equal to the number of bootstrap samples. A simple method, called the basic bootstrap interval in Davison & Hinkley (1997), is to set as the 100(1 − α)% confidence interval for K. Here, B, the number of bootstrap samples, is large, say, 999, andK * u is the u-th ordered values of the bootstrap statistic. So, for example, with B = 999 bootstrap samples, a 95% confidence interval for K is given by In our simulation studies, we use (1) to construct the confidence intervals. Standard errors forK are estimated by the standard deviation of the bootstrap estimatesK * .
While there are other methods of constructing confidence intervals from bootstrap samples (see Davison & Hinkley 1997, for example), the interest here is in the method of generating the bootstrap samples, when the data is spatial. Snethlage (1999) rightly concludes that resampling individual points do not work. If the resampled points are placed in their original positions in the observation region, there will be multiple points at single locations, which do not usually occur in most data sets. In their claim that bootstrap cannot be used for analysis for clustering, Simpson & Mayer-Hasselwander (1986) were also considering bootstrap in terms of resampling individual points.
Due to the success of bootstrap for resampling independent data, it has been extended to resample dependent data. Most of this work is for time series, but can easily be applied to spatial data in two and three dimensions. A common method is the block bootstrap: blocks of the spatial data are sampled at random, then joined together to form a new sample (Hall 1985;Künsch 1989;Liu & Singh 1992). Asymptotic arguments for the validity of the bootstrap involve limiting the range of dependence, increasing the observation region size and letting the resampling block size increase but at a slower rate than the observation region size. In this asymptotic setting, we then have many almost independent blocks of data, with each block itself containing a large subsample (see e.g. Lahiri 2003). However, the assumed conditions necessary to make the calculations tractable also means that normal approximations work well too. Some theoretical results show that the accuracy of bootstrap estimates is of a higher order than the normal approximations. Whether this difference is meaningful in actual practice is less clear. We believe that the role of non-parametric bootstrap is to serve as another objective method to obtain standard errors that do not make any model assumptions. Error estimates obtained using bootstrap can be used as a way to assess or compare with other estimates of errors. Kemball & Martinsek (2005) is a recent work in the astronomy literature that examined bootstrap for dependent data. For the non-parametric bootstrap, they focused on subsampling (Politis & Romano 1993;Politis et al. 1999), which can be considered a generalization of the jackknife procedure. In subsampling, random portions of the data are deleted, and the remaining data are treated as bootstrap samples. The standard deviation of the estimates computed from these samples serves as an estimate of the error, less a factor to adjust for the smaller subsamples and the large overlap.
When estimating correlation functions, pairs or triplets etc of points have to be counted. By joining independently resampled blocks together to form the bootstrap sample, the block bootstrap creates artifical configurations of points across the resampling blocks and distorts the dependence structure in the data. This does not matter in asymptotic arguments because the effect becomes negligible if the range of the correlation is fixed while the resampling blocks increase in size. However, Loh & Stein (2004) found that the actual coverage achieved by confidence intervals obtained using block bootstrap can be much lower than the nominal percentage level for finite samples.
In subsampling, no artificial configurations of points are created. However, while the correction weight accounts for the difference in sample sizes between the bootstrap samples and the actual data set, it does not account for the change in the boundary effects due to the different resampling regions. Since subsampling uses smaller regions as the bootstrap observation regions, boundary effects are magnified. For subsampling, there is the temptation to use large subsamples to try and retain more of the dependence structure, but like block bootstrap, theoretical justification of the method requires that the subsamples be small in size relative to the actual data set. Loh & Stein (2004) also found that subsampling can yield confidence intervals that attain very low empirical coverage. They also found that the subsampling method is sensitive to the fraction of the data used for subsampling. Loh & Stein (2004) proposed another version of spatial bootstrap, called marked point bootstrap, that reduces the effect of joining independent blocks and produces confidence intervals that achieve coverage closer to the nominal level. This is described in the next section, where we also show how it can be applied to the two-and three-point correlation function estimators commonly used in astronomy.

Improving the non-parametric bootstrap of spatial data
Suppose N points are observed in a region A. Furthermore, suppose that the quantity of interest K can be estimated using an estimator of the form Note that each point i has an associated quantity N j=1,j =i φ(x i , x j ), the inner sum of equation (2). Estimators of two-point statistics can be expressed in this form. In this case, the quantity φ(x i , x j ) will depend on the distance between x i and x j . As an aside, note that estimators of three-point statistics can be written in a similar form, with the inner sum replaced by a double sum.
With point data, the term "mark" is used to refer to some additional information associated with a point. This is usually some actual measured value. For galaxy data, for example, marks could be quantities such as luminosity, color and so on. In this paper, the bootstrap method considered uses marks associated with the points. However, these marks are not quantities such as luminosity that are directly measured. Instead they are numerical quantities that we construct and associate with the points. The actual values of these marks are not random, but are constructed so that they relate to the statistic that is of interest. If the statistic of interest is given by equation (2), then the mark associated with point At the risk of being repetitive, suppose thatΦ ≡ DD(r) = x∈D y∈D:y =x 1{|x − y| ∈ (r − dr, r + dr)}, for some r, is of interest. This quantity is used in estimators of the two-point correlation function, and is the number of pairs of points separated by (roughly) distance r. Then the mark associated with point x is y∈D:y =x 1{|x − y| ∈ (r − dr, r + dr)}, the number of points that are roughly distance r away from x. Note that the sum of all the marks gives back the value of DD(r). It is also important to note that to compute the estimate (2), the marks have to be calculated anyway. In regular applications, the algorithm doing the estimation does not individually record these marks, but keeps a running sum of the marks. In order to do the marked point bootstrap, the difference in terms of the code is that the marks now have to be stored so that they can be used in the bootstrap step.
In the block bootstrap, blocks of data points are resampled and then joined together, forming a new dataset from which K is estimated using the new configuration of points that was generated, yieldingK * . In the marked point bootstrap, blocks can be used to resample points as well. However, the crucial difference is that the bootstrap estimate is computed, not from how the resampled points are positioned, but from the marks that are associated with these points. In other words, the marked point bootstrap resamples the marks rather than the points and the bootstrap estimate is computed by summing these resampled marks.
To be more precise, suppose that N * number of points have been resampled, with the resampled points denoted by x * j , j = 1, . . . N * . Associated with each x * j is a mark m j * . We denote this mark m j * rather than m * j to emphasize the fact that these marks are sampled from the actual data, i.e. computed from the original dataset and not from the bootstrap sample. Then the bootstrap estimate of K is given by the average of the resampled marks: just likeK is given by the average of the actual marks. Note that in an actual implementation of the procedure, all that is required is keeping track of how many times each point is resampled. The step-by-step procedure for estimating and resampling the quantity (2) is as follows: 2. Obtain the estimateK, usingK = i m i /N.
3. Resample the points. This can be done by randomly placing blocks on to the observation region and keeping track of which point is resampled. Suppose point i, i = 1, . . . , N is resampled n * i times, and N * = i n * i .

The bootstrap estimate is thenK
Repeat steps 3 and 4 to get B bootstrap estimates.
A few remarks about the procedure are in order.
Remark 1 Instead of randomly placing blocks, the observation region can be divided into a number of subregions, and the regions selected randomly with replacement. This latter method is sometimes referred to as using fixed blocks as opposed to moving blocks. It is generally considered that the moving blocks bootstrap works better in terms of convergence rates in asymptotic arguments.
Remark 2 The number of blocks used is so that the total area/volume of the blocks is equal to the original area/volume of the observation region. Note that in this case N * would usually not be equal to N, though they will be of the same order of magnitude. However, this does not pose problems since the statisticK is a mean of the marks.
Remark 3 There is no real consenus on the size of the resampling blocks to use. Bühlmann & Künsch (1999) did some work on determining the optimal block size from data. Intuitively, the procedure needs large blocks so that the correlation structure is less distorted, and a large enough number of blocks so that there is enough variability between bootstrap samples. If K represents the number of blocks and N the number of data points (which is assumed to increase with the observation region size), theoretical work in e.g. Lahiri (2003) suggests that consistency is achieved as K → ∞ and N/K → ∞. Thus some trade-off is needed. A rule-of-thumb is to divide each dimension of the observation region into at least three parts, i.e. nine blocks in 2D, 27 blocks in 3D. This would ensure enough variability between bootstrap samples. Of course, for correlation functions, the maximum value of the separation distance r at which these functions are estimated would influence the decision on block size. Fortunately, Loh & Stein (2004) found that the marked point bootstrap is less sensitive to block size than block bootstrap or subsampling: they resampled an absorber catalog using slices of the sphere and found that the bootstrap errors were similar for different sizes of the slices. Our simulation results also show little difference due to different block size.
There are a few advantages to this form of spatial bootstrap over the regular block bootstrap. Since the bootstrap estimates are based on the resampled marks and not on marks recalculated from the bootstrap sample, the contribution to the bootstrap estimate is due to actual pairs of points in the original dataset. This helps to minimize the distortion of dependence structure in the dataset due to resampling.
Furthermore, for any block of resampled points, information about the points just outside the block (and therefore not sampled by this particular block) is captured by the marks associated with the points that are sampled by the block. This helps to reduce the variability of bootstrap results due to the size of the resampling blocks, compared to block bootstrap or subsampling. Also, since the resampling blocks do not need to be joined together to form a contiguous region for the bootstrap sample, there is flexibility in the choice of the shape of the resampling regions.
Lastly, the marked point bootstrap can be performed relatively quickly compared to block bootstrap or subsampling. The marks that are associated with the points are part of the actual estimator and are already computed in the estimation step. Resampling using the marked point bootstrap only involves identifying which points are resampled with each resampling region, and keeping track of how many times each point is resampled. Interpoint distances and edge correction weights do not have to be recalculated. With N data points and B bootstrap samples, block bootstrap will take roughly BN 2 computations for a statistic involving pairs of points. The marked point bootstrap will involve roughly N 2 +BN computations. The difference will be more marked for three-point computations.
Simulation studies done in Loh & Stein (2004) showed that the empirical coverages of confidence intervals obtained using the marked point bootstrap can be much closer to the nominal 95% level than those obtained with block bootstrap or subsampling.
We now describe how the marked point bootstrap can be used with estimators of the two-point correlation function. The common estimators of the two-point correlation function which are, respectively, the natural estimator (Kerscher et al. 2000), and estimators due to Davis & Peebles (1983); Hamilton (1993); Landy & Szalay (1993); Hewett (1982), where r is some distance of interest. In these expressions, dd(r) = DD(r)/N 2 , dr(r) = DR(r)/NN R and rr(r) = RR(r)/N 2 R , where DD(r) = x∈D y∈D:y =x 1{|x − y| ∈ (r − dr, r + dr)}/N 2 , DR(r) = x∈D y∈R 1{|x−y| ∈ (r −dr, r +dr)}/NN R and RR(r) = x∈R y∈R:y =x 1{|x− y| ∈ (r − dr, r + dr)}/N 2 R , R is a set of randomly generated points (i.e. Poisson) in the observation region A, and N and N R are respectively the number of points in the real and random data sets.
To apply the marked point bootstrap, assign to each point x of the dataset marks m x,1 = y∈D:y =x 1{|x − y| ∈ (r − dr, r + dr)} and m x,2 = y∈R 1{|x − y| ∈ (r − dr, r + dr)}. Bootstrap proceeds by resampling blocks of points and recording the marks associated with them. For a bootstrap sample, x * j , j = 1, . . . N * , we then have and bootstrap estimates of the two-point correlation functions are then obtained by substituting the above into (3)-(7). If each point x i of the actual data is resampled n * i times, so that N * = i n * i , DD * (r) and DR * (r) can also be written as Note that RR does not need to be resampled, since it is used as an approximation to an integral and has nothing to do with the actual data. If, as is usually the case, estimation of ξ(r) is needed for a range of values of r, then the marks m x,1 and m x,2 would be vectors, containing the relevant values for each value of r.
Estimators of the three-point correlation function can be bootstrapped in a similar way. For example, an estimator of the three-point correlation function is introduced by Peebles & Groth (1975), where ddd = DDD/N 3 , ddr = DDR/N 2 N R and rrr = RRR/N 3 R and DDD, DDR, RRR are counts of triplets of points with the desired configuration, DDD with all points from the real data set and so on. The contribution to DDD by any particular triplet of points is divided by 3 and assigned as marks to each of the three points. For any individual point, all these marks are summed together. For DDR, the contribution by each triplet is divided by 2 and assigned to the two real data points. Bootstrap proceeds by resampling the real data points and the values of DDD * and DDR * found by adding the marks of the resampled points. Substituting these into (8) gives the bootstrap estimate. Other similar estimators, such as the three-point estimator of Jing & Börner (1998) or the N-point estimators of Szapudi & Szalay (1998), can be bootstrapped in the same way.

Simulation study
We performed a simple simulation study to compare the performance of confidence intervals obtained using the marked point bootstrap with those obtained using normal approximations with Poisson errors, varying the observation region size, number density and point process model. For computational simplicity, we restrict to two dimensions. We also performed an additional study with a large observation region and approximately 50,000 points, showing the applicability of the marked point bootstrap to datasets of size comparable to current astronomy datasets.
We used the Poisson point process model and a Neyman-Scott model to generate the data points. The Neyman-Scott model is of historical interest in astronomy as a model for galaxies (Neyman & Scott 1952). It is still commonly used to model point data in other fields (Diggle 2003;Waagepetersen 2007). We chose the Neyman-Scott model as it is a model for clustered data with closed-form expressions for the two-point correlation function. The Neyman-Scott point datasets that we used are generated as follows: parent points are distributed as a Poisson point process with intensity λ p . A Poisson number with mean m of offspring points are then randomly scattered about each parent point. The collection of offspring points form the point process. We set the dispersion function of offspring points about parent points to be a bivariate normal density centered at the parent point, with standard deviation σ. This specific Neyman-Scott model is sometimes referred to as the modified Thomas model (Stoyan et al. 1995). The two-point correlation function, ξ(r), is zero for the Poisson model, while for the modified Thomas model. Thus the point pattern from a modified Thomas model is more clustered if λ p or σ is smaller. The quantity σ also controls the range of the correlation, with the range larger for larger values of σ. We used several values for λ p , m and σ in our simulation study.
For each point process model, we generated 500 realizations on the unit square. For each realization, we estimated ξ(r) for r = 0.01, . . . , 0.1. Bootstrap estimates were then produced from each realization and a nominal 95% confidence interval constructed. Thus for each point process model, we have 500 95% confidence intervals. We then checked the the empirical coverage, i.e. the proportion of these that contained the true value of ξ(r), with proportion closer to 95% being desirable. We also constructed 500 confidence intervals using the normal approximation with Poisson errors. The Poisson error e p is the inverse of the pair counts for an uncorrelated data set of the same size as the actual data, as given by Landy & Szalay (1993). The 95% confidence intervals for ξ based on the normal approximation are thus given by (ξ ± 2e p ). We then found the empirical coverage of these confidence intervals. We then repeated the procedure for the 2 × 2 and 4 × 4 squares. The results are summarized in Figures 1 to 3. Figure 1 shows plots of the empirical coverage of nominal 95% confidence intervals of the two-point correlation function for the Poisson process model, using the Hamilton (1993) estimator. Simulation results for the other estimators are similar and are not shown. The thick solid lines in the plots show the empirical coverage of confidence intervals obtained using normal approximation with Poisson errors. Note that Poisson errors are correct in this case and we find that the empirical coverage is close to 95% for all the point densities and observation region sizes considered.
The thin lines represent the empirical coverage of confidence intervals obtained from the marked point bootstrap, with the different line types representing different resampling block sizes. These were squares of lengths 0.5, 0.33 and 0.25 for the 1 by 1 regions, of lengths 1, 0.67, 0.5, 0.33 for the 2 by 2 regions and of lengths 2, 1, 0.67, 0.5 for the 4 by 4 regions (solid, dashed, dotted and dashed-dotted lines respectively for increasingly smaller blocks). The difference due to the block size used for resampling appear to be small. As mentioned, this was an advantage of the marked point bootstrap. Loh & Stein (2004) found greater variation of performance with block size for subsampling and block bootstrap.
Compared with the Poisson empirical coverage, we find that at low densities and smaller observation region sizes (plots towards the upper left of Figure 1), the bootstrap method does poorly. However, the empirical coverage of the bootstrap confidence intervals quickly increases towards 95% with increasing density (down the columns in Figure 1) and/or observation region size (across the rows in Figure 1), i.e. with larger sample sizes.  Hamilton (1993). Confidence intervals are obtained using normal approximation with Poisson errors (thick solid line) and with the marked point bootstrap using different resampling block sizes. The block sizes were squares of lengths 0.5, 0.33 and 0.25 for the 1 by 1 regions, of lengths 1, 0.67, 0.5, 0.33 for the 2 by 2 regions and of lengths 2, 1, 0.67, 0.5 for the 4 by 4 regions (solid, dashed, dotted and dashed-dotted lines respectively for increasingly smaller blocks).  -Plots of the empirical coverage of nominal 95% confidence intervals of the twopoint correlation function for the modified Thomas process model for realizations in a 2 × 2 square. The estimator used is that of Hamilton (1993). Confidence intervals are obtained using normal approximation with Poisson errors (thick solid line) and with the marked point bootstrap using different resampling block sizes (see text).   We find that when the point pattern is only weakly clustered (left plot, for the case λ p = 100, m = 2.5 and σ = 0.15), the Poisson confidence intervals had empirical coverage close to the nominal 95% level. However, as the other two plots in Figure 2 show, the empirical coverages of the Poisson confidence intervals become lower than 95% as the degree and/or range of clustering increases (i.e. with smaller σ or λ p ). On the other hand, the boostrap confidence intervals attain coverage much closer to 95% for all the cases shown, regardless of the degree of clustering. Plots of the Poisson and bootstrap error estimates are shown in Figure 3. Notice that the Poisson approximation underestimates the true error as the degree of clustering increases, while the bootstrap error estimates remain close to the true errors, even for the modified Thomas model.
Thus we find that the performance of the Poisson confidence intervals is sensitive to the degree of clustering of the point pattern. If the point pattern is Poisson, or weakly clustered, the empirical coverage of Poisson confidence intervals is close to the nominal level, even with small sample sizes. However, performance quickly deteriorates with greater degree of clustering. On the other hand, the bootstrap confidence interval does not perform well with small sample sizes. With moderate sample sizes, however, the bootstrap method performs rather well, over a wide range in the degree of clustering.
We performed an additional set of simulations using data sets of roughly 50,000 points in a 20 × 20 square and estimating ξ(r) for r = 0.01 to 2. Other than the restriction to 2D, the data size and range of r is roughly of the scale found in current astronomy data. We used the modified Thomas model and chose four sets of parameter values, varying the degree and range of clustering but with the same number density. A sample realization from each of the four models is shown in Figure 4. The models corresponding to the top row in Figure 4 have higher clustering than the models on the bottom row. For models in the same row, the strength of clustering is similar, but the model on the right has a longer correlation range. We used square resampling blocks of side length 5, 2.5 and 2 and results were very similar.
The results are summarized in Figure 5, which show the empirical coverage of confidence intervals (left column) and errors (right column) obtained from the marked point bootstrap and with Poisson errors, for each of the four Thomas models. The plots qualitatively show the same relative performance between Poisson errors and bootstrap as found in the earlier simulation study. When the range of clustering is large, the empirical coverage of confidence intervals based on Poisson errors and the normal approximation is very low (second and fourth plots on the left column of Figure 5). The coverage of the bootstrap intervals are affected too, but by much less. When the correlation is large, the Poisson errors substantially under-estimate the true errors, while the marked bootstrap errors were more realistic. At the larger values of r, especially when ξ is near 0, the bootstrap procedure appears to be somewhat conservative, while the Poisson errors become more accurate.
We made some time measurements of various sections of the algorithm: the functions computing DD and DR took 1 minute and 13 minutes respectively. Here, N R = 200, 000 and we did not use any sophisticated methods (such as tree-based algorithms) to speed up the computation. The bootstrap function, generating 999 samples and computing the estimates, took roughly 1 minute, showing the feasibility of the procedure for large data sets. The speed of the marked point bootstrap is due to the fact that the marks that are resampled have already been computed as part of the estimation. The additional computational burden of the bootstrap is sampling the points and keeping track of the number of times each point is resampled.

Discussion
In this paper, we introduced the marked point bootstrap as a method to bootstrap spatial data for estimating errors without specific model assumptions. In particular, we described how the method can be applied to estimators of the two-and three-point correlation functions. With the non-parametric bootstrap, errors are obtained from the actual data. There is no need choose a model, select parameter values or generate mock catalogs using N-body simulations. Thus errors obtained from non-parametric bootstrap can be used to compare with errors obtained from other methods with more specific model assumptions.
For non-parametric spatial bootstrap, we propose the marked point bootstrap over the more common block bootstrap or subsampling methods. There are several advantages of the marked point bootstrap. Firstly, by using information from actual pairs or triplets of points in the data, bootstrap confidence intervals using the marked point bootstrap attain better empirical coverage than confidence intervals constructed using block bootstrap (see Loh & Stein 2004 for a comparison of these two methods).
Secondly, in the marked point bootstrap, it is the marks that are used to compute the bootstrap estimate. Thus the resampled points do not have to be arranged to form a new point pattern. This makes it a lot easier to bootstrap data that are observed in irregularly shaped regions that are common in astronomy. In Loh & Stein (2004) for example, bootstrap on an absorber catalog was done using slices as well as spheres, with similar results for both types of resampling regions.
Thirdly, the marks used for resampling are part of the original estimate and are computed during the estimation step. The only additional computation required by the marked point bootstrap involves selecting points (that is, testing whether each point lies in a resampling region or not), and keeping track of the number of times each point is resampled. Unlike the block bootstrap, there is no need to re-compute from scratch the estimates for each bootstrap sample. This difference in computation is even greater for higher-order statistics. These properties of the marked point bootstrap make it a computationally feasible tool for analysis.
Our study here suggests that non-parametric bootstrap can yield valid estimates of errors under a wide range of point patterns. The lack of specific model assumptions means that the non-parametric bootstrap method, and in particular the marked point bootstrap, can serve as an alternate and complementary method for quantifying errors. Having estimates of errors obtained using Poisson approximations, parametric and non-parametric bootstrap allows one to have a better sense of the size of errors involved in an analysis.
The simulation study performed here shows that bootstrap confidence intervals do attain coverage close to the nominal level, even for the clustered point patterns where Poisson errors are known to be inaccurate, when sample sizes are large. More specifically, bootstrap performance improves with increasing number density, and also with increasing observation region size relative to the correlation length. Unfortunately, in astronomy, the correlation length may be of the same scale as the observation region. If the values of r at which the correlation function estimates are computed are small relative to the resampling blocks (and the observation region), then although the bootstrap procedure would distort the dependence structure at the large scales, it would still be valid for these smaller values of r.
If, instead, ξ(r), say, for r close to the size of the observation region is of interest, then the bootstrap procedure would start to break down, in the sense that the empirical coverage of confidence may not be close to the nominal level, and the bootstrap errors not reflect the true errors. In this case, the amount of information contained in the data is smaller and the boundary effects are magnified. With respect to the marked point bootstrap, larger blocks are needed to capture the dependence structure at this larger scale. For a fixed sample, this cannot be done without reducing the variability of the bootstrap samples. The method that might work best is parametric bootstrap, assuming that the model is correct, and that the parameter values used are close to the true values. Non-parametric bootstrap can still be useful here. Firstly, it is at least a better choice than Poisson errors, since the latter would grossly underestimate the true errors. Secondly, it can provide additional error estimates to compare with the errors obtained with the assumed model. For these most challenging instances, having a variety of methods can only be beneficial.
This research is supported in part by National Science Foundation award AST-0507687.