# Biased-estimations of the Variance and Skewness

## 1 Introduction

There has been a long history of interests, since the pioneering work of Peebles (Peebles (1980), §18) in the hierarchical amplitudes $`S_N`$, defined as the following ratio:

$$S_N=\overline{\xi }_N/\overline{\xi }_2^{N1},$$
(1)   
where $`\overline{\xi }_N`$ is the N-cumulant defined by $`\delta ^N_c`$, and $`\delta `$ is the density fluctuation smoothed on some scale. These quantities are important as a test of the gravitational instability paradigm (e.g. Fry (1984); Juszkiewicz et al. (1993); Bernardeau (1994)), a probe of possibly non-Gaussian initial conditions (e.g. Silk & Juszkiewicz (1991); Gaztañaga & Maehoenen (1996); Gaztañaga & Fosalba (1997)) as well as a measure of the galaxy-bias (e.g. Fry & Gaztañaga (1993); Frieman & Gaztañaga (1994)). However, it has also been a puzzle for quite some time that different galaxy surveys yield discordant values for $`S_N`$ (e.g. Table 1, for $`N=3,4`$). While some of the differences no doubt arise from the fact that galaxies selected in different ways might have different galaxy-biases, not all of the differences can be convincingly explained away in such a manner. For instance, a comparison between optically selected galaxy-catalogs (in Table 1) reveal a substantial and systematic difference between the measured values of $`S_N`$: redshift surveys consistently yield lower values compared to the larger angular catalogues (e.g. compare APM/LICK/EDSGC values with those from CfA/SSRS in Table 1; note that the IRAS galaxies are infrared-selected; note also an exception to this rule in the measurements by Kim & Strauss (1998)). A rather large relative galaxy-bias between these two sets of catalogues would have to be invoked to reconcile them.

Three alternative explanations are possible. The first is that redshift space distortions tend to suppress $`S_N`$, but it has been shown to be not sufficient to explain the systematic differences, especially on large scales (Fry & Gaztañaga (1994)). The second is that the local volume (sampled by the redshift surveys) just happens to have a smaller $`S_N`$ compared to the true $`S_N`$ which is presumably measured in the larger angular surveys i.e. the local universe is not a fair sample (e.g. Gaztañaga (1994)). This is related to the question of the homogeneity scale of our universe which has been a subject of some debate (see e.g. Wu et al. (1998) ). The third is to blame it on the estimator for $`S_N`$: it yields a value that is on the average biased low, with the bias (distinguish this from the galaxy-bias) getting worse as the survey becomes smaller. We will demonstrate that the third contributes to a significant fraction of the systematic differences between the different measurements. While a thorough analysis detailing exactly how much each of these factors contribute is beyond the scope of this paper, we can safely conclude that inference of large sampling fluctuations, or a large relative galaxy-bias, based on measurements of $`S_3`$, are unwarranted.

A very closely related estimation-bias has been discussed before by Colombi et al. (Colombi et al. (1994)) as a finite-volume effect.<sup>1</sup><sup>1</sup>1Related ideas have also been considered by Bromley 1998, private communication. They attributed this to the abrupt cut-off of the count probability at some finite number of particles because of the finite size of a survey. They proposed a way to correct for this bias by extending the tail of the count probability using a few phenomenological parameters calibrated from simulations (see also Fry & Gaztañaga (1994); Munshi et al. (1997)). The method works reasonably well for small smoothing scales, but not for large ones, mainly because the probability-distribution-tail becomes noisy in the latter case (Colombi 1998, private communication). Moreover, the method is feasible only with a dense-sampled survey (Bouchet et al. (1993)). More recently, Szapudi & Colombi (Szapudi & Colombi (1996)) and Colombi et al. Colombi et al. (1998) discussed how a finite survey-volume affects the errors (i.e. the variance around the mean) in the estimators of the N-cumulants $`\overline{\xi }_N`$ (or the related factorial moments), but not the mean or bias of the estimator for $`S_N`$.

However, there has been no attempt to explain quantitatively the differences in the measured $`S_N`$ from different surveys of similar selection criteria in terms of the finite-volume effect. This is, perhaps, in part due to the lack of an analytical estimate of the systematic bias of the standard estimators for $`S_N`$. As we will show, a remarkably simple statistical fact allows just such a calculation to be done, while clarifying the origin of this bias, and relating it to other known biases in large scale structure statistics, such as the integral constraint.

Consider the following elementary statement:

$$\frac{\widehat{A}}{\widehat{B}}\frac{\widehat{A}}{\widehat{B}}$$
(2)   
where $``$ denotes ensemble averaging, and $`\widehat{A}`$ and $`\widehat{B}`$ are two random variables or estimators. This statement holds generically, except in special cases such as when $`\widehat{A}`$ and $`\widehat{B}`$ are constants.

The standard method of estimating $`S_N`$ is to form estimates of $`\overline{\xi }_N`$ and $`\overline{\xi }_2`$ separately, and then take an appropriate ratio of the two. However, even in the ideal case where one has an unbiased estimator of the numerator ($`\overline{\xi }_N`$; let us call this estimator $`\widehat{A}`$, i.e. $`\widehat{A}=\overline{\xi }_N`$) and an unbiased estimator of the denominator ($`\overline{\xi }_2^{N1}`$; let us call the estimator $`\widehat{B}`$, i.e. $`\widehat{B}=\overline{\xi }_2^{N1}`$) <sup>2</sup><sup>2</sup>2We will show shortly that even this ideal case does not hold in reality., taking a ratio of the two estimators does not necessarily result in an unbiased estimator of $`S_N`$. This is captured by the statistical statement in eq. (2). We will refer to this kind of estimation-bias as the ratio-bias.

More generally, nonlinear combinations of unbiased estimators should be treated with great care. For instance, suppose $`\widehat{\overline{\xi }}_2`$ is an unbiased estimator of $`\overline{\xi }_2`$ such that $`\widehat{\overline{\xi }}_2=\overline{\xi }_2`$. It is virtually guaranteed that $`(\widehat{\overline{\xi }}_2)^2(\overline{\xi }_2)^2`$. We will refer to this kind of bias generally as a nonlinear-estimation-bias, of which the ratio-bias is a particularly simple and common example. In this paper, we will sometimes be abusing the terminology by using the two terms interchangeably.

The well-known integral constraint, in the case of measurements of the two-point function, can in fact be understood as a ratio bias. The two-point function $`\xi _2(i,j)`$, where $`i`$ and $`j`$ are two cells separated by some distance, is by definition $`\delta _i\delta _j`$ where $`\delta _i`$ is the overdensity at cell $`i`$. The catch is that one directly observes only $`n_i`$, the number of particles/galaxies in a cell. An estimate of the mean number count $`\overline{n}`$ has to be made in order to assign a value $`\delta _i`$ for each cell. This is generally taken to be $`_in_i/N_T`$ where $`N_T`$ is the total number of cells, let us call this estimator $`\widehat{\overline{n}}`$. The problem is, of course, that $`(n_i\widehat{\overline{n}})(n_j\widehat{\overline{n}})/\widehat{\overline{n}}^2\xi (i,j)`$, because of the estimator $`\widehat{\overline{n}}`$ in the denominator. As we will see, this estimation-bias is in general negative. This is what the integral constraint is about: that the measured two-point function is generally biased low because the mean number density is estimated from the same survey from which $`\xi _2`$ is being measured. This turns out to be the ratio-bias in disguise. It is also easy to see that estimates of $`\xi _N`$ would suffer from a similar bias.

Peebles Peebles (1980) first pointed out the integral-constraint-bias for estimating $`\xi _2`$, but his treatment only gave the large scale estimation-bias. Bernstein Bernstein (1994), building on earlier work by Landy & Szalay Landy & Szalay (1993), developed a perturbative approach (not perturbative in the usual sense of small density fluctuations, but perturbative in the small quantity: the average of the two-point function over the volume of the survey) to compute the full integral-constraint-bias for $`\xi _2`$, which accounted for the small-scale bias as well. (See also Kerscher (1998) for a related recent discussion.) Obviously, the integral-constraint bias affects also measurements of the one-point analogue, or the volume-average, of $`\xi _2`$ i.e. the variance $`\overline{\xi }_2`$. This integral-constraint-bias generally decreases in magnitude with increasing survey size, causing measurements of the correlation length from $`\xi _2`$ or $`\overline{\xi }_2`$ to increase with sample depth, an effect that has been observed before (Davis et al. (1988); Bouchet et al. (1993)). Hence, in the case of the two-point function or its volume average, the finite-volume-effect pointed out by Colombi et al. Colombi et al. (1994) is none other than the integral-constraint-bias.

Adopting the techniques of Bernstein Bernstein (1994), we compute analytically the biases of the standard estimators for $`\overline{\xi }_N`$ and $`S_N`$, for $`N=2,3`$. The methodology for a general $`N`$ is presented in §2. For simplicity, we illustrate how to keep track of the perturbative-ordering by going into details of the calculation for $`N=2,3`$. These cases are also of special interest because many measurements of $`\overline{\xi }_2`$, $`\overline{\xi }_3`$ and $`S_3`$ exist in the literature. We go over the calculation of the estimation-biases for these three quantities in §3. For readers not interested in the details: much of the section can be skipped; the main results are in eq. (27), (28) & (29).

We next check in §4 our analytical results using N-body simulations of the SCDM (Standard Cold-Dark-Matter) and LCDM (Lambda or Low-Density Cold-Dark-Matter) models. Ensemble averages of the standard estimators for $`\overline{\xi }_2`$ and $`S_3`$ are computed by using $`10`$ realizations for each chosen model and for various sample-volumes. The overall agreement is excellent. We also introduce a way to correct for the analytical estimates when the estimation-bias becomes so large that the perturbative approach in §3 breaks down.

In §5, we present a first step towards a re-evaluation of existing measurements of the variance and skewness from the CfA/SSRS/APM surveys. We study simulated CfA/SSRS/APM catalogues with the appropriate sizes, which include the effects of redshift distortions as well as sparse-sampling. We then consider a preliminary correction of the existing measurements of the variance and skewness from these surveys, based on our findings in §3 and §4. The correction is necessarily model-dependent (dependent on the power spectrum and the galaxy-biasing assumed), but it appears that most of the systematic differences between these surveys can be explained by the estimation-biases, under reasonable assumptions. A thorough analysis of the remaining differences would require a careful study of, among other things, projection effects and redshift distortions. We leave this for future work.

Finally, we conclude in §6 with a discussion of methods for measuring $`S_N`$ that might be subject to a less severe estimation-bias. We list other large scale structure statistics which might also suffer from this type of a nonlinear-estimation-bias. We also discuss applications of our findings outside conventional galaxy surveys, such as the Lyman-alpha forest, high-redshift Lyman-break galaxy surveys and weak-lensing maps.

