Accounting for Correlations When Fitting Extra Cosmological Parameters

In general it is not straightforward to infer the correlation between parameters A and B directly from Markov Chain Monte Carlo (MCMC) fits to the ΛCDM+A and ΛCDM+B models that have already been performed by experimental collaborations like Planck. To make progress, we assume that the data likelihood is Gaussian, with a covariance matrix ${\boldsymbol{ \mathcal C }}$ that does not depend on the cosmological parameters, and that the posterior distribution for the cosmological parameters is also Gaussian.

The assumption of Gaussianity of the likelihood is made widely across different cosmological measurements, including CMB power spectra (e.g., Planck Collaboration XI 2016; Louis et al. 2017; Henning et al. 2018), weak lensing shear (e.g., Krause et al. 2017; Wright et al. 2018; Hikage et al. 2019), galaxy clustering (e.g., Percival et al. 2014; Alam et al. 2017), supernovae distance moduli (e.g., Scolnic et al. 2018; Abbott et al. 2019), and others. The likelihood can almost always be made more Gaussian through compression of the data, often with negligible loss of cosmological information (e.g., combining CMB power spectra over a range of multipoles into bins). Neglecting the cosmological parameter dependence in the covariance (e.g., assuming a fixed fiducial model and set of parameters for computing the cosmic variance contribution to the errors) has also been demonstrated to be a suitable approximation for Planck, and other current experiments (e.g., Hamimeche & Lewis 2008; Krause et al. 2017). The fiducial model as well as the covariance are often estimated from some iterative process of fitting the actual data.

For the Planck data, cosmological parameter posterior distributions can be well approximated as Gaussian for ΛCDM, as well as for many one- and two-parameter extensions, although there are also cases (e.g., involving curvature or varying the dark energy equation of state), where Planck alone does not provide Gaussian posteriors. Other experiments only constrain a portion of the ΛCDM parameter space sufficiently well to produce Gaussian constraints in one or two parameters. Many parameters are also physically required to be positive, truncating the available parameter space and causing departure from Gaussianity when the data do not constrain the parameter to significantly differ from zero. This is the case for current cosmological constraints on the neutrino mass, for example (see Figure 30 of Planck Collaboration XIII 2016). We return to the handling of non-Gaussian cases in Section 4.3.

In certain cases, the Bayesian posterior parameter distribution and the distribution of the ML parameter values from many realizations of the data are approximately equal. Specifically, ${\mathbb{P}}({\boldsymbol{\theta }}| {\boldsymbol{d}})$, the Bayesian posterior distribution of parameters ${\boldsymbol{\theta }}$ sampled by MCMC given the experimental data, can well approximate ${\mathbb{P}}({{\boldsymbol{\theta }}}^{\mathrm{ML}}({{\boldsymbol{d}}}^{\mathrm{sim}})| {{\boldsymbol{\theta }}}^{\mathrm{fid}})$. The latter is the distribution of frequentist ML parameter estimation based on realizations of ${{\boldsymbol{d}}}^{\mathrm{sim}}$ from a fiducial model ${{\boldsymbol{\theta }}}^{\mathrm{fid}}$. The choice of ${{\boldsymbol{\theta }}}^{\mathrm{fid}}$ is usually physically motivated and is based on fits from actual data. For a detailed discussion, see e.g., Chapter 4 and Appendix B of Gelman et al. (2013).

We show in Appendix A that this correspondence is mathematically exact, if we: (1) assume the Gaussianity of both the data likelihood and the posterior distribution of parameters; (2) impose uninformative parameter priors that are far less constraining than the likelihood; (3) in the MCMC computation, replace the experimental data ${\boldsymbol{d}}$ with ${\boldsymbol{\mu }}({{\boldsymbol{\theta }}}^{\mathrm{fid}})$, the theory prediction of the fiducial model. In other words, ${\mathbb{P}}({\boldsymbol{\theta }}| {\boldsymbol{d}}={\boldsymbol{\mu }}({{\boldsymbol{\theta }}}^{\mathrm{fid}}))$ from MCMC and ${\mathbb{P}}({{\boldsymbol{\theta }}}^{\mathrm{ML}}({{\boldsymbol{d}}}^{\mathrm{sim}})| {{\boldsymbol{\theta }}}^{\mathrm{fid}})$ from simulations are the same mathematically. In this work, the fiducial model is from fitting the Planck TT data. We will show in Section 3.2 that the exact choice of ${{\boldsymbol{\theta }}}^{\mathrm{fid}}$ is not very important.

This correspondence allows us to estimate the correlation between extension parameters A and B using simulated data (i.e., frequency sampling of the likelihood) in the following steps:

• 1.  
  Generate many simulated data sets (in our case simulated Planck-like CMB power spectra) drawn from the likelihood in the form of a Gaussian distribution ${ \mathcal N }({\boldsymbol{\mu }}({{\boldsymbol{\theta }}}^{\mathrm{fid}}),{\boldsymbol{ \mathcal C }})$ for some choice of fiducial ΛCDM model parameters, ${{\boldsymbol{\theta }}}^{\mathrm{fid}}$, and using the covariance matrix ${\boldsymbol{ \mathcal C }}$ provided by the experiment collaboration.
• 2.  
  Calculate the maximum-likelihood parameters, ${{\boldsymbol{\theta }}}^{\mathrm{ML}}$, for the ΛCDM+A and ΛCDM+B models for each simulated data set, and
• 3.  
  Estimate the covariance between A in the ΛCDM+A fit and B in the ΛCDM+B fit using the sample covariance from the simulations.

Alternatively, in the special case of Gaussian parameter posteriors, we can run MCMC chains to estimate the correlation between extension parameters directly. In Appendix B, we show that all the information on the correlation between extension parameters A and B can be estimated from three sets of chains: ΛCDM+A, ΛCDM+B, and ΛCDM+A+B. Specifically, we found that the correlation between A and B when they are fitted separately is equal to minus one times the correlation between A and B when they are fitted together. We refer to this property as correlation equivalence.

A complication of estimating correlation between extension parameters by performing MCMC on the real data in this way is that we clearly do not have control over the "underlying" cosmological model in the way we do when generating simulations. It is also possible that the real data have imperfections or systematic biases that are not correctly accounted for in the likelihood. We therefore instead performed MCMC computations replacing the real data vector with a theory prediction computed using the same ${{\boldsymbol{\theta }}}^{\mathrm{fid}}$ as for the simulations.

