Methods for Bayesian power spectrum inference with galaxy surveys

Methods for Bayesian power spectrum inference with galaxy surveys


We derive and implement a full Bayesian large scale structure inference method aiming at precision recovery of the cosmological power spectrum from galaxy redshift surveys. Our approach improves over previous Bayesian methods by performing a joint inference of the three dimensional density field, the cosmological power spectrum, luminosity dependent galaxy biases and corresponding normalizations. We account for all joint and correlated uncertainties between all inferred quantities. Classes of galaxies with different biases are treated as separate sub samples. The method therefore also allows the combined analysis of more than one galaxy survey. In particular, it solves the problem of inferring the power spectrum from galaxy surveys with non-trivial survey geometries by exploring the joint posterior distribution with efficient implementations of multiple block Markov chain and Hybrid Monte Carlo methods. Our Markov sampler achieves high statistical efficiency in low signal to noise regimes by using a deterministic reversible jump algorithm. This approach reduces the correlation length of the sampler by several orders of magnitude turning the otherwise numerical unfeasible problem of joint parameter exploration into a numerically manageable task. We test our method on an artificial mock galaxy survey, emulating characteristic features of the Sloan Digital Sky Survey data release 7, such as its survey geometry and luminosity dependent biases. These tests demonstrate the numerical feasibility of our large scale Bayesian inference frame work when the parameter space has millions of dimensions. The method reveals and correctly treats the anti-correlation between bias amplitudes and power spectrum, which are not taken into account in current approaches to power spectrum estimation, a 20 percent effect across large ranges in k-space. In addition, the method results in constrained realizations of density fields obtained without assuming the power spectrum or bias parameters in advance.


With the advent of large three dimensional galaxy surveys the three dimensional large scale structure has become a major source of knowledge to understand the homogeneous and inhomogeneous evolution of our Universe. In particular, two-point statistics of the three dimensional matter distribution, are important tools to test our current picture of inflation and evaluate cosmological models describing the origin and evolution of the Universe. For example, detailed knowledge on the overall shape of the matter power spectrum can provide constraints on neutrino masses or the primordial power spectrum and break degeneracies in cosmological parameter estimation from Cosmic Microwave Background (CMB) data when measuring the parameter combination [?]. Since the physics governing the baryon acoustic oscillations (BAO) in the power spectrum is well understood [?], precise measurements of the BAO will allow us to establish this length scale as a new standard ruler to measure the geometry of our Universe through the redshift distance relation [?]. Inferring the BAOs from observation therefore constitutes an important task for determining the detailed nature of a possible dynamic dark energy component, which may explain the currently observed accelerated expansion of the Universe.

Precision cosmology of this kind requires not only the large and informative data sets which are becoming available, such as the Sloan Digital Sky Survey (SDSS), Dark Energy Survey (DES), Canada-France-Hawaii Telsecope Legacy Survey (CFHTLS) and EUCLID, but also fair and accurate data analysis methods [?]. Establishing contact between observations and theoretical predictions in large scale structure is non-trivial. This is due to the fact that generally galaxy redshift surveys are subject to a variety of systematic and statistical uncertainties that arise in the observational process such as, noise, survey geometry, selection effects, and close pair incompleteness due to fiber collisions.

Furthermore, the unknown clustering bias of any one considered galaxy sample compared to the overall mass distribution adds an intrinsic level of uncertainty [?]. Redshift space distortions can be considered a systematic for some analyses but also contain cosmological information.

A careful treatment of the survey geometry is essential for the analysis of LSS surveys [?]. Generally, the raw power spectrum, estimated from a galaxy redshift survey, yields an expectation value that is the true power spectrum convolved with the survey mask [?]. This convolution distorts the shape of the power spectrum and dramatically reduces the detectability of subtle features such as the BAOs. Also the details of galaxy clustering can have a major impact on the shape of the power spectrum [?]. The details of how galaxies cluster and trace the gravitational potentials of dark matter are complicated and not conclusively understood at present. Most likely the relation between galaxies and matter is essentially non-linear and possibly stochastic in nature, such that the inferred galaxy power spectrum is expected to deviate from that of mass [?]. The issue of the galaxy bias remains complicated even in the limit of a linear bias, since different types of galaxies exhibit different clustering behavior [?]. As an example, luminosity dependent clustering of galaxies introduces scale dependent biases, when inferring power-spectra from flux limited galaxy surveys [?].

1.1The state of the art

To address these problems, a large variety of different power spectrum inference methods have been proposed in the literature. Substantial work has been done in respect to Fourier transform based methods, such as the optimal weighting scheme. In this approach galaxies are assigned a weight, in order to reduce the error in the estimated power [?]. As alternatives to Fourier space methods, there exist methods relying on Karhunen-Loève or spherical harmonics decompositions [?]. Furthermore, there exists a rich body of research based on a variety of likelihood methods aiming at inferring the real space power spectrum [?]. To not only provide maximum likelihood estimates but also corresponding conditional errors, [?] proposed a Markov Chain Monte Carlo approach.

1.2Bayesian analysis

Inference of the three dimensional density field and the corresponding cosmological power spectrum from galaxy redshift surveys with highly structured survey geometries and statistical uncertainties is an ill-posed problem. This means, there generally exists a large range of feasible solutions consistent with observations. The Bayesian approach solves this problem by providing numerical representations of the joint three dimensional density and power spectrum posterior distribution [?]. Similar approaches have also been applied to Bayesian inference with CMB data [?].

In the age of precision cosmology the goal is to assess joint and correlated uncertainties between quantities inferred from a given galaxy redshift survey. A particularly important example are galaxy biases that, if measured incorrectly, have the potential to distort the shape of the inferred power spectrum over a large range of modes [?]. Traditionally, biases are accounted for in the following sequential pipeline process: first, measure the bias amplitudes from power-spectra of different galaxy populations fixing to a fiducial power spectrum; then infer the power spectrum from the entire catalog [?]. This approach leaves open the charge of overusing the information content of the data because biases and the power spectrum inferences are treated as independent when they are not. Alternative approaches like the optimal weighting scheme or the Bayesian power spectrum inference method assume fixed bias models [?].

In all these methods uncertainties in the inferred bias models do not feed back to the power spectrum. It is therefore important to quantify these joint uncertainties and to devise methods to account for these effects. In this paper we report on the inclusion of these effects in the inference framework described in [?].

