Estimation of the Covariance Function of Gaussian Isotropic Random Fields on Spheres, related Rosenblatt-type Distributions and the Cosmic Variance Problem

We consider the problem of estimating the covariance function of an isotropic Gaussian stochastic (cid:133)eld on the unit sphere using a single observation at each point of the discretized sphere. The spatial estimator of the covariance function is expressed in a new form which provides, on one hand a way to derive the characteristic function of the estimator, and on the other hand a computationally e¢cient method to do so. We also describe a methodology for handling the presence of the cosmic variance which can impair the results. In simulation, we use the pixelization scheme HEALPix. 12


Introduction
This paper is about the statistical analysis of a Gaussian isotropic spherical random …eld T (x) on the unit sphere S 2 = fx 2 R 3 : kxk = 1g in Euclidean space R 3 , when only one observation of the …eld is available.This perspective is relevant for the analysis of the Cosmic Microwave Background Radiation (CMBR) discovered by the astronomers Arno Penzias and Robert Wilson in 1964.It is due to the emission of black body thermal energy originating from the big bang.The spectral radiation is measured at di¤erent angles of observation of the sky, see [65].It is apparently almost isotropic.
Our goal is to estimate the covariance function of the random …eld T (x), x 2 S 2 , in a parametric setting, given a single observation at each point of the discretized sphere.As application, using a result of Veillette-Taqqu [61], we present a methodology for handling the problem of cosmic variance in this framework.The cosmic variance is de…ned in (28) below.It results from uncertainty due to the fact that one has only a single observation.
It is well known that the second order structure, i.e., either covariance or spectrum completely characterize a centered Gaussian random …eld.Therefore the estimation of these quantities is of primary importance.The estimation of the spectrum is a well-studied subject [41] [14], [13] and so is the estimation of the covariance.If the estimation is non-parametric then the cosmic variance, de…ned in (28), will prevent us in getting a good estimator unless the spectrum at low frequencies vanishes, which would be unusual.If we are dealing with a parametric problem, however, that is, if the covariance function depends on some unknown parameters, then there is a chance of getting a reasonable estimator, see [40] as well.The method would be as follows.Given observations of the random …eld T , estimate its covariance function nonparametrically and then its spectrum.Set the low frequencies to zero.This yields a modi…ed estimated spectrum.Then estimate the unknown parameters by minimizing the sum of squares of the di¤erence between the theoretical form of the spectrum and the modi…ed estimated spectrum.As indicated below, one can alleviate in this way the cosmic variance problem.On the other hand the theory underlying the estimation of the parameters of the angular spectrum using wavelets in a framework of higher frequency asymptotics can be found a number of papers, see [12], [16], [17], [15] [11], [25], [5], [5], [39], see also their references.
The paper is organized as follows.In Section 2, we review some basic notions related to isotropic random …elds, and we describe di¤erent models in Section 3. In Section 4 we focus on the covariance function C (cos ) which is a function of the angle between two points on the sphere.We estimate this covariance function using empirical covariances b C (cos ) based on a single observation at each point of the discretized sphere.The characteristic function of these empirical covariances is given in terms of cumulants.It turns out that our estimator follows a Rosenblatt type distribution for each given single angle .In Section 5 we focus on the di¤erence R = b C (cos ) C (cos ).Following results of [61], we obtain the distribution of R and related properties.In Section 6 we discuss the problem of cosmic variance, namely the e¤ects of the uncertainty due to the fact that only one realization is observed.To alleviate this e¤ect one can approximate R by R M which does not involve the low frequencies and for which the cosmic variance is negligible.We show that R M tends to a Gaussian distribution as M ! 1.In Section 7 we provide simulations using HEALPix which is a high level pixelization of the sphere S 2 and show that M as low as 4 can su¢ce.Our results can be generalized to higher dimensions d 3, by using Gegenbauer (ultraspherical) polynomials (C n ) instead of the Legendre ones (C 1=2 `= P `).Theorems 2 and 4 are of particular interest.Theorem 2 provides the characteristic function of b C (cos ).It is a theoretical result, but with a clear statistical meaning, since it speci…es the distribution of the empirical covariance.For example, one could estimate the unknown parameters of the covariance function of a spherical random …eld using nonlinear regression, and thus having information about errors is useful when applying the existing methodology.Such information would also be needed when testing hypothesis on the unknown parameters of the nonlinear regression.Theorem 4 is important because it gives a normal approximation for the tail R M .
Section 7 contains a conclusion.An appendix contains examples, a brief description of white noise analysis on the sphere, the Thorin class and measure and formulae used in the paper.