Here we provide a worked example using Planck data for both the simulation and MCMC approaches discussed above. At the time of writing, the Planck 2018 likelihood code is not yet available. Therefore, we perform a simple three-parameter test using the Planck ℓ ≥ 30 2015 TT data from the plik_lite likelihood¹ (as described in Planck Collaboration XI 2016). This simplified likelihood includes only CMB information, marginalizing over foreground template amplitudes and other nuisance parameters.

We assume the fiducial model to be the best-fit ΛCDM model of ℓ ≥ 30 Planck TT plik_lite data with the optical depth fixed at τ = 0.07, the Planck calibration parameter calPlanck equal to 1 and other nuisance parameters marginalized. The ℓ < 30 likelihood is pixel based rather than power spectrum based. For simplicity we only include the power spectrum based likelihood of the ℓ ≥ 30 data. Moreover, because the ℓ ≥ 30 TT spectrum only well constrains the parameter combination A_se^−2τ, τ is fixed to break the strong degeneracy with A_s. As for the calibration parameter calPlanck, a Gaussian prior with mean equal to 1 and standard deviation 0.0025 was originally imposed on it. Its variation has minimal impact on other parameters. For simplicity, calPlanck is fixed at 1.

For the standard ΛCDM parameters, we use

$$\begin{eqnarray*}&&\begin{array}{l}\{{{\rm{\Omega }}}_{b}{h}^{2},{{\rm{\Omega }}}_{c}{h}^{2},100{\theta }_{{MC}},\tau ,\mathrm{log}{A}_{s},{n}_{s}\}{}^{\mathrm{fid}}\\ \ =\ \{0.02216,0.1211,1.0407,0.07,3.0777,0.9601\}.\end{array}\end{eqnarray*}$$

In addition, the three extension parameters of interest and their fiducial values for this example are ${\{{A}_{L},{n}_{\mathrm{run}},{Y}_{P}\}}^{\mathrm{fid}}\,=\{1,0,0.2453\}$. As mentioned in Section 1, A_L is a phenomenological parameter that artificially scales the lensing power spectrum. It is worth investigating as Planck has a curious preference for A_L > 1. Running of the spectral index, n_run ≡ dn_s/d ln k, and the primordial Helium mass fraction Y_P are an interesting pair of extension parameters to test because of their high correlation (∼0.9), that is, they produce similar changes in the power spectrum damping tail. Even a small deviation from their degeneracy direction in their expected 2D distribution would be noticeable. In ΛCDM, Y_P is calculated from Ω_bh² and the CMB temperature through big bang nucleosynthesis (BBN) predictions (Planck Collaboration XIII 2016). In ΛCDM+Y_P, Y_P is independent of BBN and decouples from Ω_b h². We impose flat priors on all model parameters with the default CosmoMC (Lewis & Bridle 2002) bounds.

The fiducial parameters from fitting the Planck data have uncertainties due to cosmic variance and experimental noise. A slightly different set of fiducial parameters, which are still consistent with the measured power spectrum, may produce different results in the covariance of the extension parameters as well as the significance of their deviations from the fiducial ΛCDM values. We found this effect to be small for the Planck 2015 data. See Section 3.2 for details.

We run simulations to sample the distribution of the extension parameters. We draw 1000 binned power spectrum samples with mean equal to the binned power spectrum of the fiducial model and covariance equal to the plik_lite CMB band power covariance. We use the same binning scheme as Planck 2015 to convert power spectra computed by CAMB² to band powers. We replace the data spectrum in the plik_lite likelihood with our samples, thus forming the simulated likelihoods. Next, using the ML finding algorithm in CosmoMC, setting τ = 0.07 and the Planck calibration parameter calPlanck = 1, we maximize the simulated likelihoods to obtain the best-fit parameters of a specific model for each realization. The models we explore are: the standard ΛCDM, ΛCDM+A_L, ΛCDM+n_run and ΛCDM+Y_P.

To quantify the overall shift between the values of the extension parameters estimated from the actual data and their fiducial ΛCDM values, we calculate the χ² and its probability to exceed (PTE) for a χ² distribution with degrees of freedom equal to the number of extension parameters:

$$\begin{eqnarray}&&{\chi }^{2}={\left({{\boldsymbol{\theta }}}_{\mathrm{dat}}^{\mathrm{ML}}-{{\boldsymbol{\theta }}}^{\mathrm{fid}}\right)}^{T}{{\boldsymbol{\Sigma }}}_{\mathrm{ex}}^{-1}({{\boldsymbol{\theta }}}_{\mathrm{dat}}^{\mathrm{ML}}-{{\boldsymbol{\theta }}}^{\mathrm{fid}}),\end{eqnarray} \tag{ 1 }$$

where ${{\boldsymbol{\Sigma }}}_{\mathrm{ex}}$ is the covariance matrix for the extension parameters only, ${{\boldsymbol{\theta }}}^{\mathrm{fid}}$ is an array of three elements equal to ${\{{A}_{L},{n}_{\mathrm{run}},{Y}_{P}\}}^{\mathrm{fid}}$, and ${{\boldsymbol{\theta }}}_{\mathrm{dat}}^{\mathrm{ML}}=(1.1245,0.00773,0.2306)$ is an array whose values are obtained from fitting models ΛCDM+A_L, ΛCDM+n_run, and ΛCDM+Y_P, respectively, on the actual Planck plik_lite TT data with τ fixed.

To verify that our priors on the parameters can be considered as uninformative and thus do not impact the degree of freedom for the χ², we compute the χ² values and their PTE assuming three degrees of freedom using the best-fit parameter values from all simulations. Then we compare the distribution of PTE values with the expected one with three degrees of freedom and find that they are consistent.

As mentioned in Section 2.2, we can make use of correlation equivalence to estimate the expected covariance ${{\boldsymbol{\Sigma }}}_{\mathrm{ex}}$ between extension parameters fitted separately by running MCMC chains on the data likelihood, with the experimental power spectrum replaced by the fiducial one. Again we set τ = 0.07 and the Planck calibration parameter calPlanck = 1.

The diagonal elements in ${{\boldsymbol{\Sigma }}}_{\mathrm{ex}}$ are estimated from the variances of the specific extension parameters from running the MCMC chains with the modified Planck plik_lite TT likelihood, on the model ΛCDM+A_L, ΛCDM+n_run and ΛCDM+Y_P.