Our approach accounts for joint uncertainties by exploring the joint posterior distribution of the three dimensional density field, the power spectrum, galaxy biases and corresponding normalizations conditioned to observations. Numerical and statistical efficient exploration of the corresponding high dimensional parameter space is achieved via new implementations of multiple block Markov chain and Hybrid Monte Carlo methods. The schematic outline of the iterative block sampling scheme is given in figure ?. Running through the analysis cycle generates a sampled representation of the full joint posterior distribution, which naturally accounts for joint and correlated uncertainties involved in the inference process.

This paper also describes radical methodological changes in the implementation of the sampling scheme underlying the exploration of this joint posterior density. Using a combination of an efficient Hybrid Monte Carlo sampling method for all the parameters, together with the reversible jump algorithm described by [?] we obtained an algorithm that proves to sample all the parameters efficiently both in the high and in the low signal-to-noise regimes.

The paper is structured as follows. In Section 2, we summarize the notation used in this work. Section 3 describes the large scale structure posterior distribution, employed to infer the large scale three dimensional density field, power spectrum, galaxy biases and normalizations from galaxy observations. Here we also describe how the complex task of sampling the full joint posterior distribution can be dissected into a sequence of smaller sub problems within the framework of multiple block Metropolis-Hastings sampling. In Section 4, we discuss the numerical implementation of our method. The following Section 5 describes the generation of an artificial galaxy mock observation which emulates characteristic features of the Sloan Digital Sky Survey Data Release 7 [?]. In Section 6, we perform a test of our method by applying it to artificial galaxy data. These tests particularly aim at estimating the efficiency of the algorithm in a realistic scenario. Some results of the performed test are discussed in Section 7. In particular, here we show, that the posterior correlations between inferred power-spectra and galaxy biases are strongly anti-correlated over large ranges in Fourier-space. The paper is concluded in Section 8 by a discussion of the results obtained in this work and some concluding remarks.


In this section, we describe the basic notation used throughout this work. Let the quantity be the field amplitude of the three dimensional field at position . Then the index has to be understood as a multi index, which labels the three components of the position vector:

where is the th component of the th position vector. Alternatively one can understand the index as a set of three indices so that for an equidistant grid along the three axes the position vector can be expressed as:

with , and being the grid spacing along the three axes. With this definition we yield:

Also note that any summation running over the multi index is defined as the three sums over the three indices , and :

Further, we will frequently use the notation , which denotes the set of field amplitudes at different positions . In particular:

where is the total number of position vectors or voxels in the three dimensional volume.

3The Large scale structure Posterior

In this section we will describe the Large scale structure posterior distribution employed in this work. As will be demonstrated, the multiple block sampling approach permits to dissect the complex problem of exploring a joint posterior distribution into a set of smaller subtasks.

3.1The joint posterior distribution

The aim of this work is to update and complement the previously presented power spectrum sampling framework described in [?] to be numerically more efficient and also to account for additional systematics and statistical uncertainties of the galaxy distribution from which the density field and the power spectrum shall be inferred. In particular, the previously developed algorithm accounted for systematics such as arising from survey geometry and selection effects as well as statistical uncertainties due to noise in the galaxy distribution and cosmic variance [?]. Additional uncertainties arise from the galaxy bias and the correct normalization for the galaxy density contrast. Both of these effects can greatly effect the recovery of the large scale power spectrum and may yield erroneous results when not accounted for accurately. In particular, a galaxy sample consisting of many galaxy types, all tracing the underlying density field differently, does not provide a homogeneous sample of the underlying density field. This is because luminous galaxies will dominate the galaxy sample at larger distances while fainter galaxies are found in regions closer to the observer. Due to the different clustering behavior of luminous and faint galaxies, this observational effect introduces a spatially varying bias, which translates to a scale depended bias in Fourier-space [?]. In addition, a wrong assumption on the galaxy density contrast normalization will yield a non-vanishing mean resulting in erroneous large scale power when estimating power-spectra. It is therefore generally crucial to account for these effects when inferring the density field and power spectrum from realistic galaxy surveys.

Traditional approaches generally measure the power spectrum and the galaxy bias parameters in separate steps by either measuring it from the same data set or employing fitting functions calibrated on other galaxy catalogs [?]. The risk of such approaches is that they may over-use the data by ignoring joint and correlated uncertainties for inferred biases and power-spectra, leading to wrong conclusions on the accuracy of inferred quantities and possibly biased results because not all parameters are explored jointly.

In principle there may also be cases where a joint inference of power spectrum and galaxy biases may be mutually beneficial. Traditional approaches would miss such opportunities.

In this paper we aim for the full characterization of the joint posterior of different tracer populations labeled by the index . We will achieve this by simultaneously exploring the three dimensional density field , the corresponding power spectrum , the galaxy biases and the galaxy density contrast normalizations given the data – galaxy observations in the form of three dimensional number counts .

3.2Dissecting the posterior distribution

Direct analysis of the joint posterior is challenging: it is non-Gaussian and very high-dimensional. From a technical point of view, it is also not advisable to directly generate random variates from the joint posterior distribution since this may result in a numerically unfeasible algorithm preventing rapid exploration of the parameter space. In such a situation one usually has to rely on Metropolis-Hastings algorithms, which reduce the problem of numerical parameter space exploration to the design of suitable proposal distribution for the generation of candidate solutions. A particular important theorem on Metropolis Hastings block sampling permits to dissect the high dimensional sampling problem into a sequence of lower dimensional problems once conditional distributions for the subtasks can be formulated [?]. It is therefore possible to subdivide the exploration of the full joint parameter space into the following sequence of conditional sampling procedures:

where labels the sample number. Iterating these individual sampling steps will then provide samples from the full joint posterior distribution [?]. Also see figure ? where we have illustrated this multiple block sampling procedure, also clarifying the flow of information in subsequent iterations. In the following we will discuss the individual contributions to the full joint posterior distribution as listed in equation (Equation 6).

3.3The density posterior

We now derive the posterior distribution to infer the three dimensional density field [?].

The density prior distribution

This work mainly concerns the inference of the largest scales, where the density field can be reasonably well described by Gaussian statistics. A Gaussian prior will therefore be an adequate description of the statistical behavior of the large scale density field. It should nevertheless be remarked, that a Gaussian prior resembles the least informative prior, once the mean and the covariance of the density field are specified. For this reason, from a Bayesian perspective, the Gaussian prior is well justified approximation even for inference of the density field in the non-linear regime. Therefore, in the following, we will assume the prior distribution for the density contrast to be a multivariate normal distribution with zero mean and the Fourier transform of the covariance matrix being the cosmological power spectrum . With these definitions the prior distribution for the density contrast can be written as:

For further discussions the interested reader is referred to [?].

The density likelihood and the data model

As already discussed above, it is crucial to account for different clustering behavior of galaxy populations in the large scale structure sample from which density fields and power-spectra are inferred. In order to do so, we split the galaxy sample into sub samples for which we can treat the respective systematic and statistical uncertainties. The aim of this work is to present a method which can extract joint information on the three dimensional density contrast field , the power spectrum , the bias and the mean number of galaxies in the survey from a set of galaxy samples, all tracing the underlying density field differently, while properly accounting for the individual systematic and stochastic uncertainties. We present the method in a general fashion that includes joint cosmological analyses from two or more different galaxy surveys or even different probes of the large scale structure, while accounting for their respective systematic and stochastic uncertainties. In the case of sub samples or galaxy observations, a linear data model for each of the subsets can be formulated as:

where are the observed number counts of galaxies in the th sample and in the th volume element, is the corresponding expected number of galaxies, is the corresponding survey response function, which accounts for survey geometry and selection function, is a linear galaxy bias, is the linear growth function and is a random noise component. Also note, that here we are mainly interested in recovering the largest scales and so a linear bias should suffice. We will defer the inclusion of more complicated non-linear and non-local biases to a future work. The survey response operator is given by the product of the angular mask or completeness function and the radial selection function at the th volume element:

By following the arguments in [?], we will assume the random noise component to follow a Gaussian distribution with zero mean and noise covariance matrix , given as:

where is the Kronecker delta. Equation (Equation 9) resembles the Covariance matrix of a Poisson process when averaged over realizations of the density contrast , since the average of vanishes [?]. Alternatively, the assumed Gaussian uncertainty model can be derived by a Taylor-expansion of the Poissonian log-likelihood to second order in the density field. This demonstrates, that the Gaussian likelihood model is correct at large scales, where density amplitudes are small, while it describes an approximation to small scale uncertainties. Although this Gaussian model assumes independence of the noise from the underlying density field, the number counts of the tracer populations are correlated in spatially overlapping regions as they all refer to the same underlying three dimensional density field . In this fashion, the data model given in equation () emulates the expected behavior of galaxy formation models. In this work, we are focusing at inferring the largest scales from observations, for which our data model is adequate. Note, that precise inference of the large scale structure in the non-linear regime while accounting for the full Poissonian noise statistics can be performed with the method presented in [?] and [?].

With these definitions the likelihood for the dataset can be expressed as a multivariate normal distribution:

The joint likelihood from all datasets is then obtained as:

where we assumed conditional independence of the individual likelihoods, once the density contrast , the mean galaxy numbers and the biases are given.

3.4The Power-spectrum posterior distribution

As described above, sampling from the Wiener posterior distribution will provide realizations of full three dimensional density fields conditioned on the observations. Basing on these density fields, in subsequent steps realizations of the corresponding power spectrum can be generated by sampling from the conditional power spectrum posterior distribution:

where we assumed conditional independence of the data on the power spectrum, once the full three dimensional density field is given, and in the last line we used the density prior given in equation (Equation 7). The result provided by equation (Equation 12) demonstrates that once a realization of the three dimensional density field is given, inferring a realization of the power spectrum is conditionally independent on the data. This fact was first pointed out by [?] in case of CMB analysis and was later adapted for large scale structure inference by [?]. This is a particularly important result for the power spectrum sampling procedure. Since all observational systematics and uncertainties have already been accounted for in inferring a density field realization, this considerably simplifies conditional power spectrum inference [?]. Also note, that equation (Equation 12) is a standard result for Bayesian density and power spectrum inference with Wiener posteriors in the case of CMB and large scale structure analyses [?]. It is interesting to remark that equation (Equation 12) provides a direct sampling procedure to generate power spectrum realizations. Following [?], equation (Equation 12) can be re expressed as an inverse Gamma distribution, when employing a power-law prior for the power spectrum amplitudes. Thus, inserting a power law prior in equation (Equation 12) and imposing the correct normalization, the conditional power spectrum posterior can be written as [?]:

where is the power spectrum amplitude of the th bin corresponding to the Fourier mode , is the number of modes in that bin, is the power-law index of the power-law prior, is the Gamma function and is given as the squared isotropic average over the mth mode bin in Fourier space :

with being the Fourier transform of the density field. Note, that the power-law prior with describes a “flat” prior on a unit scale, while yields a Jeffreys Prior. Jeffreys prior is a solution to a measure invariant scale transformation, and hence is a scale independent prior, as different scales have the same probability. For this reason, Jeffreys prior is optimal to infer cosmological power-spectra, which constitute scale measurements, as it does not introduce any bias on a logarithmic scale. For a more detailed derivation and discussion of the Inverse Gamma distribution in case of LSS inference the reader is referred to our previous work presented in [?].

3.5The bias posterior distribution

Accounting for galaxy biases is particularly important when analyzing galaxy samples exhibiting widely varying clustering behavior across galaxy types. These galaxy biases give rise to additional systematics, which when ignored, have the potential to distort the shape of inferred power-spectra. We extend the power spectrum sampling framework described in [?] with a galaxy bias sampling procedure. To be maximally agnostic about the functional shape of the galaxy bias as a function of galaxy properties such as luminosity or color, we will assign a galaxy bias to each of the individual galaxy sub samples. Another approach would be to assume a functional shape for the galaxy bias as a function of some galaxy properties such as luminosity and just fit the corresponding function parameters. Such an approach would be a posteriori and hence overuse the data since knowledge of such a functional shape is generally extracted from the galaxy sample under consideration itself. We will instead pursue an approach that only assumes that bias properties remain approximately constant within luminosity bins.

Realizations of the bias factors can then be obtained by sampling from the conditional bias posterior distribution given as:

where the bias factors are conditionally independent of the power spectrum, once the full three dimensional density field is known. Using Bayes theorem we can rewrite the conditional bias posterior as:

where we introduced a prior factorizing in the different galaxy samples and assumed conditional independence . Since the probability distribution given in equation (Equation 14) factorizes in the galaxy sub samples, one can estimate the bias factors for each galaxy sample separately by sampling from the following distribution:

In order to follow the maximum ignorance principle also with regard to the amplitude of the bias, we assume a flat prior . With these assumptions the conditional bias posterior distribution is simply proportional to the density likelihood given in equation (Equation 10):

As can be seen from equation (Equation 15) and as demonstrated in Appendix B, this distribution turns into a univariate normal distribution for the bias factor given as:

with given as:

and the mean given as:

A random realization for the galaxy bias can thus be obtained by simply sampling from a univariate normal distribution. Note, that while previous methods such as presented in [?] measure the relative bias from estimated power-spectra, our method infers the bias parameters directly from the relation between the data and the three dimensional density field. Further, since the density field contains information inferred jointly from all galaxy sub samples, inference of individual bias factors depends on all other biases via the density field, enabling us to measure relative biases and thus the shape of the power spectrum up to an overall normalization.

3.6The posterior distribution

Of similar importance as galaxy biases are the mean numbers of galaxies for the sub samples. These mean numbers of galaxies are crucial to define the true underlying density contrast of the galaxy distribution. A false value of will yield a non-vanishing mean and result in erroneous large scale power in the inferred power-spectra. Since the are not known a priorly, they will also have to be inferred jointly, in order not to miss possible cross-correlations and interdependencies. The parameters for the sub samples can then be inferred from a conditional distribution given as:

where, in a similar fashion as described in Section 3.5, we assume conditional independence of the power spectrum once the full three dimensional density field is given. Further, following similar arguments as described in Section 3.5, it is possible to infer each parameter separately, since the corresponding probability distribution factorizes. As a result we obtain the following probability distribution for the parameter :

Note, that also here we followed a maximum ignorance approach by setting the prior . As can be seen the resultant probability distribution for the parameter is proportional to the density likelihood given in equation ( ?). Following the discussion presented in Appendix A, the probability distribution for the given in equation (Equation 40) is a generalized inverse Gaussian distribution (GIG) given as:

where is a modified Bessel function of the second kind, , is given as:

and being:

The number of observed voxels is given by:

where the Heaviside function counts the non zero elements of the three dimensional survey mask. Random realizations for the parameters can then be provided by sampling from the GIG distribution. It is interesting to note, that the inferred also depend on the galaxy bias parameters and the three dimensional density field. This indicates, that a joint estimation of all these parameters is required in order to correctly account for joint information and correlated uncertainties. Algorithms to sample from the GIG distribution are described in literature [?].

4Numerical implementation

In this section we will describe the numerical implementation of the large scale structure sampler.

4.1Sampling the density field

As is manifest from the above discussion, the three dimensional density sampling step resembles the core of our method. The Wiener posterior for the density sampling step described in Section 3 is a multivariate normal distribution, and as such possesses a clear sampling procedure, by first calculating the Wiener-mean and adding a variance term corresponding to the Wiener-variance [?]. Nevertheless, calculating the Wiener mean and variance requires matrix inversions for large matrices with typically on the order of entries. These matrix inversions can be performed with efficient operator based implementations of conjugate gradient methods [?]. Alternatively, the Wiener posterior in case of CMB data has also been explored using a Hybrid Monte Carlo (HMC) algorithm [?]. Furthermore, the HMC has been successfully applied for the inference of three dimensional density fields in the non-linear regime [?]. Since these references have found the HMC algorithm to be very efficient in large scale structure problems we implement an HMC algorithm for the exploration of the large scale structure density field. Compared to the direct sampling approach based on iterative techniques, HMC is less sensitive to the effects of numerical inaccuracies, e.g. due to insufficiently strict tolerances in the inversion scheme, because the final accept/reject step of the Metropolis Hastings sampler ensures correctness of the target density. Any such inaccuracies would only impact efficiency but not correctness.

In the following we will review the basics of the HMC algorithm required for this work.

The HMC algorithm

The numerical efficiency of the Hybrid Monte Carlo algorithm originates from the fact that it exploits techniques developed to follow classical dynamical particle motion in potentials to deterministically provide new proposals to the Metropolis-Hastings algorithm [?]. Assuming we wish to generate random variates from a probability distribution , with being a set consisting of the elements , we may interpret the negative logarithm of this distribution as a potential :

Further, by introducing a ’momentum’ variable and a ’mass matrix’ , one can add a kinetic term to the potential in order to obtain a Hamiltonian function, which, similar to ordinary mechanics, describes the dynamics in a high dimensional parameter space. Such a Hamiltonian can be written as:

The new target probability distribution is then obtained by exponentiating the Hamiltonian given in equation (Equation 19):

As can be seen from equation (Equation 20), the form of the Hamiltonian was chosen such, that the new target distribution is separable into a Gaussian distribution in the momenta and the target distribution . This immediately clarifies that the two sets of variables and are independent and marginalization over the artificial momenta yields the original target distribution .

Given an initial point in this high dimensional phase space it is possible to deterministically propose a new sample from the target distribution by following persistent trajectories for some fixed time according to Hamilton’s equations:


A new point in phase space will be accepted according to the standard Metropolis-Hastings acceptance rule:

Since Hamiltonian dynamics conserve the Hamiltonian, the HMC approach results in an acceptance rate of unity. In practical applications, however, the acceptance rate may be lower due to numerical inaccuracies in the integration scheme. Once a new sample has been accepted the momentum variable is discarded and the process restarts by randomly drawing a new set of momenta. The individual momenta will not be stored, and therefore discarding them amounts to marginalizing over this auxiliary quantity. Hence, the Hamiltonian sampling procedure basically consists of two steps. The first step is a Gibbs sampling step, which yields a new set of Gaussian distributed momenta. The second step, on the other hand amounts to solving a deterministic dynamical trajectory on the posterior surface.

Equations of motion for the Wiener Posterior

As described above, in this work we will employ a Wiener posterior for the inference of the three dimensional density field. In this case the Hamiltonian potential can be written as:

where we omitted terms corresponding to normalization constants. Variation of with respect to the three dimensional density field then yields the Hamiltonian forces which govern the dynamics in parameter space:

With this definition the equations of motion can be given as:


As can be seen from equation (Equation 27), in the framework of the HMC the quadratic form of the Wiener posterior yields a multidimensional coupled harmonic oscillator whose zero points are defined by the data. Also note, that equations (Equation 26) and (Equation 27) can be combined to yield a second order differential equation:

This clarifies that an ideal choice for the mass matrix would be given by:

which would decouple the harmonic oscillators and permit a trivial analytic solution to the equations of motion for a set of independent harmonic oscillators. Nevertheless, this approach would again require matrix inversions and hence be identical to a Gibbs sampling approach, which we try to avoid here [?]. In order to provide samples from the target distribution, we will therefore integrate the equations of motion (Equation 26) and (Equation 27) numerically over a pseudo time interval via a leapfrog scheme. In order to avoid resonant trajectories, the pseudo time interval is drawn from a uniform distribution. It is also important to remark, that the leapfrog scheme is a symplectic integrator, and as such ensures detailed balance of the Markov chain. For a detailed description of the numerical implementation of the HMC in general and in cosmology in particular the reader is referred to the literature [?]. In general the choice of the mass matrix does not affect the validity but only the numerical efficiency of the sampling method. We choose a diagonal form consisting of the diagonal elements of the full mass matrix, given in equation (Equation 29), in its Fourier representation. This choice enhances performance while avoiding matrix inversions, analogous to using a diagonal preconditioner in iterative conjugate gradients solvers.

4.2Sample the power spectrum

The conditional power spectrum posterior, given by the inverse gamma distribution in equation (Equation 13) , facilitates a direct sampling scheme [?]. In particular, by performing a change of variables one obtains the -distribution

where . Consequently, generating random variates for the power spectrum reduces to the task of drawing independent random samples from the -distribution for each Fourier mode . Let be independent, normally distributed random variates with zero mean and unit variance. Then [?]:

is -distributed, with a element vector with elements . A sample of the power spectrum amplitude at Fourier mode is then obtained by:

It should be remarked that each is a positive quantity, ensuring the requirement that the signal covariance matrix be a positive definite quantity.

4.3Exploring the low signal to noise regime

In a Poissonian picture of galaxy formation the signal-to-noise ratio can be expressed as the square root of the number of galaxies found in a given volume . Consequently, splitting a galaxy sample into a number of sub samples, to accurately account for the respective systematics will yield a lower signal-to-noise ratio for each of these galaxy sub samples. Although, the Likelihood given in equation (Equation 11) correctly accounts for the joint recovery of information from all galaxy samples by correctly weighting them corresponding to their respective noise properties, the overall signal-to-noise ratio may still be low, depending on the amount of sub samples. This effect may influence the numerical efficiency for the power spectrum sampling procedure as described above.

While the approach described in sections Section 4.1 and Section 4.2 provably converges to the target posterior on all scales probed by the data, it may take prohibitive computational time to generate a sufficient amount of independent samples in low signal to noise regime [?]. In particular, the variations in subsequent power spectrum samples are solely determined by cosmic variance, whereas the full posterior distribution is governed by cosmic variance and noise. This permits a rapid exploration of parameters at the largest scales, where cosmic variance is dominant, but yields poor statistical efficiency for the small scale regime, where cosmic variance is typically much smaller than the noise variance. This results in a prohibitively long correlation length of the sequence of spectra in the low-signal to noise regime, requiring unfeasible long Markov Chains to generate sufficient numbers of independent samples [?]. The above discussion shows that the variables chosen to explore the parameter space are optimal for the cosmic variance dominated regime, but it may also suggest that a change of variables is advantageous for the exploration of the low signal to noise regime.

In fact, a change of variables can greatly improve the mixing behavior of the Markov Chain. Following [?], we introduce a change of variables which corresponds to the following deterministic transition for the three dimensional density field in Fourier space:

while, as demonstrated in Appendix C, a new power spectrum sample has to be generated from the distribution:

with and labeling the sample number in the Markov chain. Unfortunately there exists no immediate way to generate power spectrum samples from this distribution. While [?] propose a simple Metropolis-Hastings sampler to generate samples from the distribution given in equation (Equation 34), we will employ a HMC to solve the problem. As already discussed above, the HMC has the advantage that it greatly improves the acceptance probability compared to standard implementations of a Metropolis-Hastings algorithm. Furthermore, the entire procedure of sampling a power spectrum from the distribution given in equation (Equation 34) followed by the deterministic transition, described by equation (Equation 33), for the density field is numerically as expensive as the sampling procedure outlined in sections Section 4.1 and Section 4.2. Alternating these two sampling procedures will thus not decrease the overall numerical efficiency of the algorithm. On the other hand, alternated sampling in these two variable bases greatly improves the statistical efficiency and hence the overall performance of our algorithm.

5Generating artificial observations

Above we discussed the derivation and numerical implementation of our method. In this section we will describe the generation of an artificial galaxy survey, used to test our method. Following closely the description in [?], at first a realization for the density contrast is generated from a normal distribution with zero mean and a covariance matrix corresponding to a cosmological power spectrum. The density field will be stored on a three dimensional equidistant Cartesian grid with , corresponding to volume elements, and a co moving box length of . The cosmological power spectrum for the underlying matter distribution, with baryonic wiggles, is calculated according to the prescription described in [?] and [?]. Further, a standard CDM cosmology with a set of cosmological parameters (, , , , , ) is assumed.

On top of this three dimensional density contrast realization , artificial galaxy sub samples will be generated corresponding to the data model described in section ?. For each subsample the noise term is generated from a multivariate normal distribution with zero mean and a diagonal covariance matrix .

In particular, the aim of the present work is to test our method in a realistic scenario with largely inhomogeneous galaxy samples with respect to their clustering behavior. As already discussed above, such a sample would not represent a spatially homogeneous sample of the three dimensional density field, and hence introduces scale dependent biases in Fourier space [?]. This situation is usually faced when analyzing real galaxy observations [?]. Subdividing the galaxy sample into sub samples permits to account for the respective systematic effects while jointly extracting information on the cosmological power spectrum from a galaxy survey.

In setting up a realistic test scenario, we emulate the main features of the Sloan Digital Sky survey data release 7 as closely as possible [?]. In particular, we use the completeness mask of the Sloan Digital Sky Survey data release 7 depicted in the right panel of figure ?. This mask was evaluated with the MANGLE code provided by [?] and has been stored on a HEALPIX map with [?]. Further, we assume that the radial selection function follows from a standard Schechter luminosity function with standard r-band parameters ( , ), and we only include galaxies within an apparent Petrosian r-band magnitude range and within the absolute magnitude ranges to . The corresponding radial selection function is then obtained by integrating the Schechter luminosity function over the range in absolute magnitude:

where is given in Appendix D. The selection functions for the galaxy sub samples are depicted in the left panel of figure ?. The product of the survey geometry and the selection function at each point in the three dimensional volume yields the survey response operator:

where and are the right ascension and declination coordinates corresponding to the th volume element, and is the corresponding redshift. The functions and are the continuous survey completeness and selection function, respectively, sampled at the th grid position. Further, we will model a luminosity dependent bias, in terms of absolute magnitudes M, as:

with the fitting parameters , and [?]. Since equation (Equation 35) formulates only the relative bias with respect to a reference point in absolute magnitude , we will arbitrarily set and for our test scenario. Finally we have to define a artificial set of mean galaxy numbers . In the case studied here, the different are directly related to the Schechter luminosity function by integrating over the width of the th absolute magnitude bin:

The individual galaxy observations are then obtained by evaluating Equation 8 with the quantities as specified above.

6Testing the algorithm

The aim of this section is to conduct a series of tests to judge the performance of our method in a realistic scenario by applying the algorithm to the simulated galaxy survey described in the previous section. We particularly focus on studying the burn-in and convergence behavior of our method, which determines its overall numerical feasibility and efficiency in real-world applications.

6.1Testing burn-in and statistical efficiency

Metropolis Hastings samplers are designed to have the target distribution as their stationary distribution [?]. The sampling process will provide samples from the specified joint large scale structure posterior distribution after a sufficiently long burn-in phase. The theory of Metropolis Hastings sampling by itself does not provide any criterion to determine the length of the burn-in phase, so we study it in numerical experiments.

Typically, burn-in manifests itself as a systematic drift of sampled parameter values towards their true underlying values from which the artificial data was generated. This behavior can be monitored by following the evolution of parameters in subsequent samples [?]. To test this initial burn-in behavior we will arbitrarily set the initial values for the parameters and to unity , while the initial power spectrum amplitudes for are assumed to be two orders of magnitude larger than their true underlying values. One can then observe the successive drift towards the true values of the respective parameters. The 1000 first successive samples of the Markov chain are presented in figure ?, which clearly demonstrates the drift towards the respective true solutions.

Figure ? already illustrates the interplay between bias and power spectrum samples. While bias parameters approach their true values from the top, the power spectrum approaches its fiducial value from below. This effect originates from the anti-correlation between biases and power-spectra for given data. In particular high bias values will result in lower power spectrum amplitudes and vice versa.

More quantitatively the initial burn-in behavior is studied by comparing the th sample of a parameter in the chain to its true underlying value via:

where for our case, the parameter sample has to be replaced by one of the parameters sampled by the Markov chain. The results for the , and are shown in the lower panels of figure ?. It can be seen that the algorithm drifts towards the preferred region in parameter space and starts exploring the parameters corresponding to their joint uncertainties. Also note, that figure ? already indicates the degree of scatter which can be expected in the final analysis of the Markov run.

In this test we do not observe any particular hysteresis for the poorly constrained large scale modes, meaning they do not remain at their initially set values but efficiently explore the parameter space. This reflects the ability of our algorithm to account for artificial mode coupling as introduced by the survey geometry, selection effects and biases.

These tests show that the algorithm requires an initial phase of on the order of samples before providing typical samples from the joint posterior. This initial burn-in period can be shortened by initializing the algorithm with a set of parameters closer to the true solutions, contrary to the test presented here where the initialization was intentionally chosen to test the burn-in behavior.

Statistical efficiency is a major design goal of Markov Chain Monte Carlo algorithms. Generally, successive samples in the chain are not independent but correlated with their predecessors. Consequently the degree of correlation between successive samples determines the amount of independent samples which can be drawn from a chain with given length. To study this correlation effect we follow a similar approach as described in [?] or [?].

Assuming all parameters in the Markov chain to be independent of each other, the correlation between subsequent samples for the parameter can be quantified in terms of the autocorrelation function

where is the distance in the chain measured in iterations [?]. The results for this analysis are presented in figure ? which shows the auto-correlation functions for the parameters , and .

One can define a correlation length of the Markov sampler as the distance in the chain beyond which the correlation coefficient has dropped below a threshold of . As can be seen in figure ? the correlation length is typically less than 200 samples, except for a few power-spectrum modes in the Nyquist regime of the box, which are correlated across a few hundred samples. This demonstrates the overall high statistical efficiency of the sampler. Note, that these Nyquist modes are under-sampled due to the finite resolution of the Fast Fourier Transform. Even though our method treats these eight Nyquist frequencies correctly, they will not be considered for a typical cosmological analyses where usually only modes with frequencies less than seventy percent of the Nyquist frequency are considered [?].

These tests demonstrate the numerical feasibility of our method to explore the full joint posterior distribution despite the high dimensionality of the problem presented here.

6.2Joint uncertainties and correction of systematics

The previous section demonstrates the numerical and statistical efficiency of our algorithm to explore the joint large scale structure posterior in high dimensional spaces. As will be shown in this section, the algorithm also accurately accounts for systematics and joint stochastic uncertainties in the inference of the different parameter. Of particular importance for the measurement of cosmological parameters from power-spectra are the effects of survey geometry and selection effects. Survey geometry and selection effects introduce artificial mode coupling to the Fourier modes of the three-dimensional density field, and can greatly reduce the visibility of features in the power spectrum, such as the baryonic acoustic oscillations [?]. Correcting these effects is hence of utmost importance for any method to estimate the cosmological power spectrum from observations.

The information on how well our method corrects for survey geometry effects is naturally contained in the correlation structure of the individual power spectrum samples. To examine this effect, we calculate the correlation matrix of the power spectrum samples

where is the power spectrum amplitude at the mode . The result of this exercise is presented in figure ?, where the ensemble average has been taken over power spectrum samples. As can be clearly seen, the resultant correlation matrix exhibits a well defined diagonal structure, as expected from theory. Also note the slight red band of anti-correlation around the diagonal, particularly at the highest frequencies in figure ?. This anti-correlation indicates that the power spectrum frequency resolution is higher than supported by the data. Note that these effects are accurately accounted for by our method, since it provides the full probability distribution for each mode. If one wishes, it is possible to reduce frequency resolution in a post-processing step until anti-correlation vanishes [?].

It should be noted, that generally the posterior distributions for the individual are non-Gaussian. Consequently, two-point statistics of the power spectrum, as presented in figure ?, do not provide conclusive information on the full statistical behavior of the . Since our method provides a numerical representation of the large scale structure posterior higher order statistics for the power spectrum can naturally be taken into account. To illustrate this fact, we plot two dimensional marginalized posterior distributions for a few arbitrarily chosen pairs of s in figure ?.