Preliminaries
Let ( ; z; P) be a probability space, and S 2 = x 2 R 3 : kxk = 1 be the unit sphere centered at the origin.We consider a real-valued random …eld T (!; x) = T (x), ! 2 ; x 2 S 2 , with ET (x) = 0.This random …eld is said to be second-order weakly isotropic or (simply) isotropic, if ET (x) 2 < 1; and ET (x) T (y) =ET (gx) T (gy) for any g 2 SO (3), x; y 2 S 2 , where SO (3) denotes the three dimensional rotational group under composition.The orthogonal system on S 2 is given by the complex-valued spherical harmonics Y m `, where `= 0; 1; 2; : : :, and m = `; `+1; : : : 1; 0; 1; : : : ; ` 1; `.Their expression is given in (50).The Euler angles (#; ') de…ne the position x (#; ') = (sin # cos '; sin # sin '; cos #) of a point on the sphere with colatitude3 # 2 [0; ] and longitude ' 2 [0; 2 ].The colatitude measures the north-south position and the longitude the east-west position.We suppose that T (x) is mean square continuous, and hence it admits a series expansion ( [41]), in terms of the complex-valued spherical harmonics Y m `, with coe¢cients given by where (dx) = sin #d#d' is the Lebesgue measure of surface area on S 2 , and where star denotes the complex conjugate.The series (1) converges in L 2 ( ; R) for all x 2 S 2 : If the coe¢cients a `m are independent and for …xed `are identically distributed then the covariance function C 2 (x 1 ; x 2 ) =ET (x 1 ) T (x 2 ), (ET (x) = 0) depends on the angular distance between x 1 and x 2 only.This angle results from the inner product x 1 x 2 = cos .The covariance function depends on this central angle between locations and has the form where P `denotes the Legendre polynomial, see (49) The coe¢cient f `in (3) de…nes the angular spectrum and satis…es f ` 0, see [41], [66].We assume …nite variance, and since the Legendre polynomials are bounded, jP `(y)j 1, and P `(1) = 1; we get The fact that C 2 (x 1 ; x 2 ) = C (cos ) indicates that C 2 (x 1 ; x 2 ) is invariant under the group of rotations.The random …eld T (x) said to be linear if a `m are independent and if for …xed `, they are identically distributed.We work with Gaussian random …elds, they happen to be linear and linear …elds are automatically Gaussian, see [6], [58].
From now on we assume Assumption: T (x) is Gaussian, with …nite variance.
We can obtain the angular spectrum f `from the covariance through the relation For a given T (x) we have the inversion (2) and a `; m = ( 1) m a `m, since the …eld T (x) is real-valued and since Y m `= ( 1) m Y m `.The orthogonal random 'measure' a `m is a triangular array, we have m = `; `+ 1 : : : ; ` 1; `, i.e. rows contain 2`+ 1, i.i.d Gaussian random variables a `m with Ea `m = 0; Ea `ma kn = f ` `;k m;n : In particular a `m is normal with mean 0 and variance E ja `mj 2 = f `.
Remark 1 A function de…ned by (3) with the coe¢cients f `is strictly positive de…nite if and only if f `is strictly positive for in…nitely many even and in…nitely many odd integers `, see [26] for details.
For instance a class of covariance functions on spheres can originate from covariance functions of some homogenous and isotropic random …elds on Euclidean spaces since the restriction of the …eld to the sphere yields an isotropic …eld on the sphere.In this case consider two locations x 1 and x 2 on the sphere with angle 2 [0; ].Then the distance r = kx 1 x 2 k between them expressed in terms of the angle is 2 sin ( =2), see Figure 1, and the inner product is x 1 x 2 = cos , which gives a direct correspondence between the original covariance function C 0 (r), in the Euclidean space and the covariance function on the sphere.This holds for any dimension of the Euclidean space.The disadvantage of using C 0 is that it depends on the chordal distance between two points on the sphere instead of the grand-circle (spherical or geodesic) distance, which is not practical.
One can consider a more natural Laplacian model de…ned directly on the sphere S 2 when the distance is measured using the grand-circle distance.
Example 1 Laplace-Beltrami model on S 2 .We consider the stochastic model on the sphere S 2 for an isotropic random …eld T B on S 2 , (the index B is for "Beltrami") satisfying the equation in the L 2 sense, where @W B denotes the white noise with variance 2 , see the Appendix A for the de…nition of @W B .The Laplace-Beltrami operator is A direct calculation leads to the spectrum for T B , and the covariance function is given by formula (3).This form of covariance function obtained via white-noise-driven damped di¤usion equations for modeling global temperature …elds by [43], see also [31].The rigorous theory can be developed in the same line as it is done in [32] see again Appendix A for more details and references.
The methodology described in this paper applies to some more examples, see the Appendix. .

Empirical covariances
We have de…ned C (cos ) in (3), and now we suppose that an observation of the …eld T (x), is given on the whole unit sphere S 2 and ET (x) = 0.
Consider a location x on the sphere S 2 and let an angle 2 [0; ] be given.Consider all locations y with angular distance to x, so that x y = cos .Locations y form a circle C (x; ) with center x and radius sin , see Figure 1.Now de…ne a rotation g (x; ) 2 SO (3) which rotates the sphere S 2 around x by an angle .The point x being the center will not be moved but any location y on the circle C (x; ) will be moved to some new location denoted y (x; ) = g (x; ) e y .The y (x; ) has the property x y (x; ) = cos since the rotation preserves the angular distance between two points.
In practice the data T (x) is discretized, for instance when T (x) measures the Cosmic Microwave Background anisotropies, the measurements are given on a high resolution pixel structure called HEALPix of the sphere S 2 and therefore (8) can be approximated with high precision.The calculation of (8) involves summation of products of the data as is the usual case for covariance estimators.The usual estimator of the covariances used in cosmology, due to [48], involves estimating the spectrum (through Ea 2 `m = f `) …rst, then using (3) next.We shall use a discretized version of (8) for actual computation of the estimate but formula (9) in the next Theorem will be used to obtain the distribution of the estimator.
where b a `m are independent and identically distributed normal random variables with Eb a `m = 0; and E jb a `mj 2 = f `: Proof.We denote the North pole N = (0; 0; 1) since it is at colatitude # = 0 and longitude ' = 0 and since the radius equals 1.For each location x one can …nd a rotation g such that x = gN , that is it maps the North pole to x.The inverse g 1 of the rotation g does not change the angular distance between two points hence g 1 x g 1 y (x; ) = cos .The rotation g 1 maps x to the North pole g 1 x = N , and the circle C (x; ) to the circle C (N; ).The points on that circle are g 1 y (x; ) = z (N; ) = (sin cos ; sin sin ; cos ), 2 [0; 2 ].Now in (8), the integral on C (x; ) becomes an integral from 0 to 2 and the integral on S 2 becomes to the integral on SO (3) according to the Haar measure, see (46), so We apply the series expansion (1) to both T (x) = T (gN ) and T (gz (N; )), where b a `m = a `m, are calculated in (2), they are independent and identically distributed normal random variables with Eb a `m = 0; E jb a `mj 2 = f `; (10) and D
The next Theorem gives the marginal and joint characteristic function of b C (cos ).See also [48].8) has the form (9) with characteristic function Let m 2 [0; ], m = 1; 2; : : : j, be given angles, then the joint characteristic function of and notice that the coe¢cients jb a `mj 2 f `are Hermite polynomials of degree 2, see Appendix, Formulae for details.Let H 2 (b a `m) denote jb a `mj 2 f `= jb a `mj 2 E jb a `mj 2 , for simplicity and rewrite ( 14) in terms of Hermite polynomials We now use the cumulants of the Hermite polynomials (see (44) and ( 45)), in particular the variance Var b In formula (15), all H 2 (b a `m) are independent for all `= 0; 1; 2; : : :, and m = 0; 1; : : `, hence we obtain from ( 15) The convergence of ( 16) follows from (4), i.e. f `< o ` 1 , and from the fact that P `is bounded by one for any `.Similarly for general k, we use the higher order cumulants ( 44) and ( 45) of the Hermite polynomials and obtain Hence by ( 15) The cumulant characteristic function of b C (cos ) follows: (2`+ 1) Now, from the identity we obtain (2`+ 1) ln (1 izf `P`( cos ) =2 ) ; which leads to Consider now the joint cumulant.Using the relation Cum 2 (aX; bY ) = ab Cum 2 (X; Y ) and ( 17), we get (2`+ 1) f 2 `P`( cos 1 ) P `(cos 2 ) : (2`+ 1) In general the joint cumulant is given by hence the characteristic function is ' (z 1 ; z 2 ; : : : 4 Distribution of the error We shall focus here on This is the error we make by using b C (cos ) instead of C (cos ).Recall that E b C (cos ) = C (cos ) by (13).
Theorem 3 The random variable R in (20) can be represented as where d = means equality in distribution and where U 2`+1 = (2`+ 1) is Gamma distributed with parameters (2`+ 1) =2 and 2= (2`+ 1).The characteristic function of R is = exp which is a Rosenblatt type characteristic function.It is in…nitely divisible and selfdecomposable with Lévy density (2`+ 1) exp x 8 f `P`( cos ) : Moreover ' (z) belongs to the Thorin class T (R), with Thorin measure given by (2`+ 1) 1=b `(x) ; where b `= 8 f `P`( cos ) : Proof.The characteristic function ( 22) and ( 23) of R follows from (18) and (12).Now, rewrite R given in (14) in the following form Since b a `0 is real therefore jb a `0j 2 =f `has 2 1 distribution 2 jb a `mj 2 =f `has 2 2 distribution and they are independent.A simple consequence of this is that the random variables 2`+1 (2`+ 1) distributed and independent.Hence the characteristic function of R can be expressed as Hence a Rosenblatt type characteristic function [46], [47] shows up, see also [62] and [34].
Let ID(R); SD(R) be the classes of in…nitely divisible and selfdecomposable distributions correspondingly.We next de…ne the Thorin class on R (see [59], [9], [30], [37]) as follows: We refer to the product u as an elementary gamma random variable if u is nonrandom non-zero vector in R; and is a gamma random variable on R + : Then, the Thorin class on R (or the class of extended generalized gamma convolutions), denoted by T (R); is de…ned as the smallest class of distributions that contains all elementary gamma distributions on R; and is closed under convolution and weak convergence.It is known that and inclusions are strict, [30].Since any selfdecomposable distribution on R is absolutely continuous (see, for instance, Example 27.8 of [50]) and is unimodal (by [67], see also Theorem 53.1 of [50]), then, any selfdecomposable distribution has a bounded density function.Thus the distribution with characteristic function (23) has a bounded unimodal density.Also (see Leonenko et al. [34] for details ) that the distribution with characteristic function (23) belongs to the Thorin class T (R), with Thorin measure given by where b `is given in the statement of the theorem.
Remark 2 Theorem 3 shows the similarities and di¤erences between the behavior of the estimation errors of an unknown covariance function for isotropic random …eld and the results for stochastic processes or time series, in which only asymptotic distributions are known.Surprisingly in our case we can obtain the explicit distribution, in terms of characteristic function, of the approximation error and even the rate of convergence.

Dealing with the cosmic variance problem
Consider a sample path of the …eld All information contained in this single sample path about the coe¢cients f `in the series expansion of the covariance function C (cos ), see (3), is expressed through the random variables a `m.Although a `m can be inverted with high precision for every indices `, m, see (2), for small 'frequencies' `, a `m has little information useful for estimation.For instance if `= 0, b f 0 = jb a 00 j 2 is a single value realization of a 00 which tells almost nothing about f 0 = E ja 00 j 2 .If `is large then we have b a `m, m = `; `+1 : : : ; ` 1; `, i.e. a 2`+1 'sample' for estimating and since b a `m is normal with mean 0 and variance f `.We can now de…ne the cosmic variance as see [64], p. 138.This quantity does not depend on the actual values of the spectrum and is decreasing with `.It underlines the uncertainty of statistical methods associated with the estimation of either the spectrum or the covariance function.Therefore reducing the cosmic variance is of primary importance.How to decrease the in ‡uence of the cosmic variance?Since the main contribution to that variance comes from f `with small values of `, we should try to ignore these f `by truncating the di¤erence R = b C (cos ) C (cos ), given in ( 20) and ( 21).Consider then the case when f `= 0, `= 0; 1 : : : ; M 1, in R, see (24), i.e. the remainder where the sum starts at `= M .Since R M is associated to the sample path T M (x) de…ned as and since b a `m are good approximations of the current values which are generating the observed random …eld T (x), see ( 2), (not like the estimation of f `), therefore T M (x) is a good approximation to the remainder …eld with f `= 0, `= 0; 1 : : : ; M 1, and with covariance function f `(2`+ 1) P `(cos ) : The asymptotic distribution and Berry-Esseen bound for the remainder R M given in the next Theorem can be obtained as in Theorem 3.1 of [61].
where is the standard normal CDF, and 3;M denotes the third cumulant (skewness) of R M Proof.The theorem follows from the Theorem 3.1 of [61], provided one has where `and r `are de…ned in (26) where in the nominator we applied the inequality (x + y) a x a + y a .valid for any 0 x; y 1, and any 0 < a < 1. Summarizing these we obtain Edgeworth expansion for the distribution function of R M is also given in [61] and the assumptions are satis…ed in our case as well.We will not include them here.Our simulation and the numerical example of [61] show that the speed of convergence is really fast and M can be chosen to be larger than 5, which is satisfactory for cosmology ( [64], p. 138) and as we shall see in the next section our estimator of the covariance function gives very good results even when M = 4.An earlier attempt in this direction has been made in [8].

Simulations
When dealing in practice with random …elds, we do not use for b C (cos ) its expression (8) or equivalently (9) instead we use a discretization of integral (8), namely (31) below, which corresponds to the time domain estimator of covariances in time series analysis.We consider a discretized unit sphere S 2 .The discretization, called HEALPix (Hierarchical, Equal Area and iso Latitude Pixelization), is applied.For a detailed description see [27].This pixelization of the sphere contains quadrilaterals (pixels), in our case the total number of equal-area spherical pixels equals to N pix = 49152; with area pix = 4 =N pix , since 4 is the surface of the unit sphere.
The integral for …xed angle where x j denote locations of pixel centers, T = (1=N pix ) P x T (x), is the mean and n x is the number of all possible pairs of x and x j , such that x x j = cos .Hence the covariance estimator is b Since is the angular distance and for a given location x all the locations with angular distance constitute a circle, in practice instead of a circle we considered a ring with a very narrow belt so that (30) contains all the pixel centers from this belt.The reason is that the pixel structure provides some speci…c angular distances only and we collect all the pixel centers having angular distance close enough to .One may consult with [4], [33] for properties of the above approximation.
From now on we shall scale the covariance function such that C (cos 0) = C (1) = 1 resulting in the correlation function.We do so because the parameter estimation requires that there be a unique correspondence between the model and its parameters.Thus no multiplicative constant will appear in the parameter estimation.
For simulation purposes, following (1), we generated a truncated …eld where K = 42 , and a `m, `= 0; 1; 2; : : : ; K, m = 0; 1; 2; : : : ; `, complex-valued Gaussian random numbers a `m, Ea `m = 0; Ea `ma kn = f ` `;k m;n , such that a `; m = ( 1) m a `m.The spectrum f `, `= 0; 1; 2; : : : ; K, is calculated using the Laplace-Beltrami model (Example 1), with parameters 2 = 2, c = 2.We get f `by applying (7).For all pixel centers x i on the sphere we generated the …eld T (K) (x i ), i = 1; 2; : : : ; N pix using (32).Now we consider the covariance function of the model with the values of f `truncated up to L, and we are going to estimate C (L) (cos ) using T (K) (x).The correspondence between the function C (L) (cos ) and the spectrum ff `gL 0 is given by the integral The exact value of that integral can be calculated via the Gauss-Legendre quadrature, as follows where the nodes y 1 ; : : : ; y L+1 are the roots of the Legendre polynomial P L+1 (x), while w 1 ; : : : ; w L+1 are the corresponding weights of the formula.In this case the quadrature is exact for polynomials up to order 2L + 1, [44], [56].Note that the order (highest degree of the polynomial) of C (L) (y) P `(y) is not larger then 2L, for any `.Given real data, L is not known, one chooses L and L + 1 angles 1 ; : : : ; L+1 , with L large enough so as to ensure that the estimator b C (cos ) in ( 31) provides good results.The number of terms in the second summation of (31) depends on the angular distance and if this number of terms is too small, one may end up with a bad estimate for C (cos ).The number L can therefore be considered as a bandwidth in this estimation.We used here angles 1 ; : : : ; L+1 corresponding to the roots of the Legendre polynomial P L+1 (x).Here we set L = 42.
After estimating the covariance function namely, by replacing C (L) in ( 35) by b C (y i ) obtained by using (31).The Laplace-Beltrami model with given 2 depends on an unknown parameter c, see Example 1 and has spectrum f `(c) given in (7).We estimated the parameter c in two steps.
First we used the estimated covariance function b C (cos ) in (31) to estimate the spectra b f by (36).Then we …tted f `= f `(c) to b f `, `= 1; 2; : : : ; L = 42, by the nonlinear least squares method, i.e. minimizing the and derived the estimate b c.Secondly, we obtained b c M with M = 4, in order to reduce the cosmic variance.To do so we estimate b a `m from T (K) (x) and remove the corresponding terms by We get a new model where a `m = 0, `= 0; 1; 2; 3, m = 0; 1; 2; : : : ; `, without modifying the other a `m.This yields the sample path with spectrum f `(c) = 0, `= 0; 1; 2; 3. Now we re-estimated the correlations b C (cos k ) in (31) replacing the …eld T by T (K;4) .We obtain a far better estimate of c because the cosmic variance has been reduced, setting its values corresponding to `= 0; 1; 2; 3, to zero.To get b c M , we repeated the least squares estimation by setting f `= 0, for `= 0; 1; 2; 3, and f `= f `(b c M ), for `= M; M + 1; : : : ; 42.Putting the estimated value b c M in (7) we get new values f `(b c M ), `= 0; 1; 2; : : : ; L. We use them in (33) to obtain a new estimated b C (cos ).We did 100 iterations.In each iteration we computed b c and b c M , then b f `using (7) and then the corresponding estimated covariance curve using (33).The results then were averaged including the covariance curves.Note that the displayed Figures 2, 3, 4 have di¤erent vertical scales.
The true value of c was 2 and the average of b c over 100 iterations was 1:1866 with standard deviation 0:5329, while the average of b c M was 1:8063 with standard deviation 1:3672.The average of b c M is closer to the true value then the average of b c though variance.
The plot in Figure 2 contains from top to bottom the following covariance functions obtained as follows: (a) estimated using (33) with b c, (b) estimated using (33) with b c M , (c) theoretical using (33) with c = 2, (d) estimated using (31).
Note that the curve (b) using b c M is the closest to the theoretical one (c). 4 The plot in Figure 3 shows the theoretical covariance function (continuous blue line with asterisks, computed from the spectrum (7) of the Laplace-Beltrami model using ( 33)), the average of the 100 estimations of the covariances b C (cos k ) (red dashed line with circles), upper and lower 95% con…dence intervals (red dashed line with asterisks), each of the 100 estimations b C (cos k ) (black points).It is seen that even the average of 100 estimations is not a good estimate of the covariance function mainly because of the cosmic variance.We obtain better results using b c as seen in Figure 4.
We conclude from Figure 4 that using the model T (K;4) (x) with reduced cosmic variance gives good estimates for the corresponding correlation function.As a consequence the updated estimate of the original covariance function by b c M provides a better estimator not only for the covariance function but for the spectrum as well.Figure 2 yields the same conclusion.easier to estimate the covariance function if low frequencies of the spectrum are zero. 4Removing the sample mean from the data in (30) results in ramoing a00Y 0 0 .Hence the estimated covariance function does not contains f0, that is b f0 = 0.It is true that a00Y 0 0 is random and not a constant, since a00 is Gaussian, but for real data we have only a single realization which means that we have a single value of a Gaussian variable.Thus, after estimating b c we added f0 (b c) to (33).In other words we proceeded as in the estimation of b cM but for M = 0.

Conclusion
We considered the problem of estimating the covariance function of an isotropic …eld on the sphere, at a …xed time as is the case of CMB data for instance.We derived the distributional properties of the nonparametric estimator of the covariance function.
The problem of estimating either the correlation function or the spectrum for the CMB data has a wide literature [18].Most of the methods considered, like Pseudo-C `estimators, NRML (maximum-likelihood using Newton-Raphson algorithm), QML (quadratic ML) hybrid estimator, su¤er from cosmic variance.The paper [19] considers estimates of the correlation function based on methods used in [18] paying attention to the cosmic variance as well.The aim of those investigations include checking Gaussianity, isotropy, modeling using six-parameter in ‡ationary CDM cosmology etc.There is a common agreement that "the analytic approximations at low multipoles are useless for any quantitative application such as parameter estimation" ( [18]).One of our aims was to reduce the cosmic variance in parametric models.
Theorems 2, 3 and 4 are connected.They describe and quantify the distribution of errors between the empirical covariance and the theoretical covariance function of spherical random …elds.Theorem 4 is of interest in statistical inference since it quanti…es the errors of approximation after truncation.It is known that the cosmic variance prevents us to getting a good estimator mainly because of the problem of estimating the spectrum f `at low frequencies, `= 0; 1; : : : ; M 1, say.Since the cosmic variance a¤ects mainly the spectral coe¢cients f àt low values of `, we ignore these f `setting them to 0. We then reestimate the covariance function and use it to estimate the unknown parameters of our parametric model, for example, the parameter c in the Laplace-Beltrami model.
The steps described above change the model, but this modi…ed model now provides better estimates of the covariance hence, better estimation of the spectrum and the unknown parameters.Using these estimated parameters we estimate f `for a low `as well.
We carried out simulations for a Laplace-Beltrami model.In practice when a set of observations is given, we have to decide how to choose the level of truncation which we denoted by L. The truncation parameter L depends on the number of observations and is connected to the number of angles where C (cos ) will be estimated.We used the Gauss-Legendre quadrature for calculating the spectrum from the covariance and vice versa.It involves L + 1 angles, actually the roots of the Legendre polynomial P L+1 (x).We chose the method of nonlinear least squares for estimation of the parameters c from the estimated spectra f `(c).Other methods like weighted least squares, MCMC and likelihood and noisy data are the subject of further investigations.

A Examples
The following are examples of homogenous and isotropic random …eld T 0 (x), x 2 R 3 restricted to the sphere S 2 .Then the covariance function C (r) of the stochastic …eld T (x) on the sphere S 2 equals the covariance function C 0 (r) of the original …eld T 0 restricted to the sphere.At the same time the power spectrum f `of the …eld T (x) is de…ned by the power spectrum of the …eld T 0 (x) through a formula (38) called Poisson formula, see [55], VII. 2.
Example 2 For an homogeneous isotropic random …eld on R 3 we have the spectral representation of the covariance function, where (d ) is some spectral measure, see [66], and where j m is the Spherical Bessel function of the …rst kind, j 0 (r) = sin r=r, see [2].If we consider two locations x 1 and x 2 on the sphere S 2 with angle 2 [0; ], then we obtain the covariance function C (cos ) = C 0 (2 sin ( =2)) on the sphere S 2 with spectrum where J `+1=2 denotes the Bessel function of the …rst kind, see [2].More generally, in case of R d , d 3, and the corresponding spectrum on the sphere is Example 3 Laplace model restricted to the sphere.The following covariance function corresponds to an homogeneous isotropic random …eld T on R d satisfying the equation in the L 2 sense, where , denotes the Laplace operator on R d , d 3, and @W is the white noise in R d .The covariance function of spherical random …eld T (x), x 2 S d 1 , restriction of the homogeneous isotropic random …eld T (x), x 2 R d , into the sphere where K is the modi…ed Bessel (Hankel) function of the second kind.Here 2 d 2 > 0, is the smoothness parameter which controls the continuity, and c controls the regularity [20], [54].Note K (r) ( ) (r=2) =2 if r !0. The correlation function on the sphere S d 1 is The corresponding spectrum on S 2 is cos ) P `(cos ) sin d : Note that C (cos ) belongs to the Matérn Class of Covariance Functions [42], [51].Also in [32] one can …nd a proof of the form of the covariance function of Matern class from the fractional Helmholtz equation based on the theory of generalized random …elds [21].In particular for = 1, d = 3, we have with spectrum and the spectrum is Z 0 e 2jcj sin( =2) P `(cos ) sin d : The preceding example treated an homogeneous isotropic random …eld on R d and then specialized it to the sphere.Another possible construction of covariance functions is based on the following.The covariance function C 2 (x 1 x 2 ) in ( 3) is strictly positive de…nite if all f `are 0, and only …nitely many of them are zero ( [53], [52] ).Therefore if the series (3) is …nite and only …nitely many f `= 0, then one can construct a Gaussian …eld with covariance function C (cos ) which is nonnegative de…nite, see Remark 1 also.In the case where …nitely many f `> 0, then C (cos ) is still nonnegative de…nite but will not be necessarily strictly positive.if 0 < z < 1, (see [66]).Since it is positive de…nite it can be considered as a covariance function on S 2 , in this case the spectrum f `is not given by some explicit formula.Some more examples of this type can be constructed applying formulae of series of Legendre polynomials with positive coe¢cients.
There is an other application of series of Legendre polynomials in probability theory, namely in directional statistics.A probability density of a rotational symmetric distribution on the sphere has a series expansion in terms of Legendre polynomials, see [45], [38].Now if the coe¢cients (actually the characteristic function of the distribution) of the series expansion are positive then the same function may also serve as a covariance function.The following example is one of the basic density on the sphere.

Example 5
The Fisher probability density function f on the sphere (see [45], [38], [10], [63]) The variance 2 in these examples corresponds usually to some additional noise …elds on top of the homogeneous isotropic …eld considered here.Since it is a multiplicative constant, it will not in ‡uence our results and therefore we will set 2 = 1 from now on.

B White noise analysis on the sphere
Here we outline a way to make precise the derivation of the spectrum given in Example 1.Details will appear in a forthcoming paper.Recall the notations x 2 S 2 , x (#; ') = (sin # cos '; sin # sin '; cos #), (dx) = sin #d#d'.

Figure 1 : 2 Z
Figure 1: The sphere S 2 with the circle C (x; ) The empirical covariance b C (cos ) for an angular distance will be given in two steps.First we …x a location x and superpose T (x) T (y (x; )) d =2 over all y (x; ) on the circle C (x; ) by varying , then secondly, we integrate over all x on the sphere S 2 , yielding b C (cos ) =

Figure 3 :
Figure 3: Listed from the middle (cos = 0), from top to bottom: theoretical (33) with c = 2; upper con…dence curve, estimated (31), lower con…dence curve.The vertical dots are the results of the individual simulations from 1 to 100.

Figure 4 :
Figure 4: Estimation of the correlation function when f `= 0, `= 0; 1; 2; 3. The estimated (31) is closed to the theoretical.The vertical dots are the results of the individual simulations from 1 to 100.