For the off-diagonal elements, first we obtain the correlation between two extension parameters (again, denoting them generically as A and B) from the MCMC runs that vary both of them in the ΛCDM+A+B model with the same data set. Then we calculate the covariance between A and B, given their variances as described in the last paragraph:

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{Cov}{\left(A,B\right)}_{\mathrm{fitted}\mathrm{separately}}=-\mathrm{corr}{\left(A,B\right)}_{{\rm{\Lambda }}\mathrm{CDM}+A+B}\\ \quad \ \times \ {\left(\mathrm{Var}{\left(A\right)}_{{\rm{\Lambda }}\mathrm{CDM}+A}\times \mathrm{Var}{\left(B\right)}_{{\rm{\Lambda }}\mathrm{CDM}+B}\right)}^{\tfrac{1}{2}}.\end{array}\end{eqnarray} \tag{ 2 }$$

This way, we are able to estimate all of the elements in ${{\boldsymbol{\Sigma }}}_{\mathrm{ex}}$ without running simulations at all, and can calculate the χ² defined in Equation (1) and its PTE to quantify the shifts of the extension parameters from their fiducial values.

In Table 1, we show the χ² of the difference, as defined in Equation (1), their corresponding PTE values and the level of consistency with the assumed model in terms of the number of σ. Our results imply no significant discrepancy (0.9σ) between the values of the three additional parameters estimated from the actual data in one-parameter extensions and their fiducial ΛCDM values. This is consistent with expectations if the standard ΛCDM is the correct model.

Taking a closer look at elements in the estimated covariance, we show in Table 2 the estimated uncertainties of one-parameter extensions to the base ΛCDM model, along with the estimated parameter correlations from the two methods. These results are consistent up to at most 5%, in the Y_P uncertainty. The differences between using the two methods are likely due to numerical noise and have only a minimal impact on our main conclusion.

We also visualize the shifts of the ML extension parameters from their fiducial values compared with their expected distributions in Figure 2. Notice in the 2D contour plots, how the ML points lie along the correlation/degeneracy direction of the estimated distributions, as expected if the base ΛCDM is the true model. Looking only at single extensions, the discrepancy between the ML value of A_L inferred by the plik_lite TT data and its fiducial value is only significant at 1.8σ (from the MCMC method) or 1.9σ (from simulations), while the deviations of n_run and Y_P from their fiducial values are only significant at 0.8σ and 0.7σ, respectively, from both running MCMC and simulations.

Zoom In Zoom Out Reset image size

We do not attempt to explain or deemphasize the tension in A_L. Rather, we note that there is no extra sign of discrepancy when we take into account its correlation with n_run or Y_P. The significance of the χ² of difference for the three parameters is not only reduced by the very small deviations of Y_P and n_run, but is also reduced because the Planck plik_lite TT data prefers the shifts in the parameters along their expected degeneracy directions. Considering the A_L–n_run or the A_L–Y_P pair for example, their correlations are nonnegligible: 0.31 and −0.46, respectively. This information is not contained in the single extension tests. Table 1 and Figure 2 together show that the correlation of A_L with n_run and Y_P are as expected if ΛCDM is the true model, even though A_L deviates moderately with its fiducial value. Had we found inconsistency between the ML parameters from the data and the expected ΛCDM joint distribution, it could indicate the presence of systematic error or unknown physical effects producing spurious parameter correlations.

In addition, numerical differences between the simulated contours and those of MCMC are small in this figure, again providing confidence to both methods.

As mentioned in Section 2.4, the fiducial parameters are from fitting the actual plik_lite data and therefore have uncertainties. So the standard ΛCDM model we define is not just one point in the parameter space, but an ensemble of points whose shape is described by the parameter posterior distribution inferred from the data.

To assess the impact of different fiducial values on our results, we randomly draw 1000 sets of parameters from the MCMC chains fitted to the plik_lite data, with the base ΛCDM model, τ = 0.07 and calPlanck = 1. For each set of parameters, we approximate the parameter covariance matrix by computing the C_ℓ derivatives and Fisher matrices (see Equation (14) in Appendix B). Then, we use correlation equivalence to calculate the extension parameter covariance ${{\boldsymbol{\Sigma }}}_{\mathrm{ex}}$. We find that varying parameters causes less than 3% scatter in the matrix elements.

We also plug ${{\boldsymbol{\Sigma }}}_{\mathrm{ex}}$ into Equation (1) and find that the resulting χ² ranges from 4.2 to 5.1 and the PTE from 0.17 to 0.24, corresponding to consistency within 0.7–1σ. These results show that the uncertainties in the fiducial model have very minimal impact on results from our consistency test. The stability of the results reflects the constraining power of the Planck data on the model parameters in ΛCDM.

In the ideal case where all extension parameters of interest are Gaussian, we can estimate their expected distribution from MCMC chains. We denote the individual extension parameter as θ_ex,i with i runs over 1 to n_ex, the total number of extension parameters. The suggested recipe is as follows:

• 1.  
  Define a fiducial model (e.g., a fit from existing data) and calculate its theory prediction.
• 2.  
  In the data likelihood of interest, substitute the experimental observables, e.g., power spectra in the CMB case, by the fiducial prediction.
• 3.  
  Explore parameter space around the fiducial point by running MCMC chains on the modified data likelihood(s), fitting models ΛCDM+θ_ex,i and ΛCDM+θ_ex,i + ${\theta }_{\mathrm{ex},{\rm{j}}\ne i}$.
• 4.  
  Calculate the variance of θ_ex,i from the one-parameter extension, and the correlation between θ_ex,i and ${\theta }_{\mathrm{ex},j\ne i}$ from the two-parameter extension.
• 5.  
  Using results from previous step, construct a covariance matrix ${{\boldsymbol{\Sigma }}}_{\mathrm{ex}}$ for ${{\boldsymbol{\theta }}}_{\mathrm{ex}}$, with the signs of all the correlations flipped.
• 6.  
  Calculate χ² defined in Equation (1) and its PTE to quantify deviation of experimental values to fiducial ones. Additionally, one can plot confidence ellipses such as in Figure 2 to visualize the deviations.
• 7.  
  Check validity of the input fiducial model, for example by using a Fisher Forecast (e.g., Heavens 2009) to approximate parameter covariance and estimating the shifts in the χ² when varying fiducial parameter values.

As an example, on a computer cluster with 12-core 2.5 GHz processors,³ one MCMC run with 8 parallel chains usually take ≲24 hr to converge for the 2015 plik_lite likelihood. The CPU time for running one MCMC job is at the order of 100 hr and the total computing time is ∼5 × 10⁴ hr, for running all one-parameter and two-parameter extension fits for e.g., n_ex = 11, which is the number of extra parameters fitted in the publicly available 2015 Planck chains.

We can use simulations as an alternative method to estimate the expected distribution of additional parameters around the fiducial point. As described in Sections 2.1 and 2.5, one can perform the following procedure:

• 1.  
  Generate simulated observables of the fiducial model using data covariance.
• 2.  
  For each simulation, estimate the best-fit ΛCDM+${{\boldsymbol{\theta }}}_{i}^{\mathrm{ex}}$ model for each extension parameter.
• 3.  
  Calculate the covariance of extension parameter ${{\boldsymbol{\Sigma }}}_{\mathrm{ex}}$, the χ² of difference defined in Equation (1) and its PTE to quantify the significance of difference. Confidence ellipses can also be plotted.

For a rough estimate of computing time for the simulation method, using CosmoMC and the same computer cluster with 12 cores per CPU, we find that the running time of the best-fit-finding algorithm is approximately 1 hr. For n_ex = 11, n_sim = 2000 and running 4 parallel best-fit-finding jobs for each simulation (to reduce numerical noise and avoid obtaining results from local minima), the total computing time is 9 × 10⁴ hr.

To discuss the case of non-Gaussianity, we first need to clarify from what exactly non-Gaussianity arises. Recall that throughout this paper, we assume a Gaussian data likelihood and priors that are uninformative. If the data constrains the parameters well enough, changes in the model predictions can be treated as only linearly dependent on the parameters. This means the Taylor expansion of the log likelihood around the maximum point is significant up to quadratic terms of the parameters and therefore the parameters are Gaussian (see Appendix A).

In short, non-Gaussianity of parameters is a result of the data not being constraining enough for fitting parameters. When this is the case, there may not be a mathematical correspondence between the ML parameter distribution from the simulations and the posterior distribution of the MCMC chains, as our argument in Appendix A depends on the assumption of Gaussianity. Besides, because our proof of correlation equivalence in Appendix B also rests upon approximations of the log likelihood to only the second order derivatives, the property now becomes questionable. This means that we cannot simply estimate the correlation between parameters fit separately from the MCMC chains where they are fitted together. What is worse is that the means and the covariance matrix no longer carry all information of the parameter distributions and one might need to evaluate high order tensors (Sellentin et al. 2014). Additionally, the χ² test is inappropriate for non-Gaussian parameters and we need new ways (such as one proposed in Appendix C of Motloch & Hu 2019) to quantify the significance of the overall difference between the experimental parameters and their fiducial values.

An example of non-Gaussian parameters are those with priors that are more informative than the data, e.g., the neutrino mass, $\sum {m}_{\nu }$, and the tensor to scalar ratio, r, physically defined as non-negative. Current CMB data are not sufficiently constraining to pull these parameters away from lower bounds (BICEP2 Collaboration et al. 2018; Planck Collaboration VI 2018). Another non-Gaussian parameter given just the TT data, is the curvature density Ω_k. Curvature can only be weakly constrained as allowing it to be free worsens the existing degeneracy between the physical matter density Ω_m and the Hubble constant H₀ (Zaldarriaga et al. 1997; Percival et al. 2002; Kable et al. 2019). In the MCMC method, a two-parameter extension model with both A_L and Ω_k results in very wide and non-Gaussian distributions, as they are highly correlated: a positive curvature has similar effect as an increased lensing signal (Planck Collaboration XIII 2016).

To reduce non-Gaussianity, one can always include more data sets if available to have greater constraining power on the parameter, e.g., include the BAO data (Alam et al. 2017) in parameter fitting along with Planck. However, when extensions are used as a means to test the internal consistency of one data set, adding extra data is not an option. Fortunately, there are existing methods of transforming non-Gaussian parameters into Gaussian ones. One such method is the Box–Cox transformation (Box & Cox 1964). Joachimi & Taylor (2011) and Schuhmann et al. (2016) applied this bijective transformation to non-Gaussian cosmological parameters. Theoretically, correlation equivalence still holds for the transformed Gaussian parameters. Thus we may still apply the procedure outlined in Section 4.1 for transformed parameters, obtain the expected distribution for transformed ones, and then transform them back to the model parameters. One thing to note is that the Box–Cox transformation does not guarantee Gaussianity. So one should check for the Gaussianity of transformed parameters, e.g., calculate the skewness and kurtosis of the resulting distribution, and if needed, apply a second transformation. For consistency, it is also a good idea to compare the resulting distribution of model parameters to simulations, or compare the covariance of the transformed parameters with predictions from Fisher Forecast, keeping in mind that Fisher Forecast assumes Gaussianity.

Further understanding of non-Gaussian scenarios is left for future efforts.

With ${\mathfrak{L}}\equiv \mathrm{log}L({\boldsymbol{\theta }})$, the log likelihood of parameters given a specific sample of ${\boldsymbol{d}}$ can be written as

$$\begin{eqnarray}&&{\mathfrak{L}}=-\displaystyle \frac{1}{2}{\left({\boldsymbol{d}}-{\boldsymbol{\mu }}({\boldsymbol{\theta }})\right)}^{T}{{\boldsymbol{ \mathcal C }}}^{-1}({\boldsymbol{d}}-{\boldsymbol{\mu }}({\boldsymbol{\theta }}))\end{eqnarray} \tag{ 12 }$$

up to some constants. ${\boldsymbol{\mu }}({\boldsymbol{\theta }})$ is the theory prediction of the data given some set of parameters ${\boldsymbol{\theta }}$. For any ${\boldsymbol{d}}$, the goal is to find the ${\boldsymbol{\theta }}$ that maximize ${\mathfrak{L}}$.

As in Appendix A, at ML, the partial derivative of the log likelihood with respect to (w.r.t.) θ_i is zero (i runs over n, the number of parameters). That is