In the previous Section 6.1 we saw burn-in behavior that suggested bias parameters being anti-correlated with inferred power-spectra . This behavior is easily understood: a higher bias parameter requires a lower amplitude of the power spectrum to match observations. The details and strength of these effects generally depend on the magnitude bin for the bias parameter and the mode of the power spectrum. To quantify these effects we now study the posterior correlations between the inferred bias parameters , the and the power spectrum . Figure ? shows the correlation matrices

where and have to be replaced by the parameter under consideration. The scale dependent anti-correlation between bias parameters and power-spectra is clearly visible in the left panel of figure ?. Note, that we are measuring relative biases and for that reason the bias amplitude at a magnitude is fixed to 1.

The maximum anti-correlation amounts to about 20 per cent, and significant anti-correlation can be observed at all scales and for all magnitudes. This result emphasizes our claim, that even simple linear bias amplitudes cannot be estimated independently from power-spectra and vice versa.

The middle panel of figure ? shows the correlations between the and different modes of the power spectrum. Generally correlations and anti-correlations amount maximally up to a few percent. However, at the magnitude bin of the fixed bias amplitude is strongly anti-correlated with the modes of the power spectrum. This indicates, that carries additional information on the bias amplitude at the magnitude bin .

Finally, the right panel of figure ? shows the correlations between inferred and parameters. The prominent diagonal character of this correlation matrix demonstrates that the anti-correlation between inferred biases and is strongest for values belonging to the same magnitude bin. Nevertheless, there exist percent level correlations and anti-correlations between parameters in different bins. A prominent feature of this plot is also the strong band of positive correlations for the parameter at magnitude , which demonstrates strong correlations between in this bin and the bias parameters of other magnitude bins.

Since our method provides a numerical representation of the highly non-Gaussian posterior distribution, uncertainty quantification is not limited to Gaussian approximation but can take into account all involved non-Gaussianities and non-linearities. As a demonstration, in figure ?, we plot marginalized two dimensional distributions of individual, randomly selected, bias parameters and modes of the power spectrum . It can be seen, that the marginalized distributions reflect the previously observed anti-correlations between bias parameter and power-spectra. Also note, that the marginalized distributions, presented in figure ?, exhibit non-Gaussian behavior. The results of this section support the requirement for a full joint approach to cosmological power-spectrum estimation.


In the previous section we discussed the numerical feasibility of our method and demonstrated the necessity of a full joint approach. Here we will describe the final results inferred from the Markov chain. As described in Section 5, we seek to test our method in a realistic scenario by emulating characteristic features of the SDSS DR7 main sample. A galaxy sample, such as the SDSS main sample, represents an interesting and challenging test scenario for our method, since it consists of a large class of galaxy types exhibiting different clustering behavior at different redshifts and hence requires to exploit the entire methodology as described in this work. On the other hand, the SDSS main sample is a shallow galaxy sample with a median redshift of about , which is not ideal to infer tiny features at large scales in the power spectrum, such as the baryon acoustic oscillations.

Additional information, to probe the large scales, can be obtained from highly clustered luminous red galaxies (LRG) which cover larger volumes to higher redshifts. For the sake of demonstrating the capabilities of our code to infer the power spectrum on the entire range of scales covered by the survey, we ignored LRGs in this work. Note, that LRGs can be easily added as an additional data bin in the joint analysis, where they would add their information content to the inference process. It should be remarked, that the analysis volume of , used for the inference, is only about half filled by data.. This obviously increases the challenge compared to a case where the empty volume would be filled with additional LRGs. In this respect, the presented test is harder compared to a real cosmological analysis in which we would exploit all available information.

To summarize the inferred quantities, in figure ?, we show ensemble means for the power spectrum, the bias and , estimated from Markov realizations. It is immediately obvious, that the mean density parameters can be recovered very well. The middle panel of figure ? depicts the ensemble mean for the bias parameters and corresponding uncertainties. The plot demonstrates that the shape of the magnitude dependent bias function can be accurately recovered from data in this joint analysis. Since we are measuring relative biases the bias parameter at a magnitude of is held fixed at unity.

The left panel of figure ? depicts the ensemble mean power spectrum and corresponding uncertainties. The inferred ensemble mean power spectrum nicely follows the true underlying fiducial power spectrum across the entire range of modes covered by our analysis. Even on the largest scales, poorly constrained by the artificial galaxy survey dominated by sample variance, the analysis returns sensible constraints and hints at the turn-over in the power spectrum.

Visually it is not possible to observe recovered baryon acoustic oscillations. This is due to several effects. First of all, the analyzed volume is too shallow to accurately probe the Baryon acoustic oscillations. Secondly, joint inference of the power spectrum, and correctly accounts for coupled uncertainties which will lead to a blurring of tiny features in the power spectrum compared to an analysis that takes best-fit biases as inputs to power spectrum analysis. We conclude that our artificial data set, emulating the SDSS main sample data, is not sufficiently informative to provide conclusive information on the baryon acoustic oscillations. Although preliminary tests with real observations already indicate qualitatively similar behavior, a detailed analysis of the latest SDSS data will be subject to a forthcoming publication.

Finally, we would like to point out, that our method also provides measurements of the three dimensional density field and corresponding uncertainties. Individual samples of these density fields represent constrained realizations of the underlying density field cleaned from all statistical and systematic uncertainties, and as such will generally be valuable to study galaxy and three dimensional large scale structure formation. Note that in contrast to conventional constrained realizations our method obtains these reconstructions without assuming biases, number densities, or correlations a priori.

The slices through the posterior mean density shown in figure ? can be thought of as a non-linear Wiener filter. Our method also delivers an error model for this reconstruction, here summarized in terms of standard deviation per pixel. These plots illustrate that our method is able to account for arbitrary and complex survey geometries in three dimensional space.

8summary and conclusion

This work describes the derivation and numerical implementation of a full Bayesian approach to cosmological large scale structure power spectrum analysis.

It addresses the problem of joint power spectrum and galaxy bias inference from redshift surveys, and thus aims at a natural and fully self consistent treatment of joint and correlated uncertainties. Traditionally, these correlations remain generally unaccounted for in standard sequential pipeline approaches to power spectrum inference, which assume fixed bias models [?]. The method, as described in this work, naturally accounts for correlated uncertainties between all quantities to be inferred, by exploring the joint parameter space, of the three dimensional density field, the power spectrum, luminosity dependent biases and corresponding normalizations. This is achieved through the use of multiple block Metropolis Hastings and Hybrid Monte Carlo algorithms.