## 2 Biases of the Standard Estimators for $`\overline{\xi }_N`$ and $`S_N`$

### 2.1 Definitions

The standard estimator for $`S_N`$ is given by:

$$\widehat{S}_N=\widehat{\overline{\xi }}_N/(\widehat{\overline{\xi }}_2)^{N1}.$$
(3)   
We use $`\widehat{}`$ to denote estimators of quantities we are interested in.

$`\widehat{\overline{\xi }}_2`$ is the esimator for the variance:

$$\widehat{\overline{\xi }}_2=\frac{1}{N_T}\underset{i}{}(\widehat{\delta }_i)^2\widehat{\overline{\xi }}_2^{\mathrm{shot}}.$$
(4)   
Imagine that the survey is divided into many very small cells, so small that the number of particles/galaxies in each cell is either 1 or 0. The index $`i`$ above denotes such a cell, and $`N_T`$ is the total number of such cells. $`\widehat{\delta }_i`$ is an estimate of the local overdensity smoothed over some given radius $`R`$. We will assume top-hat smoothing in this paper. In other words:

$$\widehat{\delta }_i=\underset{j}{}\frac{n_j\widehat{\overline{n}}}{\widehat{\overline{n}}}W_T(i,j),$$
(5)   
where $`n_j`$ is equal to 1 if there is a galaxy and 0 otherwise, $`\widehat{\overline{n}}`$ is an estimate of the mean density of the survey, and $`W_T(i,j)`$ is the top-hat smoothing window. The estimator $`\widehat{\overline{n}}`$ is

$$\widehat{\overline{n}}=\frac{1}{N_T}\underset{i}{}\underset{j}{}n_jW_T(i,j).$$
(6)   
We have not stated explicitly how to handle edge-effects. For instance, what should be done when the top-hat window $`W_T`$ overlaps with the boundary in eq. (5)? Note that $`_jW_T(i,j)=1`$ only for $`i`$ sufficiently far away from edges. If one adopts the strategy that one picks only $`i`$’s in eq. (5) such that the top-hat does not overlap with the boundary, then the corresponding $`N_T`$ in eq. (4) should not be the total number of all infinitesimal cells in the sample, but a smaller number: the number of centers (at the infinitesimal cells) of top-hats whose windows do not overlap with the boundary. The $`j`$-index in eq. (5) should, on the other hand, range over the whole survey volume, up to the boundary.

The standard shot-noise correction is given by

$$\widehat{\overline{\xi }}_2^{\mathrm{shot}}=\frac{1}{\widehat{N_R}}$$
(7)   
where $`\widehat{N_R}`$ is the estimated mean number of particles within a smoothing top-hat of size $`R`$, in other words, it is $`\widehat{\overline{n}}V_R`$ where $`V_R`$ is the volume of the top-hat.

Phrased in the above manner, the estimator $`\widehat{\overline{\xi }}_2`$ is equivalent to the standard estimator for the variance using the counts-in-cells method, infinitely sampled (see Szapudi (1998)).

The estimator for the N-th cumulant is defined similarly:

$$\widehat{\overline{\xi }}_N=\frac{1}{N_T}\underset{i}{}(\widehat{\delta }_i)_c^N\widehat{\overline{\xi }}_N^{\mathrm{shot}}$$
(8)   
where the subscript $`c`$ denotes the connected part of the sum, and $`\widehat{\overline{\xi }}_N^{\mathrm{shot}}`$ is a shot-noise correction generalizing $`\widehat{\overline{\xi }}_2^{\mathrm{shot}}`$. It is worthwhile at this point to introduce the continuum notation. We will replace $`_i/N_T`$ by $`𝑑V_i/V_T`$, where $`dV_i`$ is the volume of the i-th cell and $`V_T`$ is the total volume, and $`W_T(i,j)`$ by the continuum top-hat $`W(i,j)`$ normalized so that $`𝑑V_jW(i,j)=1`$ for $`i`$ sufficiently far away from edges.

It is also worthwhile to note that factorial moments are sometimes used to estimate $`\overline{\xi }_N`$, which represents a convenient way to eliminate the shot-noise contribution (Szapudi & Szalay (1993)), but otherwise results in the same estimator for $`\overline{\xi }_N`$ as in eq. (8).

### 2.2 Derivation of the Estimation-Biases: an Outline

The integral-constraint-bias for $`\widehat{\overline{\xi }}_2`$ arises from the fact that the true mean density $`\overline{n}`$ is unknown, and has to estimated from the same sample from which one tries to measure the variance. To derive it, let us first write the estimator for the mean density as $`\widehat{\overline{n}}=\overline{n}(1+\alpha )`$, where $`\widehat{\overline{n}}`$ is given in eq. (6), and $`\alpha `$ is a small fluctuation from the true mean. Then one can express the estimator for the overdensity $`\widehat{\delta }_i`$ (eq. ) as

$$\widehat{\delta }_i=(\delta _i\alpha )(1\alpha +\alpha ^2\mathrm{}),$$
(9)   
where $`\delta _i`$ is now the true overdensity (smoothed with the same top-hat as $`\widehat{\delta }_i`$ is) i.e. $`\delta _i`$ is $`n_i/\overline{n}1`$, appropriately smoothed. $`\alpha `$ is by definition equal to

$$\alpha =\frac{1}{V_T}𝑑V_i\delta _i.$$
(10)   
Substituting the above into the expression for $`\widehat{\overline{\xi }}_2`$ (eq. ), one can expand in $`\alpha `$ and write down the ensemble average of $`\widehat{\overline{\xi }}_2`$ order by order. It can be easily shown that the zero-th order term (no $`\alpha `$) of $`\widehat{\overline{\xi }}_2`$ gives $`\overline{\xi }_2`$, the true variance. The rest of the terms represent the integral-constraint-bias. A key question is by what order one should stop, which we will discuss in detail in §3.1. It suffices to note here that, strictly speaking, $`\alpha `$ by itself cannot be used as an ordering-parameter, because it is a random variable (depends on the data) which gets ensemble-averaged. Let us denote the integral-constraint-bias for $`\widehat{\overline{\xi }}_2`$, or more generally, $`\widehat{\overline{\xi }}_N`$ by:

$$\widehat{\overline{\xi }}_N\overline{\xi }_N(1+\frac{\mathrm{\Delta }_{\overline{\xi }_N}}{\overline{\xi }_N}).$$
(11)   
As we will see, the fractional bias $`\mathrm{\Delta }_{\overline{\xi }_N}/\overline{\xi }_N`$ becomes small for a large enough survey. This is the limit in which we will be working. How large is large, or how small is small, is the question we would like to address.

It is easy to see that the above methodology can be adopted for computing the ensemble average of the standard estimator for the hierarchical amplitude, $`\widehat{S}_N`$. The key idea is to assume the denominator part of the estimator fluctuates about its mean, and expand in that fluctuation. In other words, for $`S_N`$, let us first assume

$$\widehat{\overline{\xi }}_2\overline{\xi }_2(1+\frac{\mathrm{\Delta }_{\overline{\xi }_2}}{\overline{\xi }_2})(1+ϵ),$$
(12)   
where $`ϵ`$ is a small fluctuation of the measured $`\widehat{\overline{\xi }}_2`$ from its mean (which is offset from the true $`\overline{\xi }_2`$ by a bias; eq. ). Note that $`ϵ`$ depends on the data, and so cannot be taken out of $``$.

Putting eq. (11) and (12) into eq. (3), one obtains the following expression for the mean of $`\widehat{S}_N`$:

$`\widehat{S}_N`$ $`=`$ $`S_N(1+{\displaystyle \frac{\mathrm{\Delta }_{S_N}}{S_N}})`$ (13)   
$`{\displaystyle \frac{\mathrm{\Delta }_{S_N}}{S_N}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }_{\overline{\xi }_N}}{\overline{\xi }_N}}(N1){\displaystyle \frac{\mathrm{\Delta }_{\overline{\xi }_2}}{\overline{\xi }_2}}(N1){\displaystyle \frac{ϵ\widehat{\overline{\xi }}_N}{\overline{\xi }_N}}+(N1){\displaystyle \frac{N}{2}}{\displaystyle \frac{ϵ^2\widehat{\overline{\xi }}_N}{\overline{\xi }_N}}+\mathrm{}.`$ (14)   
Note that all terms up to $`ϵ^2`$ are displayed above, except for the terms: $`(N1)(\mathrm{\Delta }_{\overline{\xi }_2}/\overline{\xi }_2)[ϵ\widehat{\overline{\xi }}_N/\overline{\xi }_N]`$ and $`(N1)\frac{N}{2}(\mathrm{\Delta }_{\overline{\xi }_2}/\overline{\xi }_2)[ϵ^2\widehat{\overline{\xi }}_N/\overline{\xi }_N]`$. As we will show later on, the terms $`\mathrm{\Delta }_{\overline{\xi }_2}/\overline{\xi }_2`$, $`ϵ\widehat{\overline{\xi }}_N/\overline{\xi }_N`$ and $`ϵ^2\widehat{\overline{\xi }}_N/\overline{\xi }_N`$ are all of the order of a small parameter in which we will be expanding: hence the dropping of products involving them.

It is instructive to divide the net bias of the estimator $`\widehat{S}_N`$ into two different contributions. One is due to the integral-constraint-bias discussed above: namely that both $`\widehat{\overline{\xi }}_N`$ in the numerator and $`\widehat{\overline{\xi }}_2`$ in the denominator (see eq. ) are biased estimators. This gives the first two terms on the right hand side of eq. (14). The other contribution arises from the fact that the estimator $`\widehat{S}_N`$ is the ratio of two other estimators (to be accurate, it is in fact some nonlinear combinations of $`\widehat{\overline{\xi }}_N`$ and $`\widehat{\overline{\xi }}_2`$). This is the rest of the terms in eq. (14). Let us give these two kinds of terms explicit names:

$`{\displaystyle \frac{\mathrm{\Delta }_{S_N}^{\mathrm{int}.\mathrm{constr}.}}{S_N}}={\displaystyle \frac{\mathrm{\Delta }_{\overline{\xi }_N}}{\overline{\xi }_N}}(N1){\displaystyle \frac{\mathrm{\Delta }_{\overline{\xi }_2}}{\overline{\xi }_2}},`$ (15)   
$`{\displaystyle \frac{\mathrm{\Delta }_{S_N}^{\mathrm{ratio}}}{S_N}}=(N1){\displaystyle \frac{ϵ\widehat{\overline{\xi }}_N}{\overline{\xi }_N}}+(N1){\displaystyle \frac{N}{2}}{\displaystyle \frac{ϵ^2\widehat{\overline{\xi }}_N}{\overline{\xi }_N}},`$ (16)   
the integral-constraint-bias and the ratio-bias of $`\widehat{S}_N`$ respectively. We are slightly abusing the terminology here because the integral-constraint-bias is of course itself a form of a ratio-bias. As we will show below, it turns out that the terms contributing to the integral-constraint-bias partially cancel each other, leaving the ratio-bias to be the dominant contribution to the net bias of $`\widehat{S}_N`$ on large scales.

Note that we have implicitly assumed $`\mathrm{\Delta }_{\overline{\xi }_N}/\overline{\xi }_N`$, $`(N1)\mathrm{\Delta }_{\overline{\xi }_2}/\overline{\xi }_2`$ and $`(N1)\widehat{\overline{\xi }}_Nϵ/\overline{\xi }_N`$ and $`(N1)Nϵ^2\widehat{\overline{\xi }}_N/2\overline{\xi }_N`$ are all small, and that terms we have ignored in eq. (14) are somehow all of higher order. What is the correct ordering parameter in which we are expanding? Again, the quantity $`ϵ`$ cannot be used directly to keep track of the ordering, because it is data-dependent. For instance, after taking the expectation values, it could happen that terms that contain $`ϵ^2`$ are actually comparable to terms linear in $`ϵ`$. We will see that this is indeed the case in the next section.

eq. (13) implies the bias of $`\widehat{S}_N`$ depends in general on the $`M`$-point correlation functions. The strategy we adopt in this paper is to assume the following hierarchical relation: $`\xi _M\xi _2^{M1}`$, which is motivated by perturbation theory but has been observed to hold in the highly nonlinear regime as well. We do not need to assume anything, however, about the configuration or scale dependence (or independence) of the hierarchical amplitudes. Using this relation, it can be shown that (we will demonstrate this explicitly for $`S_3`$) the terms we have kept in $`\mathrm{\Delta }_{S_N}/S_N`$ (eq. ) all contain terms linear in the following quantity:

$$\overline{\xi }_2^L\frac{1}{V_T^2}𝑑V_i𝑑V_j\xi _2(i,j)$$
(17)   
where $`\xi _2`$ is a smoothed version of the 2-point function defined as

$$\xi _2(i,j)=𝑑V_k𝑑V_l\xi _2^{\mathrm{usmth}.}(k,l)W(k,i)W(l,j)$$
(18)   
where $`W`$ is the top-hat smoothing window of some radius $`R`$, and $`\xi _2^{\mathrm{usmth}.}`$ is the unsmoothed 2-point function. $`\overline{\xi }_2^L`$ is the 2-point function averaged over the whole survey (of size $`L`$). For a survey to be of any use at all, this quantity has to be much smaller than 1. Because of the relatively large coefficients that will be multiplying it in $`\mathrm{\Delta }_{S_N}/S_N`$, as will be shown below, the fractional bias in $`\widehat{S}_N`$ could be non-negligible, even for a relatively large survey volume.

One should therefore view the derivation of the bias in $`\widehat{S}_N`$ in this paper as an expansion in the small parameter $`\overline{\xi }_2^L`$. As we will see, the same applies to $`\widehat{\overline{\xi }}_N`$. Indeed, it can be shown that all the terms we have kept in eq. (13) contain terms linear in $`\overline{\xi }_2^L`$, and terms we have ignored are all of higher order ($`[\overline{\xi }_2^L]^2`$, etc). We will demonstrate the reasoning with an example: $`N=3`$.

## 3 Estimation-Biases for the Variance, the Third Moment and the Skewness

### 3.1 The Integral-Constraint-Bias for $`\widehat{\overline{\xi }}_2`$

The integral-constraint-bias for the two-point function has been known for a long time (see e.g. Peebles (1980); Bernstein (1994); Tegmark et al. (1998)). Our treatment follows most closely that of Bernstein Bernstein (1994), and the emphasis is on techniques that can be generalized to $`\widehat{\overline{\xi }}_3`$ and $`\widehat{S}_3`$.

Following the strategy outlined in §2.2, namely combining eq. (9), (10) with (4), we obtain

$`\widehat{\overline{\xi }}_2=`$ $`{\displaystyle \frac{1}{V_T}}{\displaystyle 𝑑V_i\delta _i^2}\widehat{\overline{\xi }}_2^{\mathrm{shot}}`$   
$`{\displaystyle \frac{1}{V_T^2}}{\displaystyle 𝑑V_i𝑑V_j\delta _i\delta _j}{\displaystyle \frac{2}{V_T^2}}{\displaystyle 𝑑V_i𝑑V_j\delta _i^2\delta _j}`$   
$`+{\displaystyle \frac{3}{V_T^3}}{\displaystyle 𝑑V_i𝑑V_j𝑑V_k\delta _i^2\delta _j\delta _k}+{\displaystyle \frac{4}{V_T^3}}{\displaystyle 𝑑V_i𝑑V_j𝑑V_k\delta _i\delta _j\delta _k}+\mathrm{}`$   
where we display all terms up to $`O(\alpha ^2)`$ (see eq. ).

As we have explained before, $`\alpha `$ cannot really be used as an ordering-parameter, because it is a random variable which gets averaged over. The key to the above expansion is instead to keep only terms up to linear order in $`\overline{\xi }_2^L`$ (eq. ).

Ignoring shot-noise for now, the first term on the right gives us the true variance $`\overline{\xi }_2`$. This is the zero-th order term. The rest of the terms represent the integral-constraint-bias for $`\widehat{\overline{\xi }}_2`$, and are all of order $`\overline{\xi }_2^L`$ or higher. Let us check.

1. The first term on the second line of eq. (3.1) is none other than $`\overline{\xi }_2^L`$ itself (eq. ).

2. The next term, again ignoring shot-noise for now, gives $`(2/V_T^2)𝑑V_i𝑑V_j\xi _3(i,i,j)`$. Applying the hierarchical relation $`\xi _M\xi _2^{M1}`$, i.e. $`\xi _3(i,i,j)\xi _2(i,i)\xi _2(i,j)+\mathrm{}`$, we can see that the term $`\xi _2(i,i)\xi _2(i,j)=\overline{\xi }_2\xi _2(i,j)`$, when integrated over $`i`$ and $`j`$, will give rise to a term proportional to $`\xi _2^L`$. Hence, the term $`(2/V_T^2)𝑑V_i𝑑V_j\xi _3(i,i,j)`$ contains a linear piece and should not be thrown away. Note that we need not make any assumption about the configuration or scale dependence of the hierarchical amplitudes. In other words, we should keep the term $`(2/V_T^2)𝑑V_i𝑑V_j\xi _3(i,i,j)`$ as is, rather than as, say $`(2/V_T^2)\overline{\xi }_2𝑑V_i𝑑V_j\xi _2(i,j)`$. The hierarchical relation is used strictly for keeping track of the ordering at this stage.

3. The integrand in the next term, 3 $`\delta _i^2\delta _j\delta _k`$, can be broken up into several disconnected pieces: $`3\delta _i^2\delta _j\delta _k=3[\xi _2(i,i)\xi _2(j,k)+2\xi _2(i,j)\xi _2(i,k)+\xi _4(i,i,j,k)]`$. It is easy to see that the second piece gives something that is second order in $`\overline{\xi }_2^L`$, when integrated over, and so does the $`\xi _4(i,i,j,k)`$ piece, assuming again the hierarchical relation. The only piece that survives is then, after integration, $`3\overline{\xi }_2\overline{\xi }_2^L`$.

4. Lastly, the $`(4/V_T^3)𝑑V_i𝑑V_j𝑑V_k\delta _i\delta _j\delta _k`$ term in eq. (3.1) is second order in $`\overline{\xi }_2^L`$, again by applying the hierarchical relation.

More generally, it can be seen that all terms, including those not explicitly displayed, in the expansion in eq. (3.1) are of the form $`V_T^m𝑑V_{i_1}\mathrm{}𝑑V_{i_m}\delta _{i_1}^\gamma \delta _{i_2}\mathrm{}\delta _{i_m}`$, where $`\gamma `$ is 1 or 2. Our arguments above show that only terms with $`m=1`$, or $`m=2`$, or $`m=3`$ and $`\gamma =2`$, contain pieces linear in $`\overline{\xi }_2^L`$. The reader can convince himself or herself that all other terms are of higher order.

Putting everything together, ignoring shot-noise, and using the definition of the fractional bias in eq. (11), we obtain to linear order in $`\overline{\xi }_2^L`$,

$`{\displaystyle \frac{\mathrm{\Delta }_{\overline{\xi }_2}}{\overline{\xi }_2}}=`$ $`{\displaystyle \frac{1}{\overline{\xi }_2V_T^2}}{\displaystyle 𝑑V_i𝑑V_j\xi _2(i,j)}{\displaystyle \frac{2}{\overline{\xi }_2V_T^2}}{\displaystyle 𝑑V_i𝑑V_j\xi _3(i,i,j)}`$   
$`+{\displaystyle \frac{3}{V_T^2}}{\displaystyle 𝑑V_i𝑑V_j\xi _2(i,j)}.`$   
The above result is consistent with that of Bernstein Bernstein (1994) for the integral-constraint-bias of the two-point function. The first term on the right was obtained by Peebles Peebles (1980) via a different argument.

How about shot-noise? A term like $`\delta _i\delta _j`$ that shows up in eq. (3.1) includes both the cosmic 2-point correlation $`\xi _2(i,j)`$ (see eq. ) and a shot-noise contribution due to Poisson sampling (see e.g. Feldman et al. (1994)):

$$\delta _i\delta _j=\xi _2(i,j)+\frac{1}{\overline{n}}𝑑V_kW(k,i)W(k,j)$$
(21)   
It can be shown that all shot-noise contributions to $`\mathrm{\Delta }_{\overline{\xi }_2}/\overline{\xi }_2`$ can be expanded in either $`1/N_R`$ or $`1/N_L`$ where $`N_R`$ is the mean number of particles in a cell of size $`R`$, and $`N_L`$ is the mean number of particles in the whole survey. For instance, the term $`V_T^1𝑑V_i\delta _i^2`$ has a Poisson term $`1/N_R`$, if one makes use of eq. (21), together with the fact that $`W`$ is a top-hat of size $`R`$ with the normalization $`𝑑V_jW(i,j)=1`$. This gets canceled by part of the shot-noise correction term $`\widehat{\overline{\xi }}_2^{\mathrm{shot}}`$ (eq. \[3.1\]). On the other hand, a term like $`V_T^2𝑑V_i𝑑V_j\delta _i\delta _j`$ gives us a Poisson term of the order of $`1/N_L`$. In general, $`1/N_L`$ is a very small quantity, and so we can ignore all terms of order $`1/N_L`$. There are, for instance, Poisson pieces linear in $`\overline{\xi }_2^L`$ in the $`\delta _i\delta _j\delta _k`$ term in eq. (3.1), but they are also of order $`1/N_L`$, and so can be ignored.

