Measuring the effective complexity of cosmological models
Abstract
We introduce a statistical measure of the effective model complexity, called the Bayesian complexity. We demonstrate that the Bayesian complexity can be used to assess how many effective parameters a set of data can support and that it is a useful complement to the model likelihood (the evidence) in model selection questions. We apply this approach to recent measurements of cosmic microwave background anisotropies combined with the Hubble Space Telescope measurement of the Hubble parameter. Using mildly non–informative priors, we show how the 3year WMAP data improves on the firstyear data by being able to measure both the spectral index and the reionization epoch at the same time. We also find that a nonzero curvature is strongly disfavored. We conclude that although current data could constrain at least seven effective parameters, only six of them are required in a scheme based on the CDM concordance cosmology.
pacs:
98.80.Es, 02.50.r, 98.70.VcI Introduction
The quest for a cosmological standard model is being driven by an increasing amount of high quality observations. An important and natural question concerns the number of fundamental parameters of the underlying physical model. How many numbers are necessary to characterize the Universe? Or in other words, how complex is the Universe? Generally the term “complexity” is employed in a rather loose fashion: a more complex model is one with a larger number of parameters that can be adjusted over a large range to fit the model to the observations. In this paper we will try to measure the effective number of model parameters that a given set of data can support. Because of the connection to the data, we will call this the effective complexity or Bayesian complexity. Our main purpose is to present a statistically sound quantity that embodies in a quantitative way the above notion of complexity when a model is compared to the observations.
Bayesian model comparison makes use of an Occam’s razor argument to rank models in term of their quality of fit and economy in the number of free parameters. A model with more free parameters will naturally fit the observations better, but it will also be penalized for the wasted parameter space that the larger number of parameters implies. Several studies have made use of Bayesian model comparison to assess the viability of different models in the cosmological context ev_use ; trotta_evidence ; nestsamp . In this work we show that Bayesian complexity is an ideal complement to model selection in that it allows to identify the number of effective parameters supported by the data.
We start by introducing our notation and the fundamentals of Bayesian statistics and model comparison. We then present the Bayesian complexity and illustrate its use in the context of a toy model in section III. In section IV we apply it to observations of cosmic microwave background anisotropies and we conclude in section V.
Ii Bayesian model selection and complexity
ii.1 Model comparison
We first briefly review the basic ingredients of Bayesian statistics and some relevant aspects of information theory. This serves both to introduce our notation and to remind the reader of the main points. We will use a fairly compact notation where possible and refer the reader to e.g. MacKay for the exact mathematical definitions. Specifically, for an outcome of a random variable we will write for the probability distribution function (pdf), ie the probability that takes a certain value . In the case of a multidimensional parameter space we will write as a short form of the joint probability over all components of , . The conditional probability of given is written .
The starting point of our analysis is Bayes theorem:
(1) 
Here, the quantity represents a collection of all external hypotheses and our model assumptions.
Given the data and a model with free parameters , statistical inference deals with the task of determining a posterior pdf for the free parameters of the model, . The latter is computed via Bayes theorem as
(2) 
On the right–hand side of this equation, is the probability of obtaining the observed data given the parameter value . Since the observed data are a fixed quantity, we interpret as a function of and we call it the likelihood, . The prior pdf embodies our state of knowledge about the values of the parameters of the model before we see the data. There are two conceptually different approaches to the definition of the prior: the first takes it to be the outcome of previous observations (i.e., the posterior of a previous experiment becomes the prior for the next), and is useful when updating ones knowledge about the values of the parameters from subsequent observations. For the scope of this paper, it is more appropriate to interpret the prior as the available parameter space under the model, which then is collapsed to the posterior after arrival of the data. Thus the prior constitutes an integral part of our model specification. In order to avoid cluttering the notation, we will write , with the model dependence understood whenever no confusion is likely to arise.
The expression in the denominator of (2) is a normalization constant and can be computed by integrating over the parameters,
(3) 
This corresponds to the average of the likelihood function under the prior and it is the fundamental quantity for model comparison. The quantity is called marginal likelihood (because the model parameters have been marginalized), or in recent papers in the cosmological context, the evidence. In the following we shall refer to it as to the likelihood for the model ^{1}^{1}1This is an effort to be consistent with common terminology in the statistical community and to avoid confusion with the colloquial meaning of the term “evidence”.. The posterior probability of the model is then, using Bayes theorem again,
(4) 
where is the prior for the model. The quantity in the denominator on the right–hand side is just a normalization factor depending on the data alone, which we can ignore. When comparing two models, versus , one introduces the Bayes factor , defined as the ratio of the model likelihoods
(5) 
In other words, the prior odds of the models are updated by the data through the Bayes factor. If we do not have any special reason to prefer one model over the other before we see the data, then , and the posterior odds reduce to the Bayes factor. Alternatively, one can interpret the Bayes factor as the factor by which our relative belief in the two model is modified once we have seen the data.
ii.2 A Bayesian measure of complexity
The usefulness of a Bayesian model selection approach based on the model likelihood is that it tells us whether the increased “complexity” of a model with more parameters is justified by the data. But the number of free parameters in a model is only the simplest possible notion of complexity, which we call input complexity and denote by . A much more powerful definition of model complexity was given by Spiegelhalter et al dic , who introduced the Bayesian complexity, which measures the number of model parameters that the data can constrain:
(6) 
On the right–hand side, is the KullbackLeibler (KL) divergence between the posterior and the prior, representing the information gain obtained when upgrading the prior to the posterior via Bayes theorem:
(7) 
The KL divergence measures the relative entropy between the two distributions (e.g. MacKay )^{2}^{2}2Notice that the KL divergence is not symmetric (i.e., in general ), but it does have the useful properties that (Gibbs inequality), with equality only for , and that it is invariant under general transformations of the parameter space.. In Eq. (6), is a point–estimate for the KL divergence. If all parameters are well measured, then the posterior pdf will collapse into a small region around , and thus the KL divergence will approximately be given by . By taking the difference, we compare the effective information gain to the maximum information gain we can expect under the model, . The factor of 2 is chosen so that for highly informative data, as shown below. As a point estimator for we employ the posterior mean, denoted by an overbar. Other choices are possible, and we discuss this further below.
We can use Eq. (7) and Bayes theorem to rewrite (6) as
(8) 
Defining an effective through the likelihood as (any constant factors drop out of the difference of the logarithms in Eq. (6)) we can write the effective number of parameters as
(9) 
where the mean is taken over the posterior pdf. This quantity can be computed fairly easily from a Markov chain Monte Carlo (MCMC) run, which is nowadays widely used to perform the parameter inference step of the analysis.
We thus see that the effective number of parameters of the model is not an absolute quantity, but rather a measure of the constraining power of the data as compared to the predictivity of the model, i.e. the prior. Hence depends both on the data at hand and on the prior available parameter space. In fact, it is clear that the very notion of “well measured parameter” is not absolute, but it depends on what our expectations are, i.e. on the prior. For example, consider a measurement of , the total energy density of the Universe, expressed in units of the critical density. The current posterior uncertainty around is about . Whether this means that we have “measured” the Universe to be flat (i.e. ) or not depends on the prediction of the model we consider. If we take a generic prior in the range , then we conclude that current data have measured the Universe to be flat with moderate odds (a precise analysis gives odds of in favor of the flat model, see trotta_evidence ). On the contrary, in the framework of e.g. landscape theories, the prior range of the model is much narrower, say
Eq. (9) is conceptually related to the approach used for example in cmpcb to analyze how many dark energy parameters are measured by the data. In that work the decrease of the value within the confidence regions encompassing 68% and 95% of the posterior was measured and compared against the theoretical decrease of a multivariate Gaussian distribution with a given number of degrees of freedom (see e.g. chapter 15.6 of numrep ). This in turn allowed to deduce the effective number of degrees of freedom of the .
We now turn to an explicit demonstration of the power and usefulness of coupling the model likelihood with the Bayesian model complexity in model selection questions, by analyzing in some detail linear toy models.
Iii Linear models
Before applying the Bayesian complexity and model likelihood to current cosmological data, we compute them explicitly for a linear model and we illustrate their use in a toy example involving fitting a polynomial of unknown degree. We show below how the Bayesian complexity tells us how many parameters the data could in principle constrain given the prior expectations under the model.
iii.1 Specification of the likelihood function
Let us consider the following linear model
(10) 
where the dependent variable is a dimensional vector of observations, is a vector of dimension of unknown regression coefficients and is a matrix of known constants that specify the relation between the input variables and the dependent variables ^{3}^{3}3In the special case of observations fitted with a linear model , the matrix is given by the basis functions evaluated at the locations of the observations, .. Furthermore, is a dimensional vector of random variables with zero mean (the noise). If we assume that follows a multivariate Gaussian distribution with uncorrelated covariance matrix , then the likelihood function takes the form
(11) 
where we have defined and . This can be cast in the form
(12) 
with the likelihood Fisher matrix given by
(13) 
and a normalization constant
(14) 
Here denotes the parameter value that maximizes the likelihood, given by
(15) 
As a shortcut, we will write
(16) 
where .
iii.2 Model likelihood and complexity
In this section we gain some intuitive feeling about the functional dependence of the model likelihood and complexity on the prior and posterior for the simple case of linear models outlined above. The results are then applied in section III.4 to an explicit toy model, showing the model likelihood and complexity in action.
Assuming as a prior pdf a multinormal Gaussian distribution with zero mean and Fisher information matrix (we remind the reader that the Fisher information matrix is the inverse of the covariance matrix), i.e.
(17) 
where denotes the determinant of the matrix , the model likelihood (3) and model complexity (9) of the linear model above are given by Eqs. (25) and (31), respectively (see Appendix A).
Let us now consider the explicit illustration of a model with parameters, and . Without losing generality, we can always choose the units so that the prior Fisher matrix is the unity matrix, i.e.
(18) 
This choice of units is natural since it is the prior that sets the scale of the problem. The likelihood Fisher matrix being a symmetric and positive matrix, it is characterized by real numbers, which we choose to be its eigenvalues , and the elements of the orthogonal matrix that diagonalizes it (corresponding to rotation angles). Here the represent the standard deviations of the likelihood covariance matrix along its eigendirections , expressed in units of the prior width. With we have that the likelihood Fisher matrix is given by and thus Eq. (31) gives
(19)  
The complexity only depends on how well we have measured the eigendirections of the likelihood function with respect to the prior. Every well–measured direction, (i.e., one for which ) counts for one parameter in , while directions whose posterior is dominated by the prior, , do not count towards the effective complexity. This automatically takes into account strong degeneracies among parameters. Notice also that once an eigendirection is well measured (i.e. in the limit ), then the prior width does not matter anymore.
In contrast, the model likelihood (25) is given by (assuming for simplicity that the mean of the likelihood corresponds to the prior mean, i.e. )
(20) 
Finally, we remark that an important ingredient of the Bayesian complexity is the point estimator for the KL divergence. Here we adopt the posterior mean as an estimate, but other simple alternatives are certainly possible, for instance the posterior peak (or mode), or the posterior median. The choice of an optimal estimator is still matter of research (see e.g. section IV and the comments at the end of Ref. dic ). The important aspect is that the posterior pdf should be summarized by only one number, namely the value plugged into the KL estimator. This is obviously going to be a very bad description for highly complex pdf’s, exhibiting for instance long, banana–shaped degeneracies. No single number can be expected to summarize accurately such a pdf. On the other hand, for fairly Gaussian pdf’s all the different estimators (mean, median and peak) reduce to the same quantity. This clearly calls for using normal directions in parameter space kosowsky , which make the posterior as Gaussian as possible, a procedure that it would be wise to follow whenever possible for many other good reasons (e.g., better and faster MCMC convergence).
iii.3 Effective complexity as a data diagnostics
We now turn to the question of how we can use the model likelihood and complexity together as a tool for model selection. We shall see that the Bayesian complexity provides a way to assess the constraining power of the data with respect to the model at hand and to break the degeneracy among models with approximately equal model likelihood.
Let us consider two models and with different numbers of parameters, (but in general the two models need not to be nested). If the additional parameters of model are required by the data, then the likelihood of model will be larger and will have larger posterior odds, thus it should be ranked higher in our preference than model . However, even if the extra parameters of model are not strictly necessary, they can lead to over–fitting of the data and compensate the Occam’s penalty factor in (20) sufficiently to lead to a comparable marginal likelihood for both models. The effective complexity provides a way to break the degeneracy between the quality–of–fit term () and the Occam’s razor factor in the marginal likelihood, and enables us to establish whether the data is good enough to support the introduction of the additional parameters of model .
To summarize, we are confronted with the following scenarios:

: model is clearly favored over model and the increased number of parameters is justified by the data.

and : the quality of the data is sufficient to measure the additional parameters of the more complicated model, but they do not improve its likelihood by much. We should prefer model , with less parameters.

and : both models have a comparable likelihood and the effective number of parameters is about the same. In this case the data is not good enough to measure the additional parameters of the more complicated model and we cannot draw any conclusions as to whether the additional complexity is warranted.
We illustrate these cases by computing the model likelihood and effective complexity of a toy model in the next section.
iii.4 An illustrative example
As a specific example of the linear model described in section III.1, we consider the classic problem of fitting data drawn from a polynomial of unknown degree. The models that we test against the data are a collection of polynomials of increasing order, with input complexity , where is the order of the polynomial. The question is then whether our model selection can correctly recover the order of the polynomial from which the data are actually drawn.
The data covariance matrix is taken to be diagonal and with a common standard deviation for all points, , while the prior over the polynomial coefficients is a multivariate Gaussian with covariance matrix given by the identity matrix. For definitiveness, we will take the underlying model from which the data are generated to have parameters. First we draw data points with noise . We plot in Figure 1 the resulting model likelihood and effective complexity as a function of the input model complexity, . The likelihood of the models increases rapidly until and then decreases slowly, signaling that parameters are not justified. (We plot the logarithm of the logarithm of the model likelihood as the models with parameters are highly disfavored and would otherwise not fit onto the figure – an example of case (1) of the list in the previous section). The effective complexity on the other hand continues to grow until , at which point the data is unable to constrain more complex models and becomes constant. We conclude that the model with is the one preferred by data and that additional parameters are not needed, although the data could have supported them. This is case (2) in the list of the previous section.
Next we decrease the number of data points to , in which case we obviously cannot recover more than four parameters. It comes as no surprise that the model likelihood stops increasing at , see Figure 2. But the effective complexity also flattens at , which means that the data cannot deal with more than four parameters, irrespective of the underlying model! In this case, corresponding to point (3) of the list in the previous section, we conclude that the available data do not support more than 4 effective parameters. On the other hand, we recognize that the flattening of the model likelihood at does not necessarily mean that the underlying model has four parameters. We thus must hold our judgment until better data become available.
As an alternative to decreasing the number of data points, we can achieve a similar data degradation effect by keeping data points but by increasing the noise to . We obtain a result similar to the previous case, which is plotted in Figure 3.
We conclude by emphasizing once more that in general the outcome of model selection based on assessing the model likelihood and effective complexity depends on the interplay of two factors. The first is the predictive power of the model, as encoded in the prior. The second is the constraining power of the data.
Iv How many parameters does the CMB need?
We now apply the above tools to the question of how many cosmological parameters are necessary to describe current cosmic microwave background (CMB) anisotropy measurements. We make use of the following CMB data: WMAP, ACBAR, CBI, VSA and Boomerang 2003. To provide an additional regularization (especially in view of including spatial curvature) we also use the HST limits on the Hubble parameter, km/s/Mpc. We should note that this strongly increases the power of the CMB data. We use both the firstyear WMAP alone (WMAP1) as well as the data for the first three years (WMAP3).
For each set of cosmological parameters we create a converged MCMC chain using the publicly available code cosmoMC cosmomc . We then compute the Bayesian complexity from the chain through Eq. (9). The model likelihood is evaluated with the SavageDickey method (see trotta_evidence and references therein): For a model that is nested within a larger model by fixing one parameter, , to a value , the Bayes factor between the two models is given by the posterior of the larger model at (normalized and marginalized over all other parameters) divided by the prior at this point. In this way it is possible to derive all the model likelihoods in a hierarchy, starting from the most complex model (which is assigned an arbitrary model likelihood, in our case 1). Since errors tend to accumulate through the intermediate steps necessary to reach the simpler models, and as a cross–check, we additionally computed the model likelihoods with nested sampling nestsamp for the WMAP1 data. Within the errorbars, we did not find any appreciable discrepancy between the two methods.
In this analysis we use four basic cosmological parameters, namely
(21) 
where () is the baryon (matter) density relative to the critical energy density, km/s/Mpc is the fudge factor, is the ratio of the sound horizon to the last scattering angular diameter distance and , with the power spectrum of adiabatic density fluctuations at a scale Mpc. We then add three more parameters in various combinations, to study whether they are necessary and supported by the observations. The additional parameters are: The reionization optical depth , the scalar spectral index and the spatial curvature (parameterized by its contribution to the Hubble equation, ).
The scalar spectral index is either fixed to (the case of a scale invariant power spectrum of initial fluctuations), or else chosen with a Gaussian prior, . We find that . We choose the prior of the reionization optical depth to reflect our current lack of understanding of how reionization proceeds. We choose it flat within and then add an exponential falloff:
(22) 
The models without this parameter have , and . The curvature contribution is either set to zero, , if we do not include the parameter, or else is used with a flat prior . The value of the prior for a flat universe is . We find that adding curvature as a free parameter when using the WMAP 3yr data leads to a nonGaussian posterior for which the mean as a pointestimate for the KL divergence in Eq. (6) is not representative. We opt here for a slightly modified estimator, given by the average of the evaluated at the mean and the mode of the posterior. For a Gaussian posterior this reduces to the mean pointestimator but it appears to be somewhat more stable.
We quote our results in Table 1 (WMAP3) and Table 2 (WMAP1), while Figure 4 gives a graphical representation. The model likelihoods are quoted relative to the model with the most parameters. We find that using WMAP1 we can measure all the parameters for the models with four and five parameters. For , the effective complexity increases more slowly than the number of input parameters, but we can still measure at least six parameters with CMB+HST. With WMAP3 we can measure all six parameters of the model. We conclude that the new WMAP3 data augmented by the HST determination of can measure all seven parameters considered in this analysis.
Taking into account the model likelihood values, we find that the models , and (to a lesser extent) as well as are strongly disfavoured. This shows that is preferred by current data, in agreement with the result of Ref. trotta_evidence . In general, adding in a non–zero spatial curvature leads to a well measurable decrease in the model probability that, together with the increase in effective model complexity, reinforces our belief that can be safely neglected for the time being. Of course this result is partially a reflection of our choice of prior on . However, it is important to remember that had we halved the range of this prior, the likelihoods for models with non–zero curvature would have only doubled. This would not change the results significantly. Alternatively, an inflation–motivated prior of the type would render the parameter unmeasured and irrelevant. In this case adding it would not change the effective complexity or the model likelihood at all.
We also find that the basic set must be augmented by either or . The inclusion of both parameters at the same time was optional with the first year WMAP data only, but using WMAP3 we find that has a significantly higher model likelihood than all other models investigated here. Also, where with WMAP1 we only gained half an effective parameter when going from or to , we now gain one full parameter. Thus we can now measure both parameters at the same time.
Overall, we conclude that was a good and sufficient base model until a few months ago. Now should be used. A wider prior on would have only a minimal impact on the complexity of models including a tilt, since seems to be rather well–measured when considered alone (ie, not in combination with ). Inclusion of a non–vanishing curvature is discouraged by Bayesian model comparison. We find that a model with 6 parameters is sufficient to explain the current CMB data, even though all seven effective parameters can be constrained now. This analysis demonstrates that the Bayesian complexity estimator (9) works with real–world data and gives useful additional information for model comparison.
V Conclusions
In this work we introduced the Bayesian complexity as a measure of the effective number of parameters in a model. We discussed extensively its properties and its usefulness in the context of linear models, where it can be computed analytically. We showed that it corresponds to the number of parameters for which the width of the posterior probability distribution is significantly narrower than the width of the prior probability distribution. These parameters can be considered to have been well measured by the data given our prior assumptions in the model.
We also showed that for linear models the Bayesian complexity probes the trace of the posterior covariance matrix, while the model likelihood is sensitive to the determinant. We argued that the Bayesian complexity allows to test for cases where the data is not informative enough for the model likelihood to be a reliable indicator of model performance, and provided an explicit example.
Finally, we applied the combination of model likelihood and Bayesian complexity to the question of how many (and which) parameters are measured by current CMB data, complemented by the HST limit on the Hubble parameter. We limited ourselves to the family of CDM models with a power–law spectrum of primordial perturbations. We demonstrated that – in addition to the energy density in baryons and matter, the CMB peak location parameter and the amplitude of the initial perturbations – we need to consider now both the reionisation optical depth and the scalar spectral index. Non–flat models are disfavoured. The effective complexity shows that the CMB data can measure all seven parameters in this scheme.
As the Bayesian complexity is very easy to compute from a MCMC chain, we hope that it will be used routinely in future data analyses in conjunction with the model likelihood for model building assessment. It will help to determine if the data is informative enough to measure the parameters under consideration. Further work is needed to study the performance of the Bayesian complexity in situations with a strongly non–Gaussian posterior.
Acknowledgements.
We thank Andrew Liddle, Pia Mukherjee and Glenn Starkman for useful discussions. M.K. is supported by the Swiss NSF. R.T. is supported by the Royal Astronomical Society through the Sir Norman Lockyer Fellowship. D.P. is supported by PPARC. The calculations were performed on the “myri” cluster of the University of Geneva and the UK national cosmology supercomputer (COSMOS) in Cambridge.Appendix A Model likelihood and complexity in the linear case
Here we compute first the model likelihood (3) for the linear model (12), using the Gaussian prior (17). An analogous computation can be found in dic . Returning to Bayes theorem (1), the posterior pdf is given by a multinormal Gaussian with Fisher information matrix
(23) 
and mean given by
(24) 
The model likelihood (3) evaluates to
(25)  
This can be easily interpreted by looking at its three components: the quality of fit of the model is encoded in , which represents the best–fit likelihood. Thus a model that fits the data better will be favored by this term. The term involving the determinants of and is a volume factor (the so called Occam’s factor). As , it penalizes models with a large volume of wasted parameter space, i.e. those for which the parameter space volume that survives after arrival of the data is much smaller than the initially available parameter space under the model prior, . Finally, the exponential term suppresses the likelihood of models for which the parameters values that maximize the likelihood, , differ appreciably from the expectation value under the posterior, . Therefore when we consider a model with an increased number of parameters we see that its model likelihood will be larger only if the quality–of–fit increases enough to offset the penalizing effect of the Occam’s factor.
Let us now turn to the computation of the Bayesian complexity, Eq. (9). Using the posterior mean (denoted by an overbar) as a point estimator for the effective chi–square, we obtain from (16)
(26) 
The expectation value of the under the posterior is given by
(27) 
We concentrate on the second term and write . The total term in the expectation value then becomes . The first term of this expression is just the posterior covariance matrix . The last term combines with to . The cross terms vanish since .
All taken together, we obtain for the complexity
(28)  
(29) 
Using the relation (23) we can rewrite the complexity as
(30)  
(31) 
Thus while the model likelihood depends on the determinant of the Fisher matrices, the complexity depends on their trace. Another important point worth highlighting is that for linear models the complexity does not depend on the degree of overlap between the prior and the posterior, nor on the noise realization (as long as the noise covariance matrix is known).
References
 (1) A.R. Liddle, Mon. Not. Roy. Astron. Soc. 351, L49L53 (2004); T.D. Saini, J. Weller and S.L. Bridle, Mon. Not. Roy. Astron. Soc. 348, 603 (2004); B.A. Bassett, P.S. Corasaniti and M. Kunz, Astrophys.J. 617, L1 (2004); G. Lazarides, R.R. de Austri and R. Trotta, Phys. Rev. D 70, 123527 (2004); A. Niarchou, A.H. Jaffe, L. Pogosian, Phys.Rev. D69 (2004) 063515; P. Mukherjee et al, astroph/0512484 (2005); M. Beltran et al, Phys. Rev. D 71, 063532 (2005); M. Kunz et al, Phys. Rev. D 73, 023511 (2006).
 (2) R. Trotta, astroph/0504022 (2005).
 (3) P. Mukherjee, D. Parkinson and A.R. Liddle, astroph/0508461 (2005).
 (4) D.J.C. MacKay, Information Theory, Inference, and Learning Algorithms, Cambridge University Press (2003).
 (5) C.E. Shannon, Bell Sys. Tech. J. 27, 379 and 623 (1948).
 (6) D.J. Spiegelhalter et al, J.R. Statist. Soc. B 64, 583 (2002).
 (7) D. Feldmann, Information Theory, Excess Entropy and Computational Mechanics, available from http://hornacek.coa.edu/dave/
 (8) E.T. Jaynes and G.L. Bretthorst (ed.), Probability Theory: The Logic of Science, Cambridge University Press (2003).
 (9) W.H. Press et al, Numerical Recipes in C (second edition), Cambridge University Press (1992).
 (10) P.S. Corasaniti et al, Phys. Rev. D 70, 083006 (2004).
 (11) A. Lewis and S. Bridle, Phys.Rev. D 66, 103511 (2002). The code is available from http://cosmologist.info
 (12) A. Kosowsky, M. Milosavljevic and R. Jimenez, Phys. Rev. D 66, 63007 (2002); M. Chu, M. Kaplinghat and L. Knox, Astrophys. J. 596, 725 (2003).