As described in Section 3, the multiple block sampling approach permits us to dissect the computationally challenging task of joint sampling into sequential sub tasks of lesser complexity. Thus, the problem can be processed as outlined in figure ?, yielding a numerical representation of the full joint posterior distribution.

The algorithm presented in this work improves over the method presented in [?], in various respects. Most notably, the previous three dimensional density sampling procedure, relying on the solution of large systems of linear equations via Krylov space methods, has been replaced by a Hybrid Monte Carlo scheme. This has proven to be more efficient, especially in high dimensions.

We further increased the statistical efficiency and hence the number of independent samples which can be drawn from the Markov chain. To do this, we followed the deterministic transition approach proposed by [?], but incorporated an HMC transition, rather than a pure Metropolis-Hastings step, to improve the overall computational efficiency of the algorithm.

Further, we derived conditional posterior distributions for the inference of absolute magnitude dependent bias and corresponding normalizations. As described in Section 3, we followed a maximally agnostic approach by assuming all absolute magnitude bins to be independent and additionally assuming flat priors for the corresponding bias and normalization amplitudes. The resultant conditional posterior distributions have analytic expressions. In particular, realizations of magnitude dependent biases can be generated from univariate normal distributions and corresponding normalization factors are drawn from generalized inverse Gaussian distributions.

We also described tests of our method in a realistic scenario. Section 5 describes the generation of an artificial galaxy survey emulating characteristic features of the SDSS DR7, such as survey geometry, selection effects, and luminosity dependent galaxy biases. The test was particularly designed to highlight the problematics of survey geometry, selection effects and luminosity dependent biases and their impact on cosmological power spectrum inference. Specifically, the main goals of the test, described in Section 6, is to study the expected numerical behavior and efficiency of our method when applied to real data. Through a simple numerical experiment, by setting the initial values of the bias and normalization factors to unity and scaling the amplitude of the initial power spectrum up by two orders of magnitude, we were able to determine the initial burn-in phase of the Markov sampler to be on the order of 1000 samples, which demonstrates numerical feasibility of the method. Also note, that initial burn-in time can be significantly shortened by a more realistic choice of the initial parameter values.

We overcame a crucial difficulty in implementing a joint inference approach: statistical dependency between the bias parameters and the power spectrum leads to unfeasibly long correlation lengths for our block sampling scheme. Adding a deterministic transition step reduces the correlation length of typical parameters by a factor of more than 1000. As a remark, in tests without such a deterministic transitions we observed still eighty percent correlations at a lag of 10000 samples for most of the power-spectrum modes yielding an unfeasibly long convergence time of the sampler. This dramatic improvement allows practical use of our algorithm even in regimes where signal-to-noise is low.

Our method accurately accounts for artificial mode coupling due to the survey geometry and selection effects. The correlation structure of the Markov samples in Fourier space reveals that the power spectrum correlation matrix is of highly diagonal shape. Even when coupled with inference of the bias parameters, maximal residual correlations between power spectrum bins are typically less than 10 per cent, emphasizing the ability of our method to account for survey geometry effects.

Our test revealed anti-correlations at the level of 20 per cent between the bias parameters and the power spectrum across large ranges in Fourier-space. It remains to be studied to what extent this coupling has impacted past cosmological analyses of galaxy survey data that ignored it.

While a linear luminosity dependent galaxy bias is one of the simplest bias models for cosmological inference from galaxy surveys, it does not seem reasonable to assume that more complicated bias models would not suffer from this coupling effect. Still, our method is easily extensible to include more complicated bias models such as second and higher order bias.

In Section 7, we summarize means and credible regions for the parameters inferred by our method. In particular, the power spectrum can be accurately recovered from the test data.

Since we treat each luminosity class as a separate galaxy sample with its own bias parameters, we actually inferred the power spectrum from 15 different galaxy sub-samples all exhibiting different systematic and statistical uncertainties. In particular, each sub-sample comes with its own mask, average number density, bias, and selection function. Our method is therefore able to perform joint power spectrum analysis from combinations of different galaxy surveys—and, after further development, possibly even different probes of large scale structure.

Note, that the normalization parameters in the different magnitude bins are directly related to the luminosity function, which in a scenario as described in this work, is a natural byproduct of the cosmological inference. Also note, that our method further provides full three-dimensional maps of the inferred matter distributions and corresponding uncertainties. These constrained realizations of the density distribution may provide useful information to test cosmological structure formation.

The method, as described here, seamlessly integrates into the larger Bayesian large scale structure inference framework. In particular, it can be used in conjunction with the photometric redshift sampling method described in [?], to account for additional redshift uncertainties. In that combination our method would accurately treat all joint and correlated uncertainties when analyzing galaxy surveys with highly uncertain radial positions of galaxies. Note, although not discussed in this work, the present bias sampling framework seamlessly integrates with effective measures of redshift space distortions corrections, for instance as described in [?], [?] or [?], but also with information theoretically more rigorous approaches, which will be subject to future publications.

In addition, a exciting next step will be to merge our method with the initial conditions inference framework presented in [?]. This would permit inferring the power spectrum of the initial conditions of the Universe from catalogs of biased tracers.

In summary, our method jointly infers all parameters in the cosmological power spectrum inference problem and accurately accounts for correlated uncertainties and interdependencies. As a highly flexible and numerically efficient approach it seamlessly integrates into a larger body of large scale structure inference methods and allows straightforward generalization to more sophisticated biasing schemes. The application of this approach to actual data will be the subject of a forthcoming publication. Preliminary tests indicate a qualitatively similar behavior to the results presented here.


We thank Andrés Balaguera Antolinez and Franz Elsner for useful discussions and support in the course of this project. JJ is partially supported by a Feodor Lynen Fellowship by the Alexander von Humboldt foundation. BDW acknowledges support from NSF grants AST 07-08849 and AST 09-08693 ARRA, and a Chaire d’Excellence from the Agence Nationale de Recherche and computational resources provided through XSEDE grant AST100029. This material is based upon work supported in part by the National Science Foundation under Grant No. PHYS-1066293 and the hospitality of the Aspen Center for Physics.

AInverse Gaussian distribution for

As described in section the posterior distribution for the mean galaxy numbers is proportional to the data likelihood:

Expanding the square in the exponent and discarding numerical constants then yields: