Sample variance in weak lensing: how many simulations are required?

0.4 MB·

by  Andrea Petri

Sample variance in weak lensing: how many simulations are required?
Andrea Petri,1,2,∗ Zolt´an Haiman,3 and Morgan May2
1Department of Physics, Columbia University, New York, NY 10027, USA
2Physics Department, Brookhaven National Laboratory, Upton, NY 11973, USA
3Department of Astronomy, Columbia University, New York, NY 10027, USA
(Dated: March 9, 2016)
Constrainingcosmologyusingweakgravitationallensingconsistsofcomparingameasuredfeature
vectorofdimensionNb withitssimulatedcounterpart. AnaccurateestimateoftheNb×Nb feature
covariance matrix C is essential to obtain accurate parameter confidence intervals. When C is
6 measured from a set of simulations, an important question is how large this set should be. To
1 answer this question, we construct different ensembles of Nr realizations of the shear field, using
0 a common randomization procedure that recycles the outputs from a smaller number Ns ≤ Nr of
2 independentray-tracingN–bodysimulations. Westudyparameterconfidenceintervalsasafunction
r of(Ns,Nr)intherange1≤Ns ≤200and1≤Nr .105. Previouswork[1]hasshownthatGaussian
a noise in the feature vectors (from which the covariance is estimated) lead, at quadratic order, to
M
an O(1/Nr) degradation of the parameter confidence intervals. Using a variety of lensing features
measured in our simulations, including shear-shear power spectra and peak counts, we show that