How about the $`O(1/N_R)`$ terms in eq. (3.1)? Without going into details, it can be shown that the only $`O(1/N_R)`$ contributions are: $`\overline{\xi }_2^L/N_R`$ from the $`\widehat{\overline{\xi }}_2^{\mathrm{shot}}`$ term, $`2\overline{\xi }_2^L/N_R`$ from the fourth term on the right, and $`3\overline{\xi }_2^L/N_R`$ from the fifth term on the right. Therefore, to $`O(1/N_R)`$, the shot-noise terms miraculously cancel! It is interesting to note that this cancellation is possible only because the integral-constraint-bias in the shot-noise correction itself is taken into account i.e. $`\widehat{\overline{\xi }}_2^{\mathrm{shot}}=1/\widehat{N}_R=(1+\alpha ^2\mathrm{})/N_R`$ (see eq.  & ).

We will see in §4 that our estimations of the biases in the variance and skewness which ignores shot-noise are in fact quite accurate. With no further justification, we will ignore shot-noise terms in the rest of our derivation, which substantially simplifies our expressions. The cancellation to $`O(1/N_R)`$ here for $`\mathrm{\Delta }_{\overline{\xi }_2}/\overline{\xi }_2`$ can be taken as suggestive evidence that shot-noise is unimportant for the estimation-biases we are interested in; we will rely on the numerical simulations in §4 for further support. We should emphasize, however, we are not saying that there is no need to subtract out shot-noise when estimating the variance and skewness themselves.

Finally, how about edge-effects? It is worth emphasizing that no assumptions about the edge-effects being small need to be made in deriving eq. (3.1). One only has to be careful about volume over which the integration is done and what $`V_T`$ has to be. As we have mentioned earlier in §2.1, one way to deal with the boundary is to use only cells that do not overlap with the edges. In that case, the dummies of integration $`i`$ and $`j`$ should range over the inner part of the survey, where any cell that is centered within it does not cut into the boundary. Similarly, $`V_T`$ should be chosen to be the volume of that inner region. We will discuss how to approximate such integrals in §3.4.

### 3.2 The Integral-Constraint-Bias for $`\widehat{\overline{\xi }}_3`$

Similarly, one can derive the integral-constraint-bias for the third cumulant using eq. (8), (11) and (9) . Ignoring shot-noise, and again, keeping only terms to first order in $`\overline{\xi }_2^L`$, the integral-constraint-bias for $`\overline{\xi }_3`$ is given by

$`{\displaystyle \frac{\mathrm{\Delta }_{\overline{\xi }_3}}{\overline{\xi }_3}}=`$ $`{\displaystyle \frac{3}{\overline{\xi }_3V_T^2}}{\displaystyle 𝑑V_i𝑑V_j\xi _3(i,i,j)}{\displaystyle \frac{9\overline{\xi }_2}{\overline{\xi }_3V_T^2}}{\displaystyle 𝑑V_i𝑑V_j\xi _2(i,j)}`$   
$`{\displaystyle \frac{3}{\overline{\xi }_3V_T^2}}{\displaystyle 𝑑V_i𝑑V_j\xi _4(i,i,i,j)}+{\displaystyle \frac{9\overline{\xi }_2}{\overline{\xi }_3V_T^2}}{\displaystyle 𝑑V_j𝑑V_k\xi _2(j,k)}`$   
$`+{\displaystyle \frac{6}{V_T^2}}{\displaystyle 𝑑V_j𝑑V_k\xi _2(j,k)}`$   
Note that by essentially the same reasoning as in the case of $`\widehat{\overline{\xi }}_2`$, one only needs to consider terms up to $`\alpha ^2`$ in deriving the above.

### 3.3 The Estimation-Bias $`\widehat{S}_3`$

The integral-constraint-bias for $`\widehat{S}_3`$ can be simply read off from eq. (15), (3.1) & (3.2). The ratio-bias for $`\widehat{S}_3`$, on the other hand, follows from eq. (16). Substituting the definition of $`ϵ`$ from eq. (12), we obtain

$`{\displaystyle \frac{\mathrm{\Delta }_{S_3}^{\mathrm{ratio}}}{S_3}}=`$ $`{\displaystyle \frac{2}{\overline{\xi }_3\overline{\xi }_2(1+\frac{\mathrm{\Delta }_{\overline{\xi }_2}}{\overline{\xi }_2})}}\widehat{\overline{\xi }}_3[\widehat{\overline{\xi }}_2\overline{\xi }_2(1+{\displaystyle \frac{\mathrm{\Delta }_{\overline{\xi }_2}}{\overline{\xi }_2}})]`$   
$`+{\displaystyle \frac{3}{\overline{\xi }_3\overline{\xi }_2^2(1+\frac{\mathrm{\Delta }_{\overline{\xi }_2}}{\overline{\xi }_2})^2}}\widehat{\overline{\xi }}_3[\widehat{\overline{\xi }}_2\overline{\xi }_2(1+{\displaystyle \frac{\mathrm{\Delta }_{\overline{\xi }_2}}{\overline{\xi }_2}})]^2.`$   
Substituting eq. (3.1), (3.2) and (3.3) into eq. (13) for $`N=3`$, ignoring shot-noise terms and keeping only terms linear in $`\overline{\xi }_2^L`$, it can be shown that

$`{\displaystyle \frac{\mathrm{\Delta }_{S_3}^{\mathrm{ratio}}}{S_3}}=`$ $`{\displaystyle \frac{6}{V_T^2\overline{\xi }_3}}{\displaystyle 𝑑V_i𝑑V_j\xi _3(i,j,j)}{\displaystyle \frac{2}{V_T^2\overline{\xi }_2\overline{\xi }_3}}{\displaystyle 𝑑V_i𝑑V_j\xi _5(i,i,i,j,j)}`$   
$`+{\displaystyle \frac{3}{V_T^2\overline{\xi }_2^2}}{\displaystyle 𝑑V_j𝑑V_k\xi _4(j,j,k,k)},`$   
where $`\xi _N`$ is the connected N-point function. The first two terms here arise from the first term on the right hand side of eq. (3.3), and the last term from the second one.

The reasoning used to arrive at the above expression is again very similar to the case of $`\widehat{\overline{\xi }}_2`$. But there are a few new tips to keep in mind.

1. Consider a term like $`\delta _i^3[\delta _j^2\overline{\xi }_2(1+\mathrm{\Delta }_{\overline{\xi }_2}/\overline{\xi }_2)]`$, which arises from the first term on the right of eq. (3.3). It can be seen that a disconnected piece with the $`j`$ index all by itself, such as $`\delta _i^3\delta _j^2`$, is going to be canceled, by $`\delta _i^3\overline{\xi }_2(1+\mathrm{\Delta }_{\overline{\xi }_2}/\overline{\xi }_2)`$. More generally, it can be seen that all terms in the expansion of eq. (3.3) contain an integrand of the form $`\delta _i^3\delta _{j_1}^2\delta _{j_2}^2\mathrm{}\delta _{j_m}^2`$, any disconnected piece of which with the $`j_1`$, or … $`j_m`$ index all by itself is going to get canceled. In other words, the $`j`$-indices must be connected to each other, or to $`i`$.

2. As before, integrals over products of the two-point function, of the form $`\delta _i\delta _j\delta _j\delta _k`$ for instance, are of higher order. An example is from the second term on the right of eq. (3.3). It contains a term with $`\delta _i^3\delta _j^2\delta _k^2`$ in the integrand. This can be broken up into several disconnected pieces. Most of them got canceled because of the reason laid out in 1. above. One disconnected piece that does not get canceled is $`\delta _i^3\delta _j\delta _j\delta _k^2`$. However, applying the hierarchical relation as before, we can see that this piece, when integrated over, is second order in $`\overline{\xi }_2^L`$. Another disconnected piece is: $`\delta _i^3\delta _j\delta _k^2`$. This deserves special attention. It gives rise to a term in the estimation-bias of the order of $`V_T^2𝑑V_i𝑑V_j[\xi _2(i,j)]^2`$, which is not guaranteed to be much smaller than $`\overline{\xi }_2^L`$ in general. However, we have checked numerically that for realistic power spectra, and $`L`$ larger than about $`10\mathrm{h}^1\mathrm{Mpc}`$, they are indeed small compared to $`\overline{\xi }_2^L`$. To summarize, we can say that any products of the two-point function where the arguments are non-degenerate (i.e. $`\xi _2[i,j]`$ where $`ij`$) can be ignored.

3. Combining 1. and 2., it can be seen that all $`ϵ^3`$ terms (or higher) can be ignored from eq. (3.3). The same argument applies for arbitrary $`N`$, and so justifies the dropping of $`ϵ^3`$ terms from eq. (14) or (16).

4. Note that eq. (3.3) can be derived by ignoring $`\mathrm{\Delta }_{\overline{\xi }_2}`$ and $`\mathrm{\Delta }_{\overline{\xi }_3}`$ altogether. Since, the lowest order terms in eq. (3.3) are already of order $`\overline{\xi }_2^L`$, including $`\mathrm{\Delta }_{\overline{\xi }_2}`$ and $`\mathrm{\Delta }_{\overline{\xi }_3}`$ can only give higher order terms.

Lastly, combining eq. (3.1), (3.2) and (3.3) and substituting into eq. (14) and (16), we obtain the net bias of the estimator $`\widehat{S}_3`$:

$`{\displaystyle \frac{\mathrm{\Delta }_{S_3}}{S_3}}=`$ $`2{\displaystyle \frac{\overline{\xi }_2^L}{\overline{\xi }_2}}+{\displaystyle \frac{4}{\overline{\xi }_2V_T^2}}{\displaystyle 𝑑V_i𝑑V_j\xi _3(i,j,j)}`$   
$`{\displaystyle \frac{9}{\overline{\xi }_3V_T^2}}{\displaystyle 𝑑V_i𝑑V_j\xi _3(i,i,j)}{\displaystyle \frac{3}{\overline{\xi }_3V_T^2}}{\displaystyle 𝑑V_i𝑑V_j\xi _4(i,i,i,j)}`$   
$`{\displaystyle \frac{2}{\overline{\xi }_2\overline{\xi }_3V_T^2}}{\displaystyle 𝑑V_i𝑑V_j\xi _5(i,i,i,j,j)}+{\displaystyle \frac{3}{\overline{\xi }_2^2V_T^2}}{\displaystyle 𝑑V_j𝑑V_k\xi _4(j,j,k,k)}`$   
### 3.4 An Analytical Approximation