$$\begin{eqnarray}\begin{array}{rcl}{\left.\displaystyle \frac{\partial {\mathfrak{L}}}{\partial {\theta }_{i}}\right|}_{\mathrm{ML}} & = & {\left.\displaystyle \frac{\partial {\mathfrak{L}}}{\partial {{\boldsymbol{\mu }}}^{T}}\right|}_{\mathrm{ML}}{\left.\displaystyle \frac{\partial {\boldsymbol{\mu }}}{\partial {\theta }_{i}}\right|}_{\mathrm{ML}}\\ & = & \ {\left.-\displaystyle \frac{\partial {{\boldsymbol{\mu }}}^{T}}{\partial {\theta }_{i}}\right|}_{\mathrm{ML}}{{\boldsymbol{ \mathcal C }}}^{-1}({\boldsymbol{d}}-{\boldsymbol{\mu }}({\boldsymbol{\theta }}))=0.\end{array}\end{eqnarray} \tag{ 13 }$$

Close to the ML point in parameter space, we can Taylor expand the log likelihood to the second order in ${\rm{\Delta }}{\theta }_{i}\equiv {\theta }_{i}-{\theta }_{i}^{\mathrm{ML}}$ as shown in Equation (6), where θ_i is well described by a multivariate Gaussian distribution. The parameter means are the values that maximize the likelihood.

The parameter covariance ${\boldsymbol{\Sigma }}$ can be calculated as the inverse of the Fisher matrix ${{ \mathcal F }}^{-1}$, with the Fisher matrix, estimated from MCMC chains or calculated analytically, defined as

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal F }}_{{ij}} & = & {\left.-\displaystyle \frac{{\partial }^{2}{\mathfrak{L}}}{\partial {\theta }_{i}\partial {\theta }_{j}}\right|}_{\mathrm{ML}}\\ & = & \ {{\boldsymbol{\mu }}}_{,i}^{\mathrm{ML},{\rm{T}}}{{\boldsymbol{ \mathcal C }}}^{-1}{{\boldsymbol{\mu }}}_{,j}^{\mathrm{ML}}.\end{array}\end{eqnarray} \tag{ 14 }$$

From now on we use the comma notation to denote partial derivatives w.r.t. ${\boldsymbol{\theta }}$, that is, ${{\boldsymbol{\mu }}}_{,i}\equiv \tfrac{\partial {\boldsymbol{\mu }}}{\partial {\theta }_{i}}$.

In the following proof, we continue to work under the assumption of Gaussianity for ${\boldsymbol{\theta }}$ and expand ${\boldsymbol{\mu }}$ only to the first order of ${\boldsymbol{\theta }}$, that is, assuming the Gaussian linear model (Raveri & Hu 2019), around a chosen fiducial model ${{\boldsymbol{\theta }}}^{\mathrm{fid}}$:

$$\begin{eqnarray}&&{\boldsymbol{\mu }}({\boldsymbol{\theta }})\approx {{\boldsymbol{\mu }}}^{\mathrm{fid}}+\displaystyle \sum _{i=1}^{n}{{\boldsymbol{\mu }}}_{,i}^{\mathrm{fid}}({\theta }_{i}-{\theta }_{i}^{\mathrm{fid}}),\end{eqnarray} \tag{ 15 }$$

where we define ${{\boldsymbol{\mu }}}^{\mathrm{fid}}={\boldsymbol{\mu }}({{\boldsymbol{\theta }}}^{\mathrm{fid}})$. In the Gaussian linear model, the Jacobian ${{\boldsymbol{\mu }}}_{,i}$ may be taken as constant over changes in ${\boldsymbol{\theta }}$ and so ${{\boldsymbol{\mu }}}_{,i}^{\mathrm{fid}}$ is approximately equal to ${{\boldsymbol{\mu }}}_{,i}^{\mathrm{ML}}$. For simplicity, from here on, we just drop all the superscripts for ${{\boldsymbol{\mu }}}_{,i}$. In practice, the fiducial model is usually based on to the ML model given by the data and updated iteratively if necessary.

Next, we move on to calculate the elements in parameter covariance ${\boldsymbol{\Sigma }}$:

$$\begin{eqnarray}&&{{\boldsymbol{\Sigma }}}_{{ij}}={\left({{ \mathcal F }}^{-1}\right)}_{{ij}}.\end{eqnarray} \tag{ 16 }$$

Remember that for any invertible square matrix ${ \mathcal M }$, its inverse can be calculated as

$$\begin{eqnarray}&&{\left({{ \mathcal M }}^{-1}\right)}_{{ij}}={\left(-1\right)}^{i+j}\displaystyle \frac{| {{ \mathcal M }}_{\setminus j\setminus i}| }{| { \mathcal M }| }\end{eqnarray} \tag{ 17 }$$

where $| { \mathcal M }|$ is the determinant of ${ \mathcal M }$. We use the notation "${{ \mathcal M }}_{\setminus j\setminus i}$," to denote a smaller matrix, corresponding to ${ \mathcal M }$ with the jth row and the ith column removed (the backslash symbol is borrowed from the notation for set difference in set theory). Then $| {{ \mathcal M }}_{\setminus j\setminus i}|$ is a minor of ${ \mathcal M }$. The determinant for an n × n matrix ${ \mathcal M }$ can be written as

$$\begin{eqnarray}&&| { \mathcal M }| =\displaystyle \sum _{i=1}^{n}{\left(-1\right)}^{i+j}{{ \mathcal M }}_{{ij}}| {{ \mathcal M }}_{\setminus i\setminus j}| \end{eqnarray} \tag{ 18 }$$

for any 1 ≤ j ≤ n.

Therefore,

$$\begin{eqnarray}&&{{\boldsymbol{\Sigma }}}_{{ij}}={\left(-1\right)}^{i+j}| {{ \mathcal F }}_{\setminus j\setminus i}| /| { \mathcal F }| .\end{eqnarray} \tag{ 19 }$$

If we choose to fix the ith parameter, we can simply remove the ith row and column from the Fisher matrix, and calculate a new parameter covariance matrix from the revised Fisher matrix.

When two parameters, A and B, are both variables, we can calculate their covariance from the Fisher matrix that includes them along with the base ΛCDM parameters. So their correlation is just

$$\begin{eqnarray}&&\mathrm{corr}(A,B)=\displaystyle \frac{{{\boldsymbol{\Sigma }}}_{{AB}}}{{\left({{\boldsymbol{\Sigma }}}_{{AA}}{{\boldsymbol{\Sigma }}}_{{BB}}\right)}^{\tfrac{1}{2}}}.\end{eqnarray} \tag{ 20 }$$

From here on, for clarity we use A and B as subscripts for the rows and columns in the matrix for the two parameters of interest.

Following from Equations (20), (21) becomes