8 cubic and quartic covariance fluctuations lead to additional O(1/N2) error degradation that is not
r
negligible when Nr is only a factor of few larger than Nb. We study the large Nr limit, and find
O] thatasingle,240Mpc/hsized5123-particleN–bodysimulation(Ns=1)canberepeatedlyrecycled
to produce as many as Nr = few ×104 shear maps whose power spectra and high-significance
C peakcountscan betreated as statistically independent. Asaresult, asmall numberof simulations
h. (Ns=1 or 2) is sufficient to forecast parameter confidenceintervals at percent accuracy.
p
PACSnumbers: 98.80.-k,95.36.+x,95.30.Sf,98.62.Sb
-
o Keywords: WeakGravitational Lensing—Simulations—Methods: analytical,numerical,statistical
r
t
s
a I. INTRODUCTION This workstudies these issues further, focusing on the
[ number ofindependent N–body simulationsrequiredfor

2 Weakgravitationallensing(WL)isapromisingcosmo- an accurate estimate of the parameter constraints. Ray-
v logical probe for constraining the dark energy equation tracingsimulationsthatresolvethe cosmicstructuresre-
2 of state w, and has been considered by a range of past sponsible for lensing on arcminute scales are limited to
9 (CFHTLens [2,3], COSMOS[4]), ongoing(DES[5])and physicalsizes of hundreds ofMpc, andthus covera solid
7 future (LSST [6], Euclid [7], WFIRST [8]) experiments. angle of only O(10 deg2). As a result, many simulations
6 In an era where cosmology is data driven, accurate nu- arerequiredtotileasignificantfractionofthesky,andto
0
merical simulations of shear fields are becoming impor- make predictions for large “all-sky” surveys, such as the
.
1 tant for several reasons, including assessing baryonic ef- ones by DES, LSST, Euclid, WFIRST. In practice, this
0 fects[9–14],theutilityofnon–Gaussianstatistics[15–23] has led to the wide-spread use of “pseudo-independent”
6
and various systematic effects [24–27]. realizations, i.e. a procedure in which one randomizes
1
: Afundamentalissuewithpredictionsfromsimulations and re-cycles the output of a single 3D simulation mul-
v is that the finite number of simulations naturally intro- tiple times. In light of the forthcoming large surveys, it
Xi ducesfluctuationsintheforecasts,duetoinevitablesam- is imperative to assess the statistical validity of this ap-
plevariance[28]. Ingeneral,quantitiessuchasthe mean proach, and to ask how many times a single simulation
r
a orthevarianceofanyfeature(e.g. theshearpowerspec- can be fairly recycled. In this paper, we address these
trumatamultipoleℓ),measuredfromafinitesetofsim- questions with ensembles of up to Nr = 105 random re-
ulations, will fluctuate, and can alsosuffer a bias. While alizations, extracted from up to Ns = 200 independent
biasesintheestimatesofboththemeanandthevariance ray-tracing N–body simulations. We focus in particular
havebeen studiedextensively,the impactof fluctuations on the parameter w, and on two different statistics: the
in the variance has receivedless attention. These fluctu- (convergence)powerspectrumandthe numbercountsof
ations have been shown to have non-negligible effects on peaks.
estimatesoffeaturescovariancesandhenceonparameter This paper is organized as follows. In § II, we sum-
constraints. In particular, in the limit of Gaussian fluc- marize the shear simulation methods we utilized, and
tuations, the parameter confidence limits are degraged describe the formalism we adopted to forecast cosmo-
by a factor 1+O(1/N ) [1, 29]. logical parameter constraints. We then vary the number
r
ofsimulationsandthenumberofpseudo-independentre-
alizations, and present our main findings in § III. These
resultsarediscussedfurtherin§IV.Weofferourconclu-
∗ [email protected] sions, and suggest follow-up future work in § V.
2
II. METHODS • Repeattheaboveprocedureforeachlens-planered-
shift z .
l
A. Ray-tracing simulations of the convergence field
This randomization procedure allows us to produce an
(almost)arbitrarynumberN ofshearrealizationsγγγ (θθθ).
r r
In this section, we describe how we constructed our
However, these realizations are not guaranteed to be in-
shear field ensembles. Background galaxies at redshift
dependent,ifN isnotlargeenough. Usingthesetof200
s
z are lensed by large scale structures between z =
s independent N–body simulations, we construct different
0 and z . The shape distortions due to the cosmic
s ensembleswithdifferentchoicesofN ∈[1,200]. Eachof
shear γγγ can be computed in terms of the dark mat- s
these ensembles consists of the same number N =1000
ter gravitational potential Φ(x,z). Because the evolu- r
ofshearrealizations. We alsobuild anadditionalensem-
tion of Φ with redshift is non–linear, it needs to be
blewithN =1andN =105realizations. Foreachreal-
s r
computed with numerical simulations. We make use of
izationofeachensemble,we reconstructthe convergence
the public code Gadget2 [30], with which we run a se-
κ (θθθ) from the trace of the light-ray deflection Jacobian
r
quence of 200 independent dark–matter–only N–body
matrix,measuredfromthe differenceindeflectionangles
simulations that track the evolution of the density fluc-
between nearby light-rays [31, 32, 34].
tuations. We assume a standard ΛCDM background
We measure the κ angularpower spectrum Pκκ(ℓ) de-
universe with the parameters (Ω ,Ω ,h,w,σ ,n ) = r
m Λ 8 s fined as
(0.26,0.74,0.72,−1,0.8,0.96). We fix the comoving size
of the simulation box to 240Mpc/h, and use 5123 par- hκ˜ (ℓℓℓ)κ˜ (ℓℓℓ′)i=(2π)2δ (ℓℓℓ+ℓℓℓ′)Pκκ(ℓ) (1)
ticles, corresponding to a dark matter particle mass of r r D r
≈1010M .
⊙ As an additional summary statistic, we consider the
We assume a uniform galaxy distribution at a con- counts of local κ maxima of a certain height κ , n (κ )
0 r 0
stant redshift zs = 2 (at which the simulation box has (hereafter peak counts), with varying κ0 chosen between
an angular size of θbox = 3.5◦) and we discretize the the minimum and maximum values measured from the
mass distribution between zs and the observer at z = 0 maps (κmin,κmax) = (−0.06,0.45). Different choices of
with a sequence of 46 two dimensional lenses of thick- κ binning used in this work are outlined in Table II.
0
ness 80Mpc/h. The surface density on each lens plane ThefactthattheensembleofN realizationsisnotcom-
r
iscomputedbyprojectingthethree–dimensionaldensity pletely independent if N is not large enough can have
s
measured from Gadget2 snapshots. We then apply the an effect on the covariance estimators of both Pκκ and
multi–lens–plane algorithm (see [31, 32] for example) to n(κ ).
0
trace the deflections of n2ray = 20482 light rays arranged Tomeasurethecosmologicaldependenceoftheκpeak
onasquaregridoftotalsizeθbox,fromz =0tozs. This counts,weperformedasetofadditionalray–tracingsim-
correspondstoapixelangularresolutionof0.1′. Ourim- ulations with different combinations of the cosmological
plementation of this algorithm is part of the LensTools parameter triplet (Ω ,w,σ ). A summary of the com-
m 8
computing package we have been developing [33], and plete set of shear ensembles used in this work is listed in
havereleasedunder the MITlicense. Manydifferent real- Table I.
izations r of the same shear fieldγγγ (θθθ) can be generated
r
by picking different lens planes that lie between the ob-
serverand zs. The randomizationprocedurewe adopt is B. Cosmological parameter inference
the following (see [34] for reference):
Let dˆ be a single estimate for a feature of dimension
• For each lens-plane redshift z , select the snapshot
l N , d(p) be the true value of this feature ata point p in
b
at z from the i–th N–body simulation, where i is
l parameter space (which has a dimension N ) and C be
p
a random integer i∈[1,N ].
s the N ×N feature covariance matrix. For the purpose
b b
of this work p is the triplet (Ω ,w,σ ) and d is one of
m 8
• Chooserandomlybetween the three orthogonaldi-
the features – either a power spectrum or a peak count
rections n ,n ,n : the lens plane will be perpen-
x y z histogram–inTableII.Althoughexistingemulatorscan
dicular to this direction.
be used, in principle, to compute both d(p) and C, the
latterismoredifficult,andtypicallyonlythemean,d(p),
• Choose the position of the plane along the snap-
has been computed to date (refs. [35, 36], but see an
shot: because the lens thickness is 1/3 the size of
exception by ref. [37]). Estimating C from simulations
thebox,wecancutthreedifferentslicesofthesim-
involves generating a series of mock realizations dˆ with
ulation box for each orientation n ,n ,n . This r
x y z r =1...N and computing the sample covariance Cˆ,
givesatotalof9choicesforgeneratingalensplane r
out of a single N–body snapshot.
1 Nr
• Perform a periodic random shift of the lens plane d¯ = dˆ , (2)
r
along its two directions. Nr
r=1
X
3
Becausepˆ isestimatedusingasinglenoisydatainstance
dˆ , its estimate willbe scatteredaroundthe true value
obs
1 Nr hpˆi . Inthe followingweusethehi notationforexpec-
Cˆ = (dˆ −d¯)(dˆ −d¯)T. (3) O O
N −1 r r tationvaluestakenwithrespecttoobservations,whilewe
r
Xr=1 keepthenotationhiforexpectationvaluestakenwithre-
Assuming a normal feature likelihood, together with a spect to the simulations. Defining the precision matrix
flatpriorontheparameterspace,theparameterposterior Ψˆ = Cˆ−1, we can express the estimator of the observa-
distribution L(p|dˆ ) given an observed instance of dˆ, tional scatter in pˆ:
obs
which we call dˆ , follows from Bayes’ theorem,
obs
−2logL(p|dˆobs)=[dˆobs−d(p)]TCˆ−1[dˆobs−d(p)]. (4) Σˆ =Fˆ−1d′TΨˆh(dˆ −d )(dˆ −d )Ti Ψˆd′Fˆ−1, (8)
p 0 obs 0 obs 0 O 0
For the sake of simplicity, we approximate the posterior
as a Gaussian around its maximum. This corresponds
to Taylor-expanding the simulated feature to first order
around a point p0 (ideally the maximum of Eq. 4): Fˆ =d′TΨˆd′. (9)
0 0
d(p)≈d +d′(p−p ). (5)
0 0 0
Here we introduced the familiar Fisher matrix estimator
We chose p to be the triplet (Ω ,w,σ ) = Fˆ = d′TΨˆd′ and, for simplicity, we assumed hdˆ i =
0 m 8 0 0 obs O
(0.26,−1,0.8). Tomeasurethederivativesofthefeatures d , so that h(dˆ −d )(dˆ −d )Ti = C. When we
0 obs 0 obs 0 O
d′0 withrespecttothecosmologicalparameters,wemake perform an observation dˆobs, the parameter estimate pˆ
use of the public code Nicaea [36] for the power spec- is a random draw from a probability distribution with
trum,andweuseanindependentsimulationset(contain- variance Σˆ , which inherits noise from the simulations.
p
ing simulations with a variety of different combinations The noise in the covariance estimator (Eq. 3) and in its
of (Ωm,w,σ8), see Table I) for the peak counts. inverseΨˆ propagateallthe wayto the posterior(Eq. 4),
Wecanbuildtheestimatorfortheposteriormaximum
the parameter estimate (Eq. 6) and its variance (Eq. 8).
pˆ, given the observation dˆ , as follows:
obs Following refs. [1, 29, 38] we compute the expectation
pˆ =p +Tˆ(dˆ −d ), (6) value of Eq. (8) over simulations, hΣˆpi, up to O(1/Nr2)
0 obs 0
by expanding Eq. (8) to quartic order in the statistical
fluctuations of Ψˆ. Denoting the true parameter covari-
ance, i.e. the usual inverse Fisher matrix, as Σ = F−1,
p
Tˆ =(d′TCˆ−1d′)−1d′TCˆ−1. (7) we find the result [39]:
0 0 0
N −N (N −N )(N −N +2) 1
hΣˆ i=Σ 1+ b p + b p b p +O . (10)
p p N N2 N3
(cid:20) r r (cid:21) (cid:18) r (cid:19)
Although we truncated the expansionto second order in
1/N ,anexactexpressionforhΣˆ ihasbeenproposedby
r p
[40] N −N
hΣˆ i=Σ 1+ b p . (12)
p p
hΣˆ i =Σ Nr−2 (11) (cid:18) Nr (cid:19)
p empirical p
N −N +N −2
(cid:18) r b p (cid:19)
This is the result obtained by ref. [1].
This empirical expressionreduces to equation (10) when
expanded at order O(1/N2) but, to our knowledge, no
r
first principles proof of its correctness exists. Next, we
2. Usually the true data covariance is unknown, and
restrictourselvestothelargeNr limit,andwefurtherin- it is tempting to plug in its estimator Cˆ, measured
vestigatethebehavioroftheO(1/Nr)term. Weconsider fromthesamesimulationsetweusetocomputeΨˆ.
three cases [41]:
Thisapproachhasbeenusedbeforeintheliterature
1. If the true data covariance C is known, the esti- (e.g.[17,24]). Ifthisis donewithoutcorrectingfor
mator in eq. (8) is biased, and the dominant con- the bias in Ψˆ (see Ref. [29] and eq. 20 below), the
tribution of the bias comes from the second order parameter variance will have a contribution from
fluctuationsinΨˆ. Oncetheexpectationvaluesover both the second and first-order fluctuations in Ψˆ,
simulations are taken, the bias sums up to which now have a nonzero expectation value. In
4
Ωm w σ8 (Ns,Nr) Numberof κ ensembles
0.26 −1 0.8 (1 to 200,1024) 16
0.26 −1 0.8 (1,128000) 1
0.29 −1 0.8 (1,1024) 1
0.26 −0.8 0.8 (1,1024) 1
0.26 −1 0.6 (1,1024) 1
TABLE I. Summary of the shear ensembles used in this work. Ns and Nr refer to the number of independent N–body
simulations, and thenumberof pseudo-independentrealizations created from these simulations, respectively.
this case the bias sums up to III. RESULTS
N −N
hΣˆ i=Σ 1− b p . (13)
p p Inthissectionwepresentthemainresultsofthiswork.
N
(cid:18) r (cid:19) We show the qualitative behavior of a variety of feature
3. Ifwerepeatthesameexerciseasabove,butwecor- dˆ probability distribution functions (PDFs) in ensem-
r
rect for the bias in the precision matrix estimator, blesbuiltwithdifferentN andN . InFigure1,weshow
s r
we are left with the PDF of the power spectrum at four selected multi-
poles, spanning the linear (ℓ = 115) to the nonlinear
(ℓ = 5283) regime. In Figure 2, we shows the ensemble
1+N
hΣˆ i=Σ 1+ p . (14) mean for these power spectra, as well as for peak counts
p p
N
(cid:18) r (cid:19) ofthree different κ0 heights (correspondingto ≈2−13σ
The error degradation in each parameter p, at leading peaks), as a function of Ns. In Figure 3, we show the
order, scales as D/N , where D = N −N , N −N , variance of the power spectrum at each multipole, as a
r b p p b
and 1 + Np for cases 1, 2, and 3, respectively. Note function of Ns, in units of the variance expected if the
that in the last case, which is most relevant when fitting convergence κ was a Gaussian random field
actual data, the estimated degradation turns out to be
(Pκκ)2
too optimistic: the parameter estimate pˆ has a variance Var(Pκκ)= ℓ (16)
whose noise growslinearly with Nb (eq. 10), whereasthe ℓ Neff(ℓ).
degradation estimated via eq. (14) is constant with N .
b HereN (l)isthe numberofindependentmodesusedto
Thiscanleadtounderestimationoferrorbars,whichcan eff
estimate the power spectrum at ℓ.
be mistakenly interpreted as a parameter bias. We test
In practice, we measure Pκκ on the Fourier transform
scaling relations of the form ℓ
ofthepixelizedsimulatedmapκ (θθθ),usingtheFFTalgo-
r
hσˆ2i=σ2 (N ) 1+ D (15) rithm,andsomecaremustbetakentocountthenumber
p p,∞ s N of modes N (ℓ) correctly. Each pixel (i ,j ) in Fourier
(cid:18) r(cid:19) eff x y
space corresponds to a mode (ℓ ,ℓ ) = 2π(i ,i )/θ ,
against our simulations, in the limits of both high and x y x y box
low N . We indicate the diagonal elements of Σˆ as with ix =−nray/2,...,nray/2and iy =0,...,nray/2. Here
r p
n = 2048 is the linear number of pixels on the ray–
σˆ2 = diag(Σˆ ) and we indicate by σ2 the expectation ray
p p p,∞ tracedconvergencemaps. Wecountthenumberofpixels
valueofthe varianceofeachparameterinthe limitofan
N(ℓ)thatfallinsideamultipolebin(ℓ ,ℓ ). Becausethe
infinite number of realizations N → ∞. We call D the 1 2
r κ field is real, the modes (±ℓ ,0) are not independent.
effective dimensionality of the feature space (which, as x
Ifwe letN(ℓ,ℓ =0)be the number ofnon–independent
seen before, can be negative in some pathologicalcases). y
modes, the effective number of independent modes for
We compute the expectation values of σˆ2 (eqs. 8 and
p the variance is given by
15) by averaging over 100 random resamplings of our
shear ensembles. For the true feature covariance matrix
h(dˆ −d0)(dˆ −d0)Ti = C we use the estimated covari- N2(ℓ)
N (ℓ)= . (17)
ance from a grand ensemble built with the union of all eff N(ℓ)+N(ℓ,ℓ =0)
y
the ensembles with different N .
s
Thetrueparametervarianceσ2 (N )inprinciplecan This correctionis important at low ℓ, where pixelization
p,∞ s
depend on the number of independent N–body simula- effects are non-negligible; N (ℓ≫2π/θ )≈N(ℓ).
eff box
tions N , which appears in the randomizationprocedure In Figures 4 and 5, we show the dependence of the
s
described in § IIA above. This is because if N is not confidence range hσˆ2i on N , derived from the features
s w r
largeenough,thedifferentshearrealizationscannotbeall used in this work (see Table II for a comprehensive list).
independent,andhencethetruevarianceσ2 (N →∞) Figure 4 shows the behavior in the limit of a large num-
p,∞ s
cannot be recovered for low N even if N is arbitrarily ber N ≫ 500 of realizations, and compares it with the
s r r
large. In the next section, we present our main findings. scaling of the form in equation (15). Figure 5 shows the
5
l=115 l=344
6 1e7 5 1e8
Ns=1
5 Ns=2
Ns=5 4
Ns=50
4
κκP)l3 NNss==110,N0r=128000 κκP)l3
( (
L L2
2
1
1
0 0
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 0.0 0.2 0.4 0.6 0.8 1.0 1.2
Pκκ 1e−7 Pκκ 1e−8
l l
l=1349 l=5283
8 1e9 1.2 1e11
7
1.0
6
0.8
5
) )
κ κ
κPl4 κPl0.6
( (
L L
3
0.4
2
0.2
1
0 0.0
0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 0.6 0.7 0.8 0.9 1.0 1.1
Pκκ 1e−9 Pκκ 1e−10
l l
FIG.1. PDFoftheκpowerspectrumL(Pκκ)atfourselectedmultipolesℓ=115,344,1349,5283,fordifferentshearensembles
l
constructedfromonNs =1(black),2(blue),5(green),50(red),and100(purple)independentN–bodysimulations. Eachcurve
isbasedonNr =1024realizations. ThedashedblackcurvescorrespondtoensemblesgeneratedwithNs =1andNr =128000.
For Ns ≥ 2, the distributions appear similar to the eye; this similarity is confirmed by the comparisons in Figures 2 and 3
below.
Feature Specifications Nb Symbol Color
Power Spectrum,log binning ℓ∈[100,800] 8 × black
Power Spectrum,log binning ℓ∈[1000,6000] 7 (cid:4) black
Power Spectrum,log binning ℓ∈[100,6000] 15 • red
Power Spectrum,linear binning ℓ∈[100,2000] 15 + red
Power Spectrum,linear binning ℓ∈[2500,4500] 15 × red
Power Spectrum,linear binning ℓ∈[100,4500] 30 • green
Power Spectrum,linear binning ℓ∈[100,6000] 39 • blue
Low peaks κ0 ∈[−0.06,0.09] 15 + red
Intermediate peaks κ0 ∈[0.1,0.27] 15 ⋆ red
High peaks κ0 ∈[0.28,0.45] 15 ⋄ red
Low+Intermediate peaks κ0 ∈[−0.06,0.27] 30 × green
Intermediate+High peaks κ0 ∈[0.1,0.45] 30 (cid:4) green
Allpeaks κ0 ∈[−0.06,0.45] 45 (cid:4) magenta
TABLE II. Catalog of feature types used in this work, along with the chosen number of bands Nb and the plot legends for
Figure 4.
6
0.6 l=115 1000 500 300 200 150 Nr 100 90 70
l=1027
(N)d(200)]/C(200)siii−(cid:0)−−00000.....40422 lκκκ=000===520008...3012578 22ˆ/σσw,w∞11222(cid:0)(cid:1).....68024 ℓℓκκκℓℓℓκκℓℓκ000000∈∈∈∈∈∈∈∈∈∈∈∈∈[[[[[[[1111212[[[[[[0000505000−−−0000000...000,,,,2110008264,,,8,,...600045600,00000005000..66624000]000.,,,,4750]]0]0000N,,,5]]]]]...NNN,,,,,]b024NNNNN,bbb=975Nbbbbb===]]],,,b8=====NNN(113=l13755103bbbo150(((5(9===gllll5((oioi)llnn134giignn))505))))
di
[
1.4
−0.6
1.2
−0.8
0 50 100 150 200
N
s 1.0
0.000 0.002 0.004 0.006 0.008 0.010 0.012 0.014
1/Nr
FIG.2. Themeanvalued¯ofvariousfeatures,measuredfrom
ensembles created from different numbers Ns of simulations. FIG.4. Expectationvalueofthevarianceofwcomputedfrom
For each case, the difference compared to the mean in the
equation(8),shownasafunctionof1/Nr. Thedifferentsym-
Ns = 200 ensemble is shown, in units of the statistical error bols and colors correspond to the features listed in Table II.
measuredintheNs =200ensemble. Thecoloredcurvesrefer The dashed and this solid curves show the analytic predic-
to shear–shear power spectra measured at ℓ = 115 (black), tions from equation (10) at orders O(1/Nr) and O(1/Nr2),
1027 (cyan), and 5283 (green), and peak counts with heights
respectively. The thick solid curves show the empirical pre-
κ0 =0.05 (red), 0.17 (purple), and 0.28 (orange). The κ bin dictions from equation (11). The asymptotic variance σw2,∞
widthforthepeakcountshasbeenfixedto∆κ=0.011. The has been computed from a linear regression of hσˆw2i vs 1/Nr
dashedblacklineshowsalevelof0.1σ accuracyforreference.
for Nr > 500. The figure clearly shows that terms beyond
ForNs ≥2,themeansarestatisticallyindistinguishable(even O(1/Nr) need to be considered, unless Nr ≫Nb.
at ∼0.1σ) from those in theensemble with Ns =200.
9
Ns=1 103
8 Ns=2
Ns=5 102
7 Ns=10
κκ,2Pl6 NNss==51000 110001
/
)Nl×5 2σw,∞10-1
κκPl4 −1(cid:0)0-2
Var(3 2ˆσw1(cid:1)0-3
2 10-4
1 10-5 ℓℓ∈∈[[110000,,6205000],]N,Nbb==41(l5o(glo)g)
10-6 κ0∈[0.44,0.48],Nb=4
00 1000 2000 3000 4000 5000 6000 κ0(θG=1′)>0.15,Nb=20
l 10-7
102 103 104 105
Nr
FIG. 3. Variance of the κ power spectrum as a function of
the multipole ℓ, in units of the expected Gaussian variance FIG.5. Biasinthevarianceofw,hσˆw2i−σw2,∞,asafunctionof
from equation (16). The variance is measured from different thenumberofrealizationsNr usedtoestimatethecovariance
shear ensembles based on Ns =1 (black), 2 (blue), 5 (green), (eq.3). Thefigureshowsboththetrendmeasuredinthesim-
10 (red), 50 (purple), or 100 (orange) N–body simulations. ulations(solidlines)andtheirscalingexpectedfromEq.(15)
Non–Gaussianities of the underlying structures increase the with D = Nb−Np (dashed line). The asymptotic variance
variance on small scales, but no clear trend with Ns can be σw2,∞hasbeenestimatedtobethevaluehσˆw2i(Nr =105). Dif-
identied on any scale. ferentfeaturesareconsidered: powerspectrawithlogarithmi-
cally spaced ℓ∈[100,6000] (black), ℓ∈[100,250] (red), peak
counts in the unsmoothed maps with height κ0 ∈[0.44,0.48]
(green)andpeakcountsinthesmoothedmaps(withaGaus-
′
large N trends of the w constraint. Figure 4 illustrates sian kernel of size θG =1) with height κ0 >0.15 (blue). No
the behravior at relatively low Nr, and compares hσˆw2i dtoevNiarti≈onfsewfr×om10t4h,eexexceppetctfeodrt1h/eNlrarbgeeh-sacvailoerpaorweeorbssperevcetrduump,
measureddirectlyfromthesimulationswiththeanalytic in which case thedeviations occur much earlier (Nr ≈103).
expectationsfromequation(10). Finally,inFigure6,we
show how the w confidence limit changes with N .
s
7
of the Gaussian expectation. We find that, even with
Power spectrum log binning
1.06 Power spectrum linear binning Ns =1,weareabletorecovertheknownresultthatnon–
Peak counts
Gaussianstructuresincreasethevariancesignificantlyon
smallscales(see[34,43]forreference). Ourresultsarein
mean1.04 fact in excellent quantitative agreementwith [34], which
2σ,∞ used Ns = 400 independent N–body simulations. This
/
) result is highly encouraging, suggesting that individual
Ns1.02
( N–body runs can be recycled repeatedly. However, it is
2 ∞
σ not sufficient by itself to conclude that N does not im-
s
1.00 pactthe parameterinferences,since these dependonthe
cross band covariances.
Figure4investigatesthe parametererrors. This figure
0.98
0 50 100 150 200 shows that error degradation estimates truncated at or-
N
s derO(1/N )aretoooptimisticwhenthenumberofsimu-
r
lationsN usedtomeasurethecovarianceisonlyafactor
r
FfrIoGm.t6h.eTinhteevrcaerpiatnocfetohfewfitinσw2thvesli1m/Nitro)f,Nvarry→in∞gth(meneausmubreedr oInfftehweslearcgaesrestheaffnectthsecdoimmienngsitohneonfetxhte–tfeoa–tleuardeisnpgaocerdNerbs.
oafncsiem(uelqa.ti3o)n.sWNessuhsoedwinthtehdeeepnesnemdebnlceetoofesσtw2im,∞a(tNesth)eincouvnairtis- O(1/Nr2) become non–negligible on constraint degrada-
of the mean over the union of 16 ensembles with different tion. In particular,we find that alreadyfor Nb =30 and
Ns and Nr = 1024, for the power spectrum logarithmically Nr ∼ 100, the error degradation estimates to the next
binned (black, Nb =15,ℓ∈[100,6000]), the power spectrum leading order, O(1/Nr2), remain too optimistic. Accu-
linearly binned (red,Nb = 39,ℓ ∈ [100,6000]) and the peak rate analytic estimates in this regime would require at
counts (green,Nb =45,κ0 ∈[−0.06,0.45]). No trend with Ns least terms of order O(1/Nr3), which come from higher–
is seen for Ns ≥2, and the differences are only of order 1%. than–quartic Ψˆ fluctations.
InFigure 5, we examine how the degradationin the w
constraintdepends onthe number ofsimulationsusedto
IV. DISCUSSION
estimatethecovariance,inthelimitoflargeN . Wefind
r
excellent agreement with the expected scaling (eq. 15)
In this section we discuss our main findings and their up to N ∼few×104 when using the κ power spectrum
r
implications. Figure 1 shows that, although different in the multipole range ℓ ∈ [100,6000]. The same be-
choices of Ns do not affect the power spectrum PDF on haviorisobservedwhenconsideringthehigh-significance
large scales (top two panels), there are some qualitative peak counts (> 10σ for unsmoothed maps and > 5σ for
differences on smaller scales (bottom two panels). On 1′ smoothed maps). As the figure shows, around these
these smaller scales, shear ensembles built from Ns = 1 values of Nr the hσˆw2i−σw2,∞ curve becomes noisy and
donotproducethesamestatisticalbehaviorasensembles reaches negative values. This is a clear indication that
built with larger Ns. In particular, looking at the black the 1/Nr behavior is broken and a plateau in hσˆw2i is
curves, we see that the Ns = 1 ensembles exibit large reached. The negative values in the plot are a conse-
shifts with respect to the other PDFs to lower power, quenceofthenoiseintheestimationofthisplateauvalue
including the locationsofthe peaksofthe PDFs. We at- (or equivalently in the estimated value of σ2 ).
w,∞
tribute these offsets to large (random) statistical errors. We conclude that a single N–body simulation is suf-
Interestingly, we need as few as Ns = 2 simulations ficient to construct an ensemble of up to a few×104
to recoverthe rightPDF for the small–scalepower spec- mutually independent convergence power spectra. For
trum. Figure2showsthatmultipleindependentN–body N ≫104,theshearrealizationscannolongerbeconsid-
r
simulations Ns ≥ 2 are indeed necessary for measuring ered independent. We emphasize that the precise value
the means of feature ensembles to an accuracy corre- of this N will depend on the size of the simulation box
r
sponding to 10% of the statistical error. The number (which,inourcase,is(240Mpc/h)3,with5123 particles)
of required simulations Ns depends on the feature type andalsoontherangeofmultipolesℓusedtoconstrainthe
and ranges from a few (Ns =1 or 2) for the power spec- parameters. Figure 5 shows that when we infer w only
trum at low multipoles (ℓ . 500) to Ns ≈ 30−50 for from large–scale modes, ℓ . 250, the plateau is reached
the power spectrum at larger multipoles (ℓ & 1000) or at least an order of magnitude earlier in the number of
peak counts above a high threshold (κ0 ≈ 0.3). On the realizations. In other words, the number of independent
otherhand,relaxingthe requiredaccuracyto50%ofthe power spectra we can generate decreases as we increase
statistical error, we find Ns = 2 to be always sufficient. the spatial scales of interest. This is due to the fact
Aspointedoutby[42],theboxsizeusedfortheN–body that, because of the finite box size, the number of in-
simulations can also play an important role in the accu- dependent lens plane shifts (as described in § IIA) de-
racy of the power spectrum ensemble means. creases as the mode size approaches the size of the box.
Figure 3 shows the variance of convergence power Similarly, one may expect that the independence in the
spectrum computed from different ensembles, in units statistics of high-amplitude peaks, which are predomi-
8
nantlyproducedbysinglemassivehalos,maybecompro- • As few as one or two independent N–body simula-
mised by these halos being present repeatedly, in many tions are sufficient to forecast w error bars to 1%
of the pseudo-independent realizations. However,Figure accuracy,providedthat a sufficiently largenumber
5 shows that this is not the case: the peak count statis- N of realizations are used to measure feature co-
r
tics are shown at κ thresholds corresponding to massive variances. In particular, provided that biases in
(≈1015M )halos,yetthereisnoevidencethattheinde- the inverse covariance are corrected, percent–level
⊙
pendence of the maps breaks down until N =few×104. forecasts require N >100(N −N ) realizations.
r r ∼ b p
Apparently, randomly projected structures, which vary • Depending on the feature type used to constrain
fromrealization-to-realization,contributesignificantlyto cosmology,alargernumberofN–bodysimulations
the statistics of these high peaks. might be needed to measure accurate ensemble
Figure 6 shows how the “true” w constraint (in the means to an accuracy corresponding 10% of the
limit N → ∞; or equivalently the w constraint with statistical error. If this accuracyrequirementis re-
r
the knownN dependence factoredout),depends onN . laxed to 50% of the statistical error, we find that
r s
We find that, in the range N ∈ [1,200] the inferred w– as low as N = 2 simulations are sufficient for the
s s
variance σ2 fluctuates stochastically only by 1%, and feature types we consider in this work.
w,∞
does not show any trend with N .
s
Finally, we found that when we estimate the data co-
varianceCfromthesamesimulationsetusedtomeasure Future extensions of this work should involve extend-
Ψˆ,theeffectivedimensionalityD decreaseswithincreas- ing our analysis to a larger set of cosmological param-
ing N in the case where the Ψˆ bias is not corrected eters, and to more general feature spaces, such as the
b
(eq. 13). This N -dependence disappears when the bias ones that characterize non–Gaussian statistics (e.g. in-
b
is corrected(eq. 14). This fact that shouldbe taken into cluding higher moments of the κ field, Minkowski Func-
consideration when forecasting parameter errors purely tionals, and higher-order κ correlators). While our re-
from simulations, as the errors will otherwise be under- sults are highly encouraging, and suggest that a single
estimated. A similar conclusion was reached by [44] (al- N–body simulation can be recycled repeatedly, to pro-
though their paper did not address the impact of using duce asmany as104 independent shearpowerspectraor
the same simulation set for C and Ψˆ). peak count histograms. In order to scale our results to
largefuturesurveys,suchasLSST,itwillbenecessaryto
determine if our findings hold when challenged by larger
V. CONCLUSIONS and higher-resolution N–body simulations [45].
Inthiswork,wehaveexaminedtheeffectofforecasting
cosmological constraints based on shear ensembles gen-
eratedfromafinite numberofN–bodysimulations. Our
main results can be summarized as follows: ACKNOWLEDGEMENTS
• When the feature covariance matrix is measured
We thank Lam Hui for useful discussions. The simu-
from simulations, parameter constraints are de-
lations in this work were performed at the NSF XSEDE
graded. Thisdegradationisappreciablylargerthan
facility, supported by grant number ACI-1053575, and
the O(1/N ) computed by [1] when the number of
r
at the New York Center for Computational Sciences, a
realizations N is only a factor of few larger than
r
cooperative effort between Brookhaven National Labo-
the feature vector size N .
b
ratory and Stony Brook University, supported in part
• We can recycle a single 240Mpc/h N–body sim- by the State of New York. This work was supported in
ulation to produce an ensemble of O(104) shear part by the U.S. Department of Energy under Contract
maps whose small-scale power spectra and high- Nos. DE-AC02-98CH10886 and DE-SC0012704, and by
signficiance peak counts are statistically indepen- the NSF Grant No. AST-1210877 (to Z.H.) and by the
dent. The mean feature measured from a shear ResearchOpportunitiesandApproachestoDataScience
ensemble, though, could be inaccurate if only one (ROADS)programattheInstituteforDataSciencesand
N–body simulation is used. Engineering at Columbia University (to Z.H.).
[1] S. Dodelson and M. D. Schnei- lier, P.Simon, C. Bonnett, J. Coupon, L.Fu,J. Harnois
der, Phys. Rev.D 88, 063537 (2013), D´eraps, M. J. Hudson, M. Kilbinger, K. Kuijken,
arXiv:1304.2593 [astro-ph.CO]. B. Rowe, T. Schrabback, E. Semboloni, E. van Uitert,
[2] C. Heymans, L. Van Waerbeke, L. Miller, T. Erben, S. Vafaei, and M. Velander, MNRAS427, 146 (2012),
H. Hildebrandt, H. Hoekstra, T. D. Kitching, Y. Mel- arXiv:1210.0032 [astro-ph.CO].
9
[3] L. Miller, C. Heymans, T. D. Kitching, L. van Waer- arXiv:1309.4460.
beke, T. Erben, H. Hildebrandt, H. Hoekstra, Y. Mel- [18] L. Marian, R. E. Smith, and G. M. Bernstein,
lier, B. T. P. Rowe, J. Coupon, J. P. Dietrich, ApJ 698, L33 (2009), arXiv:0811.1991.
L. Fu, J. Harnois-D´eraps, M. J. Hudson, M. Kil- [19] M. Takada and B. Jain, MNRAS337, 875 (2002),
binger, K. Kuijken, T. Schrabback, E. Semboloni, astro-ph/0205055.
S. Vafaei, and M. Velander, MNRAS429, 2858 (2013), [20] M. Takada and B. Jain, MNRAS344, 857 (2003),
arXiv:1210.8201 [astro-ph.CO]. astro-ph/0304034.
[4] A. M. Koekemoer, H. Aussel, D. Calzetti, P. Capak, [21] M. Takada and B. Jain, MNRAS348, 897 (2004),
M. Giavalisco, J.-P. Kneib, A. Leauthaud, O. Le F`evre, astro-ph/0310125.
H. J. McCracken, R. Massey, B. Mobasher, J. Rhodes, [22] J. Berg´e, A. Amara, and A. R´efr´egier,
N. Scoville, and P. L. Shopbell, ApJS 172, 196 (2007), ApJ 712, 992 (2010), arXiv:0909.0529.
astro-ph/0703095. [23] J.P.DietrichandJ.Hartlap,MNRAS 402, 1049 (2010),
[5] The Dark Energy Survey Collaboration, ArXiv Astro- arXiv:0906.3512.
physicse-prints (2005), astro-ph/0510346. [24] M. Shirasaki, N. Yoshida, and T. Hamana,
[6] LSSTDarkEnergyScienceCollaboration,ArXive-prints ApJ 774, 111 (2013), arXiv:1304.2164 [astro-ph.CO].
(2012), arXiv:1211.0310 [astro-ph.CO]. [25] D. Bard, J. M. Kratochvil, C. Chang, M. May, S. M.
[7] J. Amiaux, R. Scaramella, Y. Mellier, B. Altieri, Kahn, Y. AlSayyad, Z. Ahmad, J. Bankert, A. Con-
C. Burigana, A. Da Silva, P. Gomez, J. Hoar, R. Lau- nolly, R. R. Gibson, K. Gilmore, E. Grace, Z. Haiman,
reijs, E. Maiorano, D. Magalh˜aes Oliveira, F. Renk, M.Hannel,K.M.Huffenberger,J.G.Jernigan,L.Jones,
G. Saavedra Criado, I. Tereno, J. L. Augu`eres, S. Krughoff, S. Lorenz, S. Marshall, A. Meert, S. Na-
J. Brinchmann, M. Cropper, L. Duvet, A. Ealet, garajan, E. Peng, J. Peterson, A. P. Rasmussen,
P. Franzetti, B. Garilli, P. Gondoin, L. Guzzo, M. Shmakova, N. Sylvestre, N. Todd, and M. Young,
H. Hoekstra, R. Holmes, K. Jahnke, T. Kitching, ApJ 774, 49 (2013), arXiv:1301.0830.
M. Meneghetti, W. Percival, and S. Warren, in [26] C. Chang, S. M. Kahn, J. G. Jernigan, J. R. Peterson,
Society of Photo-Optical Instrumentation Engineers (SPIE) ConfeYr.enAcleSaSyeyriaeds,,Z. Ahmad, J. Bankert, D. Bard, A. Con-
Society of Photo-Optical Instrumentation Engineers nolly, R. R. Gibson, K. Gilmore, E. Grace, M. Hannel,
(SPIE) Conference Series, Vol. 8442 (2012) p. 84420Z, M.A.Hodge,M.J.Jee,L.Jones,S.Krughoff,S.Lorenz,
arXiv:1209.2228 [astro-ph.IM]. P. J. Marshall, S. Marshall, A. Meert, S. Nagarajan,
[8] D.Spergel,N.Gehrels,C.Baltay,D.Bennett,J.Breckin- E. Peng, A. P. Rasmussen, M. Shmakova, N. Sylvestre,
ridge, M.Donahue,A.Dressler, B.S.Gaudi, T. Greene, N. Todd, and M. Young, MNRAS 428, 2695 (2013),
O.Guyon,C.Hirata,J.Kalirai, N.J.Kasdin,B.Macin- arXiv:1206.1378.
tosh,W.Moos,S.Perlmutter,M.Postman,B.Rauscher, [27] D. Huterer, M. Takada, G. Bernstein, and B. Jain,
J.Rhodes,Y.Wang,D.Weinberg,D.Benford, M. Hud- MNRAS 366, 101 (2006), astro-ph/0506030.
son, W.-S. Jeong, Y. Mellier, W. Traub, T. Ya- [28] Samplevarianceisabroadtermthathasbeenusedinthe
mada, P. Capak, J. Colbert, D. Masters, M. Penny, literaturetodescribearangeofphenomena.Throughout
D. Savransky, D. Stern, N. Zimmerman, R. Barry, this paper, we use it to refer to the fluctuations in an
L. Bartusek, K. Carpenter, E. Cheng, D. Content, ensemble of simulations, which represent random real-
F. Dekens, R. Demers, K. Grady, C. Jackson, G. Kuan, izations of thesame initial conditions.
J. Kruk, M. Melton, B. Nemati, B. Parvin, I. Poberezh- [29] A. Taylor, B. Joachimi, and T. Kitching,
skiy, C. Peddie, J. Ruffa, J. K. Wallace, A. Whip- MNRAS 432, 1928 (2013), arXiv:1212.4359.
ple, E. Wollack, and F. Zhao, ArXiv e-prints (2015), [30] V. Springel, MNRAS 364, 1105 (2005),
arXiv:1503.03757 [astro-ph.IM]. astro-ph/0505010.
[9] X.Yang, J. M. Kratochvil, K.Huffenberger, Z. Haiman, [31] S.Hilbert,J.Hartlap,S.D.M.White, andP.Schneider,
and M. May, Phys. Rev.D 87, 023511 (2013), A&A 499, 31 (2009), arXiv:0809.5035.
arXiv:1210.0608. [32] B. Jain, U. Seljak, and S. White, ApJ530, 547 (2000),
[10] E. Semboloni, H. Hoekstra, and J. Schaye, astro-ph/9901191.
MNRAS434, 148 (2013), arXiv:1210.7303. [33] A. Petri, “The lenstools computing package,” (2015).
[11] M. White, Astroparticle Physics 22, 211 (2004), [34] M. Sato, T. Hamana, R. Takahashi, M. Takada,
astro-ph/0405593. N. Yoshida, T. Matsubara, and N. Sugiyama,
[12] H. Zhan and L. Knox, ApJ 616, L75 (2004), ApJ 701, 945 (2009), arXiv:0906.2237 [astro-ph.CO].
astro-ph/0409198. [35] K. Heitmann, D. Higdon, M. White, S. Habib,
[13] A. R. Zentner, D. H. Rudd, and W. Hu, B. J. Williams, E. Lawrence, and C. Wagner,
Phys.Rev.D 77, 043507 (2008). ApJ 705, 156 (2009), arXiv:0902.0429 [astro-ph.CO].
[14] D. H. Rudd, A. R. Zentner, and A. V. Kravtsov, [36] M. Kilbinger, K. Benabed, J. Guy, P. Astier, I. Tereno,
ApJ672, 19 (2008), astro-ph/0703741. L. Fu, D. Wraith, J. Coupon, Y. Mellier, C. Balland,
[15] X. Yang, J. M. Kratochvil, S. Wang, F.R.Bouchet,T.Hamana,D.Hardin,H.J.McCracken,
E. A. Lim, Z. Haiman, and M. May, R. Pain, N. Regnault, M. Schultheis, and H. Yahagi,
Phys.Rev.D 84, 043529 (2011), arXiv:1109.6333. A&A 497, 677 (2009), arXiv:0810.5129.
[16] J. M. Kratochvil, E. A. Lim, S. Wang, [37] M. D. Schneider, L. Knox, S. Habib,
Z. Haiman, M. May, and K. Huffen- K. Heitmann, D. Higdon, and C. Nakhleh,
berger, Phys. Rev.D 85, 103513 (2012), Phys. Rev.D 78, 063529 (2008).
arXiv:1109.6334 [astro-ph.CO]. [38] S. Matsumoto, ArXiv e-prints (2010),
[17] A. Petri, Z. Haiman, L. Hui, M. May, and arXiv:1004.4717 [math.ST].
J. M. Kratochvil, Phys. Rev.D 88, 123002 (2013), [39] The details of the calculation are given in AppendixA.
10
[40] A. Taylor and B. Joachimi, MNRAS442, 2728 (2014), cally,anda perturbativeexpansionis necessary. Writing
arXiv:1402.6983. Ψˆ = Ψ+δΨˆ, we can expand eq. (8) in powers of δΨˆ.
[41] The details of these calculations are given in Appendix Theexpectationvalueofeachterminthisexpansioncan
B. becalculatedintermsofmomentsoftheinverseWishart
[42] L. Casarini, O. F. Piattella, S. A. Bonometto, and
distribution. Ref. [38] provides a general framework to
M. Mezzetti, ApJ812, 16 (2015),arXiv:1406.5374.
compute these moments, and give exact expressions for
[43] M. Takada and D. N. Spergel,
moments up to quartic order. First, let us expand the
MNRAS441, 2456 (2014), arXiv:1307.4399.
inverse of the Fisher matrix estimator (eq. 9) in powers
[44] J. Hartlap, P. Simon, and P. Schneider,
A&A464, 399 (2007), arXiv:astro-ph/0608064. of δΨˆ. The n–th order of this expansion will be
[45] K.Heitmann,N.Frontiere,C.Sewell,S.Habib,A.Pope,
H. Finkel, S. Rizzi, J. Insley, and S. Bhattacharya,
ApJS 219, 34 (2015), arXiv:1411.3396. δFˆ−1 =(−1)n(F−1δFˆ)nF−1 (18)
(n)
with
APPENDIX A: CUBIC AND QUARTIC
COVARIANCE FLUCTUATIONS δFˆ =d′TδΨˆd′ (19)
0 0
The goal of this appendix is to give a derivation of Using eq. (18), we can expand eq. (8) to an arbitrary
eq. (10). When the simulatedfeaturevectordˆ isdrawn order in δΨˆ, take the expectation values of the fluctu-
r
from a Gaussian distribution, the covariance estimator ations over the inverse Wishart distribution, and finally
Cˆ follows the Wishart distribution, and its inverse Ψˆ arrive at eq. (10). We use the notation ν ≡ N −1 and
r
follows the inverse Wishart distribution (see Ref. [29] γ ≡(ν−N −1)/2, and we indicate with capital letters
b
for analytical expressions for these probability distribu- pairs of matrix indices, for example I = (i ,i ), where
1 2
tions). Computing expectation values of eq. (8) over i =1..N . The main results we utilize from ref. [38] re-
a b
the inverse Wishart distribution is not possible analyti- gardingthefirstfourmomentsare(uptoorderO(1/ν2))
ν
hΨˆ i= Ψ (20)
I I
2γ
ν2Ψ Ψ +ν2γΨ Ψ
hδΨˆ δΨˆ i= I J {I J} (21)
I J 4γ2(γ−1)(2γ+1)
ν3Ψ Ψ Ψ
hδΨˆ δΨˆ δΨˆ i= {I J K} (22)
I J K
8γ(γ−1)(γ−2)(γ+1)(2γ+1)
ν4(2γ2−5γ+9)Ψ Ψ Ψ Ψ
hδΨˆ δΨˆ δΨˆ δΨˆ i= {I J} {K L} . (23)
I J K L
16γ(γ−1)(γ−2)(γ−3)(2γ−1)(γ+1)(2γ+1)(2γ+3)
Herethecurlybracketnotationisashorthandforasym- expansionofeq.8),weneedtoapplyanadditionalfactor
metrization over pair of indices: for example of (2γ/ν)n to eqs. (20–23), where n is the order of the
moment up to which we are applying the correction. If
we limit ourselvesto computing the expectationvalue of
Ψ Ψ =Ψ Ψ +Ψ Ψ (24)
{I J} i1j1 i2j2 i1j2 i2j1 eq. (8) up to order O(1/ν2), we do not need to worry
about this correction for eqs. (22–23), as the dominant
Eq. (20) expresses the bias in the Ψˆ estimator that al- termhereisalreadyO(1/ν2). Thenextstepisexpanding
readyappearsintheliterature[44]. Ifwewanttousethe eq. (8) in powers of δΨˆ up to fourth order: this is easily
bias-correctedΨˆ estimator(requiredfortheperturbative done:
4 4
Σˆ = F−1+ δFˆ−1 d′T(Ψ+δΨˆ)C(Ψ+δΨˆ)d′ F−1+ δFˆ−1 . (25)
p (n) 0 0 (n)
! !
n=1 n=1
X X
Carrying out the calculations is simpler than it looks: expansion is proportional to Σ f (N ,N )/Na, where
p a b p r
because of the structure of eq. (25), each term in the