The expressions in eq. (3.1), (3.2) & (3.3) give the exact fractional bias in $`\widehat{\overline{\xi }}_2`$, $`\widehat{\overline{\xi }}_3`$ and $`S_3`$ to first order in $`\overline{\xi }_2^L`$, excluding shot-noise. No assumption about the two-point function $`\xi _2`$ itself being small has been made. The hierarchical relation $`\xi _N\xi _2^{N1}`$ has only been used for book-keeping. We have not assumed anything about the configuration or scale dependence of the hierarchical amplitudes.

In the present form, these expressions are not very useful as the N-point functions, up to $`N=5`$, are required to compute the estimation-biases. We will approximate them as products of the two-point function (it is from this point on, that we use the hierarchical relation for more than simply book-keeping) using the following relation (see Bernardeau (1994)):

$$\delta _i^m\delta _j^m^{}_c=c_{mm^{}}\overline{\xi }_2^{m+m^{}2}\xi _2(i,j),$$
(26)   
where $`\delta _i^m\delta _j^m^{}_c`$ is the connected cosmic $`m+m^{}`$-point function (no Poisson terms), with only at most two differing indices.

Putting eq. (26) into eq. (3.1), (3.2) and (3.3) respectively, we obtain:

$$\frac{\mathrm{\Delta }_{\overline{\xi }_2}}{\overline{\xi }_2}=\frac{\overline{\xi }_2^L}{\overline{\xi }_2}+[32c_{12}]\overline{\xi }_2^L,$$
(27)   
$$\frac{\mathrm{\Delta }_{\overline{\xi }_3}}{\overline{\xi }_3}=3\frac{c_{12}}{S_3}\frac{\overline{\xi }_2^L}{\overline{\xi }_2}+\left[3\frac{c_{13}}{S_3}+6\right]\overline{\xi }_2^L,$$
(28)   
$$\frac{\mathrm{\Delta }_{S_3}}{S_3}=\left[29\frac{c_{12}}{S_3}\right]\frac{\overline{\xi }_2^L}{\overline{\xi }_2}+\left[4c_{12}3\frac{c_{13}}{S_3}2\frac{c_{23}}{S_3}+3c_{22}\right]\overline{\xi }_2^L.$$
(29)   
It is also instructive to distinguish between the two different contributions to $`\mathrm{\Delta }_{S_3}/S_3`$, as in eq. (15) & (16), one from the integral-constraint-biases of $`\widehat{\overline{\xi }}_2`$ and $`\widehat{\overline{\xi }}_3`$ themselves:

$$\frac{\mathrm{\Delta }_{S_3}^{\mathrm{int}.\mathrm{cr}.}}{S_3}=\left[23\frac{c_{12}}{S_3}\right]\frac{\overline{\xi }_2^L}{\overline{\xi }_2}+\left[4c_{12}3\frac{c_{13}}{S_3}\right]\overline{\xi }_2^L$$
(30)   
and the other from the ratio-bias due to the division of $`\widehat{\overline{\xi }}_3`$ by $`\widehat{\overline{\xi }}_2^2`$:

$$\frac{\mathrm{\Delta }_{S_3}^{\mathrm{ratio}}}{S_3}=6\frac{c_{12}}{S_3}\frac{\overline{\xi }_2^L}{\overline{\xi }_2}+\left[3c_{22}2\frac{c_{23}}{S_3}\right]\overline{\xi }_2^L$$
(31)   
In essence then, there are basically two kinds of terms in the estimation-biases, one that does not change with the smoothing scale $`R`$ (the $`\overline{\xi }_2^L`$ term), and the other that increases in magnitude as $`R`$ approaches the size of the survey $`L`$ (the $`\overline{\xi }_2^L/\overline{\xi }_2`$ term). We can write this in general as:

$$\frac{\mathrm{\Delta }_E}{E}=\alpha _1\frac{\overline{\xi }_2^L}{\overline{\xi }_2}+\alpha _2\overline{\xi }_2^L$$
(32)   
where $`\mathrm{\Delta }_E/E`$ denotes the fractional estimation-bias for the estimator $`\widehat{E}`$, and $`\alpha _1`$ and $`\alpha _2`$ are coefficients that depend on various hierarchical amplitudes, such as $`S_3`$ and $`c_{mm^{}}`$.

The relation in eq. (26) is motivated by perturbation theory, and so, strictly speaking, only holds in the weakly nonlinear regime. But we have reasons to believe (from N-body work in preparation) that the same form should work on non-linear scales, albeit with the coefficients $`c_{mm^{}}`$ slightly altered from the perturbative (tree order) values (as it is known to happen with $`S_N`$). In the same vein, we will use the tree order value for $`S_3`$ in the above estimates of the fractional biases. This is admittedly crude for small scales, but we will see this is not a bad approximation in the next section. For $`\overline{\xi }_2`$ and $`\overline{\xi }_2^L`$ on the other hand, we have used both the linear and non-linear values and found little difference in the predicted biases. We use the linear values in all the figures of this paper.

The perturbative values for the various hierarchical amplitudes are (Bernardeau (1994); ignoring galaxy-bias):

$`S_3`$ $`=`$ $`34/7+\gamma `$   
$`c_{12}`$ $`=`$ $`68/21+\gamma /3`$   
$`c_{13}`$ $`=`$ $`11710/441+61\gamma /7+2\gamma ^2/3`$   
$`c_{22}`$ $`=`$ $`c_{12}^2`$   
$`c_{23}`$ $`=`$ $`c_{12}c_{13}`$ (33)   
where $`\gamma =\gamma (R)`$:

$$\gamma \frac{d\mathrm{log}\overline{\xi }_2}{d\mathrm{log}R}$$
(34)   
is the logarithmic slope of the variance. For a power-law power spectrum with spectral index $`n`$ we have: $`\gamma =(n+3)`$.

Substituting the above into eq. (27) to (31), it can be seen that, for $`n`$ of interests, a) the overall estimation-biases in $`\widehat{\xi }_2`$, $`\widehat{\xi }_3`$ and $`\widehat{S}_3`$ are negative, i.e. $`\alpha _i<0`$ in eq. (32); b) the integral-constraint-bias and ratio-bias contributions to $`\mathrm{\Delta }_{S_3}/S_3`$ are comparable on small smoothing scales and c) the ratio-bias contribution to $`\mathrm{\Delta }_{S_3}/S_3`$ dominates on large scales. All these are illustrated in Fig. 1 which shows the coefficients $`\alpha _1`$ (continuous line) and $`\alpha _2`$ (dashed line) as a function of $`\gamma `$, for different estimators $`\widehat{E}=\widehat{\overline{\xi }_2},\widehat{S}_3`$.

Note that for a Gaussian model where $`c_{ij}=S_J=0`$, the estimation-biases are quite different. The coefficient $`\alpha _2`$ becomes $`3`$ for $`\widehat{E}=\widehat{\overline{\xi }_2}`$ so that the bias is positive on small scales. Also, $`\alpha _1=0`$ and $`\alpha _2=6`$ for $`\widehat{E}=\widehat{\overline{\xi }_3}`$ while $`\alpha _1=2`$ and $`\alpha _2=0`$ for $`\widehat{E}=\widehat{S}_3`$, which means the bias is always positive. Hence the Gaussian prediction for the estimation-bias can be quite misleading (even for models with Gaussian initial conditions).

Finally, a word on what value to assume for $`\overline{\xi }_2^L`$. As we have noted before in §3.1, care should be taken in dealing with the edge-effects. The expression for $`\overline{\xi }_2^L`$ is given in eq. (17) and (18). In obtaining $`\widehat{\overline{\xi }}_2`$ and $`\widehat{\overline{\xi }}_3`$, if one insists on using only cells that do not overlap with the boundary, then one has to restrict $`i`$ and $`j`$ in eq. (17) over an inner region of the survey where any cell centered within it does not cut the edges, and one should equate $`V_T`$ with the volume of this inner region. On the other hand, the $`k`$ and $`l`$ indices of eq. (18) should still range over the whole survey volume. We will not try to compute $`\overline{\xi }_2^L`$ exactly. Instead, we make use of the following observation: adhering to the convention for $`i`$, $`j`$, $`k`$ and $`l`$ above, we can rewrite $`\overline{\xi }_2^L`$ as $`𝑑V_k𝑑V_lf_kf_l\xi _2^{\mathrm{usmth}.}(k,l)`$, where $`f_kV_T^1𝑑V_iW(k,i)`$. The quantity $`f_k`$ varies between $`1/V_T`$, for $`k`$ sufficiently far away from the edges, to $`0`$, for $`k`$ sitting on the boundary. If one very crudely replaces $`f_k`$ by its volume-average, which is equal to $`1/V_T^{}`$ where $`V_T^{}`$ is the total volume of the survey (everything within the boundary), it can be seen that $`\overline{\xi }_2^L`$ is then simply equal to $`\overline{\xi }_2`$, except with a funny top-hat that covers the whole volume of the survey. We will further approximate this by estimating $`\overline{\xi }_2^L`$ using $`\overline{\xi }_2`$ with a spherical top-hat of size $`R_L`$ such that its volume is the same as that of the survey (e.g. if the survey is a cubical box of side-length $`L`$, then $`L^3=4\pi R_L^3/3`$). This is admittedly crude, but seems to be sufficiently accurate for the N-body experiments we study, at least for a smoothing scale $`R`$ which is not too large compared to the size of the survey. In practice, one might want to go back to the original definition of $`\overline{\xi }_2^L`$ (eq. ), and compute $`\overline{\xi }_2^L`$ more carefully.

## 4 Comparison with N-body Simulations