$$\begin{eqnarray}&&\mathrm{corr}(A,B)=\displaystyle \frac{-| {{ \mathcal F }}_{\setminus A\setminus B}| }{{\left(| {{ \mathcal F }}_{\setminus A\setminus A}| | {{ \mathcal F }}_{\setminus B\setminus B}| \right)}^{\tfrac{1}{2}}}.\end{eqnarray} \tag{ 21 }$$

The minus sign is from setting the index of B equal to that of A plus one.

For a generic set of parameters, given the data array ${\boldsymbol{d}}$ (experimental or simulated) and the Fisher matrix ${ \mathcal F }$, we can substitute (16) into (13) to obtain a relationship between the ML parameters ${{\boldsymbol{\theta }}}^{\mathrm{ML}}$ and the fiducial ones:

$$\begin{eqnarray}&&{{\boldsymbol{\theta }}}^{\mathrm{ML}}({\boldsymbol{d}})-{{\boldsymbol{\theta }}}^{\mathrm{fid}}={{ \mathcal F }}^{-1}{{\boldsymbol{\mu }}}_{,i}^{T}{{\boldsymbol{ \mathcal C }}}^{-1}({\boldsymbol{d}}-{\boldsymbol{\mu }}).\end{eqnarray} \tag{ 22 }$$

So we can express ${{\boldsymbol{\theta }}}^{\mathrm{ML}}$ as:

$$\begin{eqnarray}&&{{\boldsymbol{\theta }}}^{\mathrm{ML}}({\boldsymbol{d}})={{ \mathcal F }}^{-1}{{\boldsymbol{\mu }}}_{,i}^{T}{{\boldsymbol{ \mathcal C }}}^{-1}{\boldsymbol{d}}+\hat{{\boldsymbol{\theta }}},\end{eqnarray} \tag{ 23 }$$

where we use the vector $\hat{{\boldsymbol{\theta }}}$ to represent all the terms in (22) that are independent of data ${\boldsymbol{d}}$:

$$\begin{eqnarray*}&&\hat{{\boldsymbol{\theta }}}\equiv -{{ \mathcal F }}^{-1}{{\boldsymbol{\mu }}}_{,i}^{T}{{\boldsymbol{ \mathcal C }}}^{-1}{{\boldsymbol{\mu }}}^{\mathrm{fid}}+{{\boldsymbol{\theta }}}^{\mathrm{fid}}.\end{eqnarray*}$$

As we will show below, $\hat{{\boldsymbol{\theta }}}$ does not contribute to the correlation between A and B.

Again for clarity, we use A and B to index the rows and columns for parameters A and B, respectively, while using the generic i and j to index n base ΛCDM parameters. So A is the (n + 1)th parameter and B (n + 2)th.

When fixing B to its fiducial ΛCDM value B^fid, the expression for the best estimate for other parameters is almost the same as Equation (23), except that ${ \mathcal F }\to {{ \mathcal F }}_{\setminus B\setminus B}$ and the elements associated with B in ${{\boldsymbol{\mu }}}_{,i}^{T}$, ${{\boldsymbol{\theta }}}^{\mathrm{ML}}$ and $\hat{{\boldsymbol{\theta }}}$ are deleted. The expression for the rest of elements in ${{\boldsymbol{\theta }}}^{\mathrm{ML}}$ is:

$$\begin{eqnarray}&&\left(\begin{array}{c}{{\boldsymbol{\theta }}}_{1}^{\mathrm{ML}}\\ {{\boldsymbol{\theta }}}_{2}^{\mathrm{ML}}\\ \vdots \\ {{\boldsymbol{\theta }}}_{n}^{\mathrm{ML}}\\ {A}^{\mathrm{ML}}\end{array}\right)={\left({{ \mathcal F }}_{\setminus B\setminus B}\right)}^{-1}\left(\begin{array}{c}{{\boldsymbol{\mu }}}_{,1}^{T}{{\boldsymbol{ \mathcal C }}}^{-1}{\boldsymbol{d}}\\ {{\boldsymbol{\mu }}}_{,2}^{T}{{\boldsymbol{ \mathcal C }}}^{-1}{\boldsymbol{d}}\\ \vdots \\ {{\boldsymbol{\mu }}}_{,n}^{T}{{\boldsymbol{ \mathcal C }}}^{-1}{\boldsymbol{d}}\\ {{\boldsymbol{\mu }}}_{,A}^{T}{{\boldsymbol{ \mathcal C }}}^{-1}{\boldsymbol{d}}\end{array}\right)+\hat{{\boldsymbol{\theta }}}.\end{eqnarray} \tag{ 24 }$$

The last row of (25) gives

$$\begin{eqnarray}&&\begin{array}{l}{A}^{\mathrm{ML}}({\boldsymbol{d}})\\ \ =\ \left[\displaystyle \sum _{i=1}^{n}{\left({{ \mathcal F }}_{\setminus B\setminus B}\right)}_{{Ai}}^{-1}{{\boldsymbol{\mu }}}_{,i}^{T}{{\boldsymbol{ \mathcal C }}}^{-1}+{\left({{ \mathcal F }}_{\setminus B\setminus B}\right)}_{{AA}}^{-1}{{\boldsymbol{\mu }}}_{,A}^{T}{{\boldsymbol{ \mathcal C }}}^{-1}\right]{\boldsymbol{d}}+{\hat{\theta }}_{A}.\end{array}\end{eqnarray} \tag{ 25 }$$

Similarly, fixing A and letting B vary leads to

$$\begin{eqnarray}&&\begin{array}{l}{B}^{\mathrm{ML}}({\boldsymbol{d}})\\ \ =\ \left[\displaystyle \sum _{i=1}^{n}{\left({{ \mathcal F }}_{\setminus A\setminus A}\right)}_{{Bi}}^{-1}{{\boldsymbol{\mu }}}_{,i}^{T}{{\boldsymbol{ \mathcal C }}}^{-1}+{\left({{ \mathcal F }}_{\setminus A\setminus A}\right)}_{{BB}}^{-1}{{\boldsymbol{\mu }}}_{,B}^{T}{{\boldsymbol{ \mathcal C }}}^{-1}\right]{\boldsymbol{d}}+{\hat{\theta }}_{B}.\end{array}\end{eqnarray} \tag{ 26 }$$

The variance of A and B can be obtained from the Fisher matrix (14) and its minors (18), as follows:

$$\begin{eqnarray}&&\mathrm{Var}({A}^{\mathrm{ML}})={\left({{ \mathcal F }}_{\setminus B\setminus B}\right)}_{{AA}}^{-1}=\displaystyle \frac{\left|{{ \mathcal F }}_{\mathop{\setminus B\setminus B}\limits_{\setminus A\setminus A}}\right|}{| {{ \mathcal F }}_{\setminus B\setminus B}| }\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&\mathrm{Var}({B}^{\mathrm{ML}})={\left({{ \mathcal F }}_{\setminus A\setminus A}\right)}_{{BB}}^{-1}=\displaystyle \frac{\left|{{ \mathcal F }}_{\mathop{\setminus A\setminus A}\limits_{\setminus B\setminus B}}\right|}{| {{ \mathcal F }}_{\setminus A\setminus A}| }.\end{eqnarray} \tag{ 28 }$$

And the covariance between A^ML and B^ML is:

$$\begin{eqnarray}&&\mathrm{Cov}({A}^{\mathrm{ML}},{B}^{\mathrm{ML}})=\langle ({A}^{\mathrm{ML}}-\langle {A}^{\mathrm{ML}}\rangle )({B}^{\mathrm{ML}}-\langle {B}^{\mathrm{ML}}\rangle )\rangle ,\end{eqnarray} \tag{ 29 }$$

with the brackets here representing the averaging over all realizations of the data ${\boldsymbol{d}}$ .

Recall that the data covariance ${\boldsymbol{ \mathcal C }}$ is assumed to be fixed and with our assumption of a Gaussian Linear Model, all partial derivatives w.r.t. the parameters are also constant. Then in (26) and (27), only ${\boldsymbol{d}}$ is a variable. Inserting (26) and (27) into (30), we find that terms involving ${\hat{\theta }}_{A}$ and ${\hat{\theta }}_{B}$ cancel. In addition, all constant terms can be taken out of the brackets, leaving only $\langle {({\boldsymbol{d}}-\langle {\boldsymbol{d}}\rangle )}^{T}({\boldsymbol{d}}-\langle {\boldsymbol{d}}\rangle )\rangle$, which is just ${\boldsymbol{ \mathcal C }}$. This simplifies to

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{Cov}({A}^{\mathrm{ML}},{B}^{\mathrm{ML}})\\ \ =\ \displaystyle \sum _{i,j=1}^{n}{\left({{ \mathcal F }}_{\setminus B\setminus B}\right)}_{{Ai}}^{-1}{{\boldsymbol{\mu }}}_{,i}^{T}{{\boldsymbol{ \mathcal C }}}^{-1}{{\boldsymbol{\mu }}}_{,j}{\left({{ \mathcal F }}_{\setminus A\setminus A}\right)}_{{Bj}}^{-1}\\ \ \ \ +\ \displaystyle \sum _{i=1}^{n}{\left({{ \mathcal F }}_{\setminus B\setminus B}\right)}_{{Ai}}^{-1}{{\boldsymbol{\mu }}}_{,i}^{T}{{\boldsymbol{ \mathcal C }}}^{-1}{{\boldsymbol{\mu }}}_{,B}{\left({{ \mathcal F }}_{\setminus A\setminus A}\right)}_{{BB}}^{-1}\\ \ \ \ +\ \displaystyle \sum _{i=1}^{n}{\left({{ \mathcal F }}_{\setminus B\setminus B}\right)}_{{AA}}^{-1}{{\boldsymbol{\mu }}}_{,A}^{T}{{\boldsymbol{ \mathcal C }}}^{-1}{{\boldsymbol{\mu }}}_{,i}{\left({{ \mathcal F }}_{\setminus A\setminus A}\right)}_{{Bi}}^{-1}\\ \ \ \ \ +\ {\left({{ \mathcal F }}_{\setminus B\setminus B}\right)}_{{AA}}^{-1}{{\boldsymbol{\mu }}}_{,A}^{T}{{\boldsymbol{ \mathcal C }}}^{-1}{{\boldsymbol{\mu }}}_{,B}{\left({{ \mathcal F }}_{\setminus A\setminus A}\right)}_{{BB}}^{-1}.\end{array}\end{eqnarray} \tag{ 30 }$$

Notice that terms like ${{\boldsymbol{\mu }}}_{,i}^{T}{{\boldsymbol{ \mathcal C }}}^{-1}{{\boldsymbol{\mu }}}_{,j}$ are just elements of the Fisher matrix and we can also use Equation (18) to express terms like ${\left({{ \mathcal F }}_{\setminus B\setminus B}\right)}_{{Ai}}^{-1}$ in terms of determinants of minors of the Fisher matrix. Thus,

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{Cov}({A}^{\mathrm{ML}},{B}^{\mathrm{ML}})=\displaystyle \frac{1}{| {{ \mathcal F }}_{\setminus B\setminus B}| | {{ \mathcal F }}_{\setminus A\setminus A}| }\\ \ \ \times \ \left[\mathop{\underbrace{{\displaystyle \sum }_{i,j=1}^{n}{\left(-1\right)}^{i+n+1}\left|{{ \mathcal F }}_{\mathop{\setminus B\setminus B}\limits_{\setminus i\setminus A}}\right|{{ \mathcal F }}_{{ij}}{\left(-1\right)}^{j+n+1}\left|{{ \mathcal F }}_{\mathop{\setminus A\setminus A}\limits_{\setminus j\setminus B}}\right|}}\limits_{[1]}\right.\\ \ \ +\ \mathop{\underbrace{{\displaystyle \sum }_{i=1}^{n}{\left(-1\right)}^{i+n+1}\left|{{ \mathcal F }}_{\mathop{\setminus B\setminus B}\limits_{\setminus i\setminus A}}\right|{{ \mathcal F }}_{{iB}}\left|{{ \mathcal F }}_{\mathop{\setminus A\setminus A}\limits_{\setminus B\setminus B}}\right|}}\limits_{[2]}\\ \ \ +\ \mathop{\underbrace{{\displaystyle \sum }_{i=1}^{n}\left|{{ \mathcal F }}_{\mathop{\setminus A\setminus A}\limits_{\setminus B\setminus B}}\right|{{ \mathcal F }}_{{Ai}}{\left(-1\right)}^{i+n+1}\left|{{ \mathcal F }}_{\mathop{\setminus A\setminus A}\limits_{\setminus i\setminus B}}\right|}}\limits_{[3]}\\ \ \ +\ \left.\mathop{\underbrace{\left|{{ \mathcal F }}_{\mathop{\setminus B\setminus B}\limits_{\setminus A\setminus A}}\right|{{ \mathcal F }}_{{AB}}\left|{{ \mathcal F }}_{\mathop{\setminus A\setminus A}\limits_{\setminus B\setminus B}}\right|}}\limits_{[4]}\right].\end{array}\end{eqnarray} \tag{ 31 }$$