To test our analytical predictions in the last section, we use simulations of two different spatially flat cold dark matter (CDM) dominated models. One set of simulations is of the SCDM model, with $`\mathrm{\Omega }_0=1,h=0.5`$, and another is of the LCDM model with $`\mathrm{\Omega }_0=0.2,h=1`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.8`$. The power spectra, P(k) for these models are taken from Bond & Efstathiou Bond & Efstathiou (1984) and Efstathiou, Bond and White Efstathiou et al. (1992). The shape of P(k) is parametrized by the quantity $`\mathrm{\Gamma }=\mathrm{\Omega }h`$, so that we have $`\mathrm{\Gamma }=0.5`$ CDM and $`\mathrm{\Gamma }=0.2`$ CDM. Each simulation contains $`10^6`$ particles in a box of comoving side-length $`300h^1\mathrm{Mpc}`$ and was run using a P<sup>3</sup>M $`N`$-body code (Hockney & Eastwood (1981); Efstathiou et al. (1985)). All outputs are normalized to $`\sigma _8=1`$. The simulations are described in more detail in Dalton et al. Dalton et al. (1994) and Baugh et al.Baugh et al. (1995). In this paper we use 10 realizations of each model for computing ensemble averages, with error bars being estimated from their standard deviation.

¿From each realization we extract one subsample within a cubical box of size $`L=300h^1\mathrm{Mpc}/M`$, where $`M`$ is taken to be an integer, $`M=1,2,3,\mathrm{},7`$. So we have a set of 10 realizations of subsamples for each box-size, from $`L40h^1\mathrm{Mpc}`$ to $`300h^1\mathrm{Mpc}`$. We estimate moments of counts-in-cells (as in Baugh et al. (1995)) in each set of subsamples to study how the estimation biases of the variance and skewness vary with survey volume.

Note the importance of using subsamples of large simulations, rather than running simulations with an intrinsically small box-size. The latter introduces dynamical effects due to the missing of the large scale power, which we are not interested in for the purpose of this paper. Note also the importance of taking one subsample from each realization, rather than extracting multiple subsamples from a single realization, to ensure the statistical independence of the subsamples.

### 4.1 The Variance $`\overline{\xi }_2`$

The results for the variance of the LCDM model are shown in Figure 2. Because the LCDM has more power on large scales this model shows a more pronounced integral-constraint-bias compared to the SCDM model. Open circles show the measured variance averaged over the 10 realizations of the full box ($`L=300h^1\mathrm{Mpc}`$). Filled triangles show the mean measured variance for the smaller boxes with the box-size $`L=300/2,300/4,300/5,300/7h^1\mathrm{Mpc}`$ as indicated in each panel. In all cases, the error-bars represent $`1\sigma `$ deviations in the measured variance over the relevant 10 realizations. The solid line gives the linear perturbation theory prediction for $`\overline{\xi }_2`$. (This is not perturbation theory in the sense of §2.2, but perturbation theory in the usual sense: an expansion in $`\overline{\xi }_2`$ or the density fluctuation amplitude; to avoid confusion, we will refer to it simply as PT.) The agreement of the solid line and the open circles on large scales indicate that the measured variance from the full box does not suffered from a significant bias, for the smoothing scales shown. The dashed line is the integral-constraint-bias prediction in eq., which is in excellent agreement with the simulation results.

### 4.2 The Skewness $`S_3`$

The results for the skewness are shown in Fig. 3 & 4 for the SCDM and LCDM models. Again, because the LCDM has more power on large scales, this model shows a larger estimation-bias. Open circles show the mean $`\widehat{S}_3=\widehat{\overline{\xi }}_3/\widehat{\overline{\xi }}_2^2`$ in 10 realizations of the full box ($`L=300h^1\mathrm{Mpc}`$). Filled triangles show the mean $`\widehat{S}_3`$ over 10 realizations for each smaller box-size: $`L=300/2`$, $`300/4`$, $`300/5`$, $`300/7h^1\mathrm{Mpc}`$ as indicated in each panel. The squares, on the other hand, show the mean $`\widehat{\overline{\xi }}_3`$ divided by the mean $`\widehat{\overline{\xi }}_2^2`$, which is our way of isolating the integral-constraint-bias of $`\widehat{S}_3`$ (eq. ).

The solid line corresponds to the tree-level PT prediction for $`S_3`$. Its agreement with the open circles on large scales indicates that the measured skewness from the full box does not suffer from any appreciable estimation-bias, for the smoothing scales shown. The short-dashed line is our analytical prediction for the integral-constraint-bias of $`\widehat{S}_3`$ (eq.  i.e. no ratio-bias), and is in good agreement with the simulation results (squares). The long-dashed line is the net estimation-bias of $`\widehat{S}_3`$ (eq. ), which includes both the ratio-bias (eq. ) and the integral-constraint-bias (eq. ), and should therefore be compared with the triangles from the simulations. It can be seen that the ratio-bias of $`\widehat{S}_3`$ (eq. ) always dominates on large scales, and that the integral-constraint-bias of $`\widehat{S}_3`$ is negligible except for the smallest subsamples.

The agreement between our analytical prediction and the simulation results is good as long as the estimation-bias is not too large. Our analytical prediction breaks down when the bias becomes too large, as in the case of $`L=300/7h^1\mathrm{Mpc}`$. This is hardly surprising as the calculation was done by explicitly assuming that the estimation-bias is small. We have found a phenomenological fit to correct for this: instead of having $`\widehat{S}_3=S_3(1+\mathrm{\Delta }_{S_3}/S_3)`$, we use

$$\widehat{S}_3=S_3\mathrm{exp}(\mathrm{\Delta }_{S_3}/S_3)$$
(35)   
where $`\mathrm{\Delta }_{S_3}/S_3`$ is our linear estimate (in $`\overline{\xi }_2^L`$) for the fractional bias as before (eq. ). The higher order terms from the exponential helps partially cancel the over-prediction due to the linear term alone. This ansatz is shown as a dot-dashed line in the Fig. 3 & 4. As can be seen, there is a reasonable agreement with the simulation results of both models, indicating that this is an acceptable extrapolation. An alternative would be to go beyond linear order in $`\overline{\xi }_2^L`$, and compute the estimation-biases to second order. We will not attempt to do so here.

It should be emphasized that we have used the perturbation theory values for quantities such as $`c_{mm^{}}`$, $`S_3`$ and $`\overline{\xi }_2`$ in the analytical predictions for the various estimation-biases (eq. ,  & ). The good agreement on small scales between our analytical predictions and the numerical results above should be seen as somewhat fortuitous. Lacking actual measurements from simulations of the quantities $`c_{mm^{}}`$ on nonlinear scales, we will not attempt to do any better here. In practice, one might want to use improved determinations of $`c_{mm^{}}`$, etc at nonlinear scales, from simulations for example, in eq. (27) to (29), or even go back to their original formulations in eq. (3.1), (3.2) and (3.3). But there are a few reasons why our PT-based analytical predictions should work fairly well: 1. the terms involving $`\overline{\xi }_2^L/\overline{\xi }_2`$ are only important on large scales, and so using the PT values is adequate for these terms; 2. the true $`S_3`$ and $`c_{mm^{}}`$ change only slowly with the smoothing scale (i.e. the tree-order PT predictions are not too far off); 3. most of the terms involve ratios of $`c_{mm^{}}`$ and $`S_3`$, which are perhaps even slower functions of the smoothing scale.

## 5 Simulated Galaxy Catalogues and a First Step Towards Reconciliation

To make contact with existing measurements from actual galaxy catalogues, we add three elements of realism to our simulations: 1. use box-sizes similar to those of surveys where some of the measurements of the variance and skewness have been made; 2. introduce redshift distortions for redshift catalogues and projection for angular-catalogues; 3. allow for a realistic level of (sparse) sampling. Galaxy-biasing is not implemented here, however.

We simulate CfA/SSRS volume-limited catalogues similar to those studied in Gaztañaga Gaztañaga (1992), based on the LCDM model. Redshift distortions are modeled in the usual way, with the distant observer approximation. We consider two sets of subsamples, taken from each of the 10 full-box ($`L=300h^1\mathrm{Mpc}`$) realizations of the LCDM model as in Fig. 4. The first set has the same volume as the CfA/SSRS50 catalogues in Table 1 (a cubic box of $`L=40h^1\mathrm{Mpc}`$ on a side, with a volume equivalent to that of a sphere of radius $`R_L=25h^1\mathrm{Mpc}`$; the effective depth is $`𝒟50h^1\mathrm{Mpc}`$) and the second set has a volume similar to the CfA92 catalogue in Table 2 (a cubic box of $`L=78h^1\mathrm{Mpc}`$ on a side, with a volume equivalent to that of a sphere of radius $`R_L=48h^1\mathrm{Mpc}`$; the effective depth is $`𝒟92h^1\mathrm{Mpc}`$). We have checked that the results in this section are essentially unchanged if, instead of a cubical box, one considers a conical geometry which resembles that of the actual surveys. This is in part because of the actual solid angle substended by these surveys ($`1.6`$; see caption of Table 1).

Following Gaztañaga Gaztañaga (1992), we concentrate on counts-in-cells for spherical cells of radius $`R`$ in a range between $`422h^1\mathrm{Mpc}`$, depending on the subsample. The lower limit is chosen to avoid too much shot-noise, whereas the upper limit is picked to avoid large edge-effects. We vary the number of particles/galaxies in each catalogue to assess the effect of shot-noise on the estimation-biases.

Simulated angular catalogues that resemble the APM are taken from Gaztañaga & Bernardeau Gaztañaga & Bernardeau (1998). The power spectrum is that measured from the APM. We consider two different box-sizes $`L=378h^1\mathrm{Mpc}`$ and $`L=600h^1\mathrm{Mpc}`$.

### 5.1 The Variance $`\overline{\xi }_2`$ in Simulated CfA/SSRS Catalogues

Figure 5 shows the integral-constraint-bias for $`\widehat{\overline{\xi }}_2`$ from the simulated CfA/SSRS catalogues. The point-symbols correspond to the values measured from the simulated catalogues, while the lines correspond to either values measured from the full-box LCDM simulations (long- and short-dashed lines) or analytical predictions for the integral-constraint-bias (solid lines). The variance measured in the nearby sample (CfA/SSRS50) is significantly smaller than that in the deeper sample (CfA/SSRS90), at the smoothing scale $`R=9h^1\mathrm{Mpc}`$. (Note that for clarity, we do not show the measurements from the deeper sample at smaller scales, but they follow those from the full-box LCDM simulation rather closely.) This is due to the systematic bias introduced by the integral constraint, which generally gets worse for a smaller survey volume (eq. ; as shown also in Figure 2). This results in the phenomenon that the amplitude of the measured variance seems to increase with the sample depth, as found in several studies (e.g. Davis et al. (1988); Gaztañaga (1992); Bouchet et al. (1993)). Bouchet et al. Bouchet et al. (1993) correctly attributed (a significant part of<sup>3</sup><sup>3</sup>3Luminosity segregation also plays a role here.) this observed phenomenon to a finite-volume-effect along the lines of Colombi et al. Colombi et al. (1994). Seeing this as none other than the integral constraint allows us to predict the size of this bias analytically.

Note how, on large scales ($`R>9h^1\mathrm{Mpc}`$), the variance in redshift space (long-dashed line) is larger than that in real space (short-dashed line), as predicted by Kaiser Kaiser (1987), whereas the reverse holds on small scales, because of shell crossing and virialization (the finger-of-God). The analytical predictions (eq. \[3.1\]) for integral-constraint-bias of both catalogues are shown as solid lines. We show the predictions in real-space only, but it can be seen that the biased-estimates of the variance in real- and redshift-space are in fact quite similar – we will see this even more clearly for the skewness. As mentioned before, these predictions are only approximate for non-linear scales ($`R{}_{}{}^{}{}_{}{}^{<}\mathrm{\hspace{0.17em}8}h^1\mathrm{Mpc}`$), because we have not modeled properly the non-linear values of $`c_{mm^{}}`$ and $`S_J`$. Nevertheless there is an overall agreement between the simulation results and the predictions.

Lastly, we examine the effect of shot-noise on the estimation-bias, by sparsely sampling our catalogues (use 200 galaxies in each sub-sample instead of the $`10^4`$ in CfA/SSRS50 or $`10^5`$ in CfA/SSRS90). The effect can be seen to be small (compare closed triangles with open triangles). 

### 5.2 The Skewness $`S_3`$ in Simulated CfA/SSRS and APM catalogues

Figure 6 shows the results for $`S_3`$ in the simulated CfA/SSRS catalogues. Note how $`S_3`$ from N-body simulations is closer to the real space PT prediction (dotted line), when measured in redshift space (long-dashed line) than in real space (short-dashed line).

As before, we clearly see the variation of the estimation-bias with survey volume. The bias is more significant for CfA/SSRS50 than for CfA/SSRS90. The analytical prediction here comes from the phenomenological ansatz we introduce in §4.2 (eq. ). Only the analytical prediction for real-space is shown (solid lines). Note how the biased-measurements of the skewness yield very similar values in real- and redshift-space, even though the true $`S_3`$’s (i.e. measured from the full box; long and short-dashed lines) are quite different in the two cases, especially on small scales. In the next section, we will take advantage of this fact and attempt to perform a preliminary correction of some existing (biased) measurements of $`S_3`$ in redshift-space (which we assume to be very close to their real-space values), using our real-space analytical prediction for the estimation-bias. In other words, lacking analytical predictions for the redshift-space $`c_{mm^{}}`$, we use the real-space PT values for the various hierarchical coefficients in eq. (29) to make corrections for measurements that are actually done in redshift-space.

Figure 7 shows the measured 2D $`S_3`$ in simulations of APM-like angular catalogues. The results are a reproduction from Gaztañaga & Bernardeau Gaztañaga & Bernardeau (1998). Two box-sizes are considered. The large box with $`L=600h^1\mathrm{Mpc}`$ gives measurements which are in good agreement with the tree-level PT predictions on large scales, indicating that the estimation-bias is negligible in this case, at least for $`\theta {}_{}{}^{}{}_{}{}^{<}\mathrm{\hspace{0.17em}10}`$ deg., e.g. $`R{}_{}{}^{}{}_{}{}^{<}\mathrm{\hspace{0.17em}70}h^1\mathrm{Mpc}`$. Note that the actual APM survey has a size even bigger than this large box. The smaller box with $`L=378h^1\mathrm{Mpc}`$ results in a more appreciable estimation-bias,. Both of this results agree with our analytical prediction short-dashed line (long-dashed) line for the smaller (larger) box. The predictions can be obtained from eq.3.3 by just replacing the $`3D`$ hierarchical amplitudes by the $`2D`$ ones, as our derivation was totally general in this respect. We use the LCDM model PT predictions with the APM selection function and assume that $`c_{ij}^{2D}r_{i+j}c_{ij}^{3D}`$, with where $`r_{i+j}1`$ as given in Gaztañaga Gaztañaga (1994). These results are not very sensitive to the exact values of $`r_{i+j}`$.

### 5.3 A Preliminary Reconsideration of Some Existing Measurements from Actual Galaxy Surveys

#### 5.3.1 The Variance $`\overline{\xi }_2`$ in the Actual IRAS, CfA and SSRS Surveys

Any corrections of existing measurements of $`\xi _2`$ based on eq.(27) are necessarily model-dependent, because values for quantities such as $`c_{12}`$, $`\overline{\xi }_2^L`$ and even $`\overline{\xi }_2`$ itself need to be assumed. All of them vary with the amount and nature of galaxy-biasing. Instead of conducting a detailed analysis covering many possible models, we ask the following simpler question: assuming that the LCDM or SCDM model for the shape of $`\overline{\xi }_2`$ are the true ones, and assuming the PT-theory (real-space) values for $`c_{12}`$, what would be the corrected measurements of $`\overline{\xi }_2`$ for the IRAS, CfA and SSRS catalogues, using eq.(27), assuming no galaxy-biasing?

Fig. 8 shows the correlation length $`R_0`$, defined as $`\widehat{\overline{\xi }_2}(R_0)=1`$, as a function of $`R_L`$, the equivalent radius of the corresponding subsample. Open circles (squares) correspond to the values of the CfA (SSRS) volume limited subsamples at the bottom of Table 1. Filled squares correspond to the values in the IRAS 1.2 Jy volume limited subsamples by Bouchet etal (1993). The lines show the predictions for the measured $`R_0`$ taking into account the integral-constraint-bias for $`\widehat{\overline{\xi }_2}`$. We adopt the shape of the linear LCDM (continuous line) or the linear SCDM (dashed line) power spectrum to estimate all quantities, $`\overline{\xi }_2^L`$, $`\overline{\xi }_2(R)`$ and $`c_{12}`$, in eq.(27). This is only approximate as we do not take into account redshift distortions or non-linearities in the predictions, but note that in Fig. 5 we have found this approximation to be good.

The amplitude of the linear power spectrum is chosen to give the best fit to the data points in Fig. 8. The best fit values for the LCDM model are $`R_010.0h^1\mathrm{Mpc}`$ (where $`R_0`$ here is the $`R_0`$ we infer from the best-fit amplitude when matching our integral-constraint-prediction with the data points) for the CfA/SSRS (top continuous line) which has a joint $`\chi ^2=5.0/6`$, and $`R_06.5h^1\mathrm{Mpc}`$ for IRAS (bottom continuous line) which has a $`\chi ^2=37/20`$. The best fit values for the SCDM model are $`R_08.0h^1\mathrm{Mpc}`$ for the CfA/SSRS (top dashed line) which has a $`\chi ^2=11.4/6`$, and $`R_06.1h^1\mathrm{Mpc}`$ for IRAS (bottom dashed line) which has a $`\chi ^2=27/20`$. Note that the SCDM model gives a poorer fit to the CFA/SSRS data, while IRAS is compatible with both models. These values are to be compared with the mean values (which are usually taken to give the true amplitude): $`R_08h^1\mathrm{Mpc}`$ for the CfA/SSRS and $`R_05.5h^1\mathrm{Mpc}`$ for IRAS. Note that our correction for $`R_0`$ here ignores galaxy-biasing. Under such an assumption, the large value of $`R_010h^1\mathrm{Mpc}`$ (ie $`\sigma _81.21.3`$) found in the CfA/SSRS seems difficult to reconcile with the amplitude infered from the angular APM Galaxy catalogue: $`\sigma _8<1.08`$ (see Gaztañaga (1995)). To be strictly consistent, we should go back and allow for the effect of biasing on our correction for $`R_0`$: we would not attempt to do so here. Note also that the larger the volume limited subsample the brighter the absolute magnitudes of galaxies it contains; our correction for $`R_0`$ implicitly ignores luminosity segregation i.e. that the intrinsic clustering does not change significantly as a function of the absolute magnitude of the galaxies. A direct measurement in redshift space in the Stromlo-APM Catalogue gives $`\sigma _8=1.1\pm 0.1`$ for the brightest sample with $`M_{b_j}<20`$ (sample d. in Table 3 in Loveday et al.Loveday etal (1996)). This Stromlo-APM subsample contains galaxies with similar absolute magnitudes to those of the CfA92 sample, where $`M_{b_0}<20.3`$, given that $`b_0b_j0.3`$ Dalton & Gaztañaga (1998). Thus, our analysis indicates a small relative galaxy-bias between the CfA/SSRS and the APM galaxies that seems not attributable entirely to luminosity segregation.

#### 5.3.2 The Skewness in the Actual CfA, SSRS and APM Catalogues

Again here, it is clear that any corrections of existing measurements of $`S_3`$ based on eq. (29) are necessarily model-dependent, because values for quantities such as $`c_{mm^{}}`$, $`\overline{\xi }_2`$, $`\overline{\xi }_2^L`$ and even $`S_3`$ itself need to be assumed. All of them vary with the amount and nature of galaxy-biasing. Instead of conducting a detailed analysis covering many possible models, we ask the following simpler question: assuming different shapes for the power spectrum all normalized to $`\sigma _8=1`$, and assuming the PT-theory (real-space) values for $`c_{mm^{}}`$ and $`S_3`$, what would be the corrected measurements of $`S_3`$ for the CfA and SSRS catalogues, using eq. (29), or more appropriately, its extension in eq. (35), if no galaxy-bias is assumed?

Note that here, we are taking advantage of the finding in §5.2 (Fig. 6) that the biased-measurements of $`S_3`$ yield very similar values in real- and redshift-space, and so we can simply apply the correction formulated in real-space to the measurements in redshift-space. On large scales, this is a safe assumption because even the true $`S_3`$’s are very similar in real- and redshift-space. On small scales, this assumption remains to be further scrutinized. Obviously, this assumption must break down when the survey-size is large enough (c.f. short- and long-dashed lines in Fig. 6).

We will concentrate on the last six measurements in Table 1. These published values of $`S_3`$ (Gaztañaga (1992)) are the mean values in the corresponding range of smoothing scales shown in the 4th column. Here we will associate each with the mean value of $`R`$ in the respective range of scales i.e. $`R=5,8`$ & $`15h^1\mathrm{Mpc}`$ with $`S_3=1.8,1.7`$ & $`1.7`$ from the CfA50, CfA80 & CfA92 catalogues on the one hand, and $`S_3=1.4,1.9`$ & $`2.2`$ from the SSRS50, SSRS80 & SSRS115 catalogues on the other. The corrected values of $`S_3`$, adopting the assumptions stated above, are shown as lines in Fig. 9. The correction is larger for the smaller scales which correspond to sub-samples of a smaller size. The SCDM predictions are much lower as it has less power on large scales. 

If there is no biasing, we can assume that the APM-values for $`\overline{\xi }_2`$ is close to the true ones, which is is reasonable given the size of the APM (see e.g. Fig. 7). In this case we should concentrate on the continuous line in Fig. 9. At the larger smoothing scales, namely $`8`$ & $`15h^1\mathrm{Mpc}`$, the corrected values are quite consistent with the APM-values. However, at the smallest smoothing scale of $`5h^1\mathrm{Mpc}`$, the corrected $`S_3`$’s are significantly higher than the APM-values. Three points should be noted here. First, at this smoothing scale, the correction is so large that the validity of the ansatz expressed in eq. (35) might be called into question. Second, the implicit assumption that the biased- measurements of $`S_3`$ in real- and redshift-space yield similar values should be checked using more simulations. Third, it is in fact well known that the small-scale $`S_3`$ (of the mass) one would infer from an N-body simulation with an APM-like power spectrum and Gaussian initial conditions is larger than the measured $`S_3`$ from the APM survey Baugh & Gaztañaga  (1996). One possible interpretation is a scale-dependent galaxy-bias, which tends to diminish on large scales but becomes significant on small scales. 

On large scales $`R{}_{}{}^{}{}_{}{}^{>}8h^1\mathrm{Mpc}`$, we can safely say that most of the discrepancies in existing measurements of $`S_3`$ from the CfA/SSRS/APM catalogues can be explained by an estimation-bias. Remaining differences are attributable to a) a small relative galaxy-bias between the different surveys on large scales; b) redshift-distortions; c) deprojection effects (see Gaztañaga & Bernardeau (1998)) and d) sampling fluctuations. In fact, our analysis as shown in Fig. 8 does support the existence of a galaxy-bias between the CfA/SSRS and the APM. To be strictly consistent, we should have taken this into account in our “correction” of the CfA/SSRS values for $`S_3`$. Doing so is beyond the scope of the present paper. Nonetheless, it should be emphasized that the amount of galaxy-bias one would infer based on measurements of $`S_3`$ is reduced if the estimation-bias is taken into account. A careful assessment would require the inclusion of all the above effects, and ideally, a re-measurement of $`S_3`$ from the different catalogues using methods that are perhaps less prone to the ratio-bias. We hope to pursue these in a future paper.

## 6 Discussion

The main results of this paper are summarized in eq. (3.1), (3.2) & (3.3), with the associated useful approximations given in eq. (27), (28) & (29). Together they tell us the estimation-biases associated with the standard estimators for the variance $`\overline{\xi }_2`$, the third cumulant $`\overline{\xi }_3`$, and the skewness $`S_3`$ (eq.  & ). The calculation is based on an expansion in the small parameter $`\overline{\xi }_2^L`$ (eq. ), which is the variance smoothed on the scale of the survey, but otherwise does not assume or require the smallness of the variance on the scale of interest $`R`$.

¿From eq. (27), (28) & (29), it can be seen that the standard estimators are all asymptotically unbiased, in the sense that for a given smoothing scale $`R`$, the estimation-biases tend to zero as the survey-size increases. On the other hand, for a fixed survey-size, the estimation-biases become large as $`R`$ approaches the size of the survey.

There are two types of terms in the fractional estimation-biases, one dependent on the smoothing scale $`R`$, being proportional to $`\overline{\xi }_2^L/\overline{\xi }_2`$ where $`\overline{\xi }_2`$ is the variance smoothed on scale $`R`$, and the other not, being proportional to $`\overline{\xi }_2^L`$ only. In other words, a general form for the estimation bias of an estimator $`\widehat{E}`$ can be represented by eq. (32) where $`\alpha _1`$ and $`\alpha _2`$ are coefficients that depend on the various hierarchical amplitudes, such as $`S_3`$ and $`c_{mm^{}}`$. For reasonable choices of the parameters $`S_3`$ and $`c_{mm^{}}`$, the estimation-biases are negative (see Fig. 1). The magnitude of the estimation-biases can be surprisingly large, especially for $`E=S_3`$, because of the large coefficients multiplying $`\overline{\xi }_2^L`$ or $`\overline{\xi }_2^L/\overline{\xi }_2`$.

Our analytical predictions are borne out by numerical experiments discussed in §4. Examples can be found in Fig. 2, 3 & 4. In cases where the estimation-biases are so large that the higher order terms in our expansion in $`\overline{\xi }_2^L`$ become important, we show that a simple ansatz works reasonably well (eq. ).

In the case of $`S_3`$, we distinguish between two types of contributions to its estimation-bias. The standard estimator for $`S_3`$ is $`\widehat{S}_3=\widehat{\overline{\xi }_3}/\widehat{\overline{\xi }_2}^2`$. Part of the bias arises from the biases of the estimators $`\widehat{\overline{\xi }_3}`$ and $`\widehat{\overline{\xi }_2}`$ themselves – this is the integral-constraint-bias (eq. ). The second part arises from the particular nonlinear combination of these two estimators, which we dub to be the ratio-bias (eq. ), or more precisely, the nonlinear-estimation-bias. It turns out the second always dominates on large scales.

We present a preliminary attempt to correct some existing measurements of the variance and skewness in §5. Our main conclusions are a) the apparent increase of the measured variance with survey depth observed by some authors can be nicely explained by the integral-constraint bias (Fig. 8; see e.g. Davis et al. (1988); Gaztañaga (1992); Bouchet et al. (1993)); b) our analysis indicates a small relative galaxy-bias between the CfA/SSRS and the APM galaxies, c) the APM survey should give small estimation-biases for the standard estimator for $`S_3`$ on scales of interest (see Table 1); d) on large scales $`R{}_{}{}^{}{}_{}{}^{>}8h^1\mathrm{Mpc}`$, most of the differences between the measured skewness from the CfA/SSRS and the APM surveys can be attributed to the estimation-bias; however, a more careful analysis, taking into account redshift-space distortions, relative galaxy-bias (such as due to luminosity segregation) and deprojection effects, is necessary to access the significance of the remaining differences, and whether they are due to pure sampling fluctuations.

Two areas clearly warrant further investigations. First, for the purpose of future surveys such as the SDSS or the AAT 2dF, the $`\overline{\xi }_2`$-independent terms in the estimation biases (the terms associated with $`\alpha _2`$ in eq. ) will probably be unimportant, but the $`\overline{\xi }_2`$-dependent terms (the $`\alpha _1`$-terms) can always become important for a sufficiently large smoothing scale $`R`$. <sup>4</sup><sup>4</sup>4However, for most current models of the power spectrum, as $`R`$ increases $`\gamma `$ becomes more negative and the coefficient $`\alpha _2`$ in $`S_3`$ becomes smaller (see Fig. 1). It would therefore be good to have an idea of what that scale is, not just for the variance and the skewness, but also for $`S_N`$ where $`N>3`$. We will present results for a general N in a separate paper.

Second, we have focused in this paper exclusively on the standard estimators for $`\overline{\xi }_N`$ and $`S_N`$, where $`\overline{\xi }_N`$ is estimated by the standard counts-in-cells technique and $`S_N`$ is estimated by taking the appropriate ratio of $`\widehat{\overline{\xi }}_N`$ and $`\widehat{\overline{\xi }}_2`$. There are probably other estimators that suffer from smaller estimation-biases. For instance, Kim & Strauss Kim & Strauss (1998) recently introduced an interesting method to obtain $`S_3`$ by fitting the one-point probability distribution of counts (PDF) using an Edgeworth series. They obtained values higher than previous measurements from the same surveys (see Table 1) indicating that their method is less susceptible to an estimation-bias. An extension of their method using a PDF which is better behaved in the presence of more significant nonlinearities is worth pursuing.

Another possibility which might reduce the ratio-bias in $`\widehat{S}_N`$ due to the division of $`\widehat{\overline{\xi }}_N`$ by $`\widehat{\overline{\xi }}_2^{N1}`$: instead of dividing to obtain an estimate of $`S_N`$, fit a curve parametrized in some form to the two-dimensional plot of $`\widehat{\overline{\xi }}_N(R)`$ and $`\widehat{\overline{\xi }}_2(R)`$ at each smoothing scale $`R`$. This is akin to, for instance, how the Hubble constant is usually measured: instead of dividing some estimate of velocity by some estimate of distance for each data point followed by averaging, a linear $`\chi ^2`$ fit to all points is performed. This method might give a smaller ratio-bias, but it has to be tested. This procedure has in fact been carried out before by e.g. Gaztañaga Gaztañaga (1992) and Bouchet et al. Bouchet et al. (1993), who obtained values of $`S_3`$ close to that using the standard estimator. But a different parametrization of the relation between $`\widehat{\overline{\xi }}_N(R)`$ and $`\widehat{\overline{\xi }}_2(R)`$ (in other words, taking into account carefully the change of $`S_N`$ with $`R`$) might yield different results.

Perhaps the most important and obvious lesson of our investigation here is that nonlinear combinations of estimators should be used with caution. The only sure-fire way of avoiding an estimation-bias is through Monte-Carlo simulations, such as those performed in §4. There is nothing novel about this point, except that its importance has not been sufficiently emphasized in measurements of certain large scale structure statistics, such as $`\overline{\xi }_N`$ and $`S_N`$ discussed here.

There are also strong interests in measuring such quantities outside galaxy surveys, e.g. for the transmission distribution of quasar spectra and for the convergence distribution in weak-lensing maps. The skewness of the former, for instance, provides a test of the gravitational-instability picture of the Lyman-$`\alpha `$ forest (Hui (1998)), while the skewness of the latter is a sensitive probe of cosmology (Bernardeau et al. (1997)). Measurements in these areas require similar caution as in the case of galaxy surveys. Our methodology in deriving the estimation-bias for the standard estimator of skewness (ratio of the third cumulant to the second cumulant squared) could be adapted for such measurements (in the case of weak lensing, we learned as this work was being completed that a related calculation was done by Bernardeau et al. (1997)). An important ingredient of our calculation is the use of the hierarchical relation $`\xi _N\xi _2^{N1}`$ in keeping track of the ordering, which might have to be modified for these other applications.

Also, for these applications, as well as for less conventional galaxy surveys such as the Lyman-break galaxy surveys at high redshifts (Steidel et al. (1998)), one often has available several independent fields for which a simple and obvious method should help to reduce much of the ratio-bias of the skewness: measure the third moment and the second moment for each field, separately average each of them over all fields, and only then does one combine these averaged moments to estimate the skewness.

A related suspect of a similar estimation-bias is the measurement of $`Q_N`$, defined as the ratio of the N-point correlation to the sum of suitable permutations of products of the two-point functions (Fry (1984)). The configuration dependence of $`Q_3`$, for instance, provides an elegant test of the galaxy-bias (Fry (1994)). Common ways of estimating it, where estimators are divided by each other, are susceptible to the ratio-bias just as in the case of $`S_N`$. Another possibly problematic statistic is the ratio of the quadrupole to monopole power in redshift-space, or the ratio of the monopole power in redshift-space to that in real-space, which is often used to estimate the parameter $`\beta =\mathrm{\Omega }^{0.6}/b`$ (see review by Hamilton (1997)). An examination of published estimates of $`\beta `$ show a large scatter even from the same surveys, with maximum likelihood methods yielding consistently higher values (Hamilton (1997)), suggesting an estimation-bias of some sort might be lurking here. It is likely, however, that such attempts to measure $`\beta `$ are at least equally, if not more strongly, affected by our poor understanding of translinear distortions. We hope to pursue some of the above issues in future work.

We are grateful to Joshua Frieman for useful discussions. This work was in part supported by the DOE and the NASA grant NAG 5-7092 at Fermilab, and by a NATO Collaborative Research Grants Programme CRG970144 between Fermilab and IEEC. EG acknowledges support by IEEC/CSIC, and by DGES (MEC), project PB96-0925.