To simplify (32), first we combine term [1] and [2]:

$$\begin{eqnarray*}&&\begin{array}{l}[1]+[2]=\displaystyle \sum _{i=1}^{n}{\left(-1\right)}^{i+n+1}\left|{{ \mathcal F }}_{\mathop{\setminus B\setminus B}\limits_{\setminus i\setminus A}}\right|\\ \ \times \ \left[\displaystyle \sum _{j=1}^{n}{{ \mathcal F }}_{{ij}}{\left(-1\right)}^{j+n+1}\left|{{ \mathcal F }}_{\mathop{\setminus A\setminus A}\limits_{\setminus j\setminus B}}\right|+{{ \mathcal F }}_{{iB}}\left|{{ \mathcal F }}_{\mathop{\setminus A\setminus A}\limits_{\setminus B\setminus B}}\right|\right]\\ \ \ \ \ =\ \displaystyle \sum _{i=1}^{n}{\left(-1\right)}^{i+n+1}\left|{{ \mathcal F }}_{\mathop{\setminus B\setminus B}\limits_{\setminus i\setminus A}}\right|\\ \ \times \ \left[\displaystyle \sum _{j=1}^{n}{{ \mathcal F }}_{{ji}}{\left(-1\right)}^{(n+1)+j}\left|{{ \mathcal F }}_{\mathop{\setminus A\setminus A}\limits_{\setminus j\setminus B}}\right|+{{ \mathcal F }}_{{Bi}}{\left(-1\right)}^{2\times (n+1)}\left|{{ \mathcal F }}_{\mathop{\setminus A\setminus A}\limits_{\setminus B\setminus B}}\right|\right]\end{array}\end{eqnarray*}$$

where we have used the property that the Fisher matrix is symmetric.

Compared with Equation (18), note that the terms inside the last bracket above sum up to the determinant of an (n + 1) × (n + 1) matrix, which is the same as ${{ \mathcal F }}_{\setminus A\setminus A}$, except that its (n + 1)th column is replaced by its ith. So not all of the columns for this matrix are linearly independent, resulting in its determinant being zero. Thus, [1] + [2] = 0.

To simplify terms [3] and [4], we have

$$\begin{eqnarray}&&\begin{array}{l}[3]+[4]\\ \ =\ \left|{{ \mathcal F }}_{\mathop{\setminus A\setminus A}\limits_{\setminus B\setminus B}}\right|\left[\displaystyle \sum _{i=1}^{n}{{ \mathcal F }}_{{Ai}}{\left(-1\right)}^{i+n+1}\left|{{ \mathcal F }}_{\mathop{\setminus A\setminus A}\limits_{\setminus i\setminus B}}\right|+{{ \mathcal F }}_{{AB}}\left|{{ \mathcal F }}_{\mathop{\setminus A\setminus A}\limits_{\setminus B\setminus B}}\right|\right]\\ \ =\ \left|{{ \mathcal F }}_{\mathop{\setminus A\setminus A}\limits_{\setminus B\setminus B}}\right|\left[\displaystyle \sum _{i=1}^{n}{{ \mathcal F }}_{{iA}}{\left(-1\right)}^{i+n+1}\left|{{ \mathcal F }}_{\mathop{\setminus A\setminus B}\limits_{\setminus i\setminus A}}\right|+{{ \mathcal F }}_{{BA}}\left|{{ \mathcal F }}_{\mathop{\setminus A\setminus B}\limits_{\setminus B\setminus A}}\right|\right]\\ \ =\ \left|{{ \mathcal F }}_{\mathop{\setminus A\setminus A}\limits_{\setminus B\setminus B}}\right|\left|{{ \mathcal F }}_{\setminus A\setminus B}\right|\end{array}\end{eqnarray} \tag{ 32 }$$

Inserting (33) into (32), we have

$$\begin{eqnarray}&&\mathrm{Cov}({A}^{\mathrm{ML}},{B}^{\mathrm{ML}})=\displaystyle \frac{\left|{{ \mathcal F }}_{\mathop{\setminus A\setminus A}\limits_{\setminus B\setminus B}}\right|\left|{{ \mathcal F }}_{\setminus A\setminus B}\right|}{| {{ \mathcal F }}_{\setminus B\setminus B}| | {{ \mathcal F }}_{\setminus A\setminus A}| }.\end{eqnarray} \tag{ 33 }$$

Using Equations (28) and (29), the correlation between the best-fit values of parameters A and B estimated separately can be written as

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{corr}({A}^{\mathrm{ML}},{B}^{\mathrm{ML}}) & = & \displaystyle \frac{\mathrm{Cov}({A}^{\mathrm{ML}},{B}^{\mathrm{ML}})}{{\left(\mathrm{Var}({A}^{\mathrm{ML}})\times \mathrm{Var}({B}^{\mathrm{ML}})\right)}^{\tfrac{1}{2}}}\\ & = & \ \displaystyle \frac{\left|{{ \mathcal F }}_{\mathop{\setminus A\setminus A}\limits_{\setminus B\setminus B}}\right|\left|{{ \mathcal F }}_{\setminus A\setminus B}\right|}{| {{ \mathcal F }}_{\setminus B\setminus B}| | {{ \mathcal F }}_{\setminus A\setminus A}| }\displaystyle \frac{{\left(| {{ \mathcal F }}_{\setminus B\setminus B}| | {{ \mathcal F }}_{\setminus A\setminus A}| \right)}^{\tfrac{1}{2}}}{\left|{{ \mathcal F }}_{\mathop{\setminus A\setminus A}\limits_{\setminus B\setminus B}}\right|}\\ & = & \ \displaystyle \frac{| {{ \mathcal F }}_{\setminus A\setminus B}| }{{\left(| {{ \mathcal F }}_{\setminus A\setminus A}| | {{ \mathcal F }}_{\setminus B\setminus B}| \right)}^{\tfrac{1}{2}}},\end{array}\end{eqnarray} \tag{ 34 }$$

which is the same as Equation (22) up to a minus sign, thus completing our proof.