# 2D probe of multi-parameter functions: an MCMC alternative

Joel S. Jayson<EMAIL_ADDRESS>P.O. Box 34. Brooklyn, NY 11235

###### Abstract

We evaluate multi-parameter functions by fixing all parameter values, but two.
The temperature power spectrum of the lensed cosmic microwave background (CMB)
serves as an example in which the amplitude, $A_{s}$, nears linearity at small
deviations, reducing power spectra computation to a 1D problem. We apply the
analysis to assess the shift in the apparent value of ${H_{0}}$ as a result of
cosmic infrared background (CIB), and Poisson point sources (PS). By
iteratively cycling through parameters paired with $A_{s}$, and within a few
hundred calls for spectra, we derive values in agreement with Planck. We
adjudge that when neither variable is linear, thousands of calls are required,
still competitive with the MCMC method.

## Introduction.–

We evaluate multi-parameter functions when all parameter values, with the
exception of two are frozen. The temperature power spectrum of the cosmic
microwave background (CMB) provides a ready example for demonstration.

A compelling reason to look at CMB is that, of the six parameters needed to
define a spectrum, one of them is the amplitude, $A_{s}$. With no lensing
$A_{s}$ is proportional to the power spectra. With lensing the relationship is
not quite linear. We incorporate lensing in our analysis, but still assume
linearity, fixing $A_{s}$, and multiplying the power spectrum by a factor,
$\gamma$, which serves as the amplitude variable. As $\gamma$ approaches
unity, the introduced error approaches zero. This procedure affords us an
appreciable reduction in calls for power spectra computation, incrementing in
one dimension, rather than in a two dimensional grid. The ability to visualize
solutions is useful in two dimensional analysis, and after introducing the ln
Likelihood function and its partial derivatives, we depict pairs of these
partial derivatives. As intimated above, computed power spectra are stored,
and the various pairs formed by multiplying stored spectra by different values
of $\gamma$.

Galactic foregrounds, in particular, cosmic infrared background (CIB), and
Poisson point sources (PS) are assessed by taking advantage of the flexibility
provided in the selection of the fixed parameters. Thus, Addison, et al. [1]
evaluated the Hubble parameter from a power spectrum with maximum multipole,
1000, to be 69.7$\pm$1.7 km $\text{s}^{-1}$ $\text{Mpc}^{-1}$, as contrasted
with the Planck value computed from a power spectrum with maximum multipole,
2508, of 67.36$\pm$0.54 km $\text{s}^{-1}$ $\text{Mpc}^{-1}$. By maintaining
the fixed parameter values from the analysis with maximum multipole, 1000, for
an analysis with maximum multipole 2508, we establish a synthetic situation in
which all perturbations between multipole 1000, and multipole 2508 have been
eliminated. We then add the galactic foreground components into that gap, and
assess their effect on the value of the computed Hubble parameter.

To be able to solve for two parameter values, knowing the fixed values of the
others, is a potent tool, but its worth can most be appreciated when we
sequentially cycle through pairing each parameter with $\gamma$ in an
iterative process. By this stratagem we solve for the power spectra parameter
values. We have conducted four runs, each with a different set of initial
values, all producing similar results with small variance, and in agreement
withPlanck results [2] that were obtained using the Markov Chain Monte Carlo
(MCMC) method [3]. In two of the runs we introduced an algorithm,
$A_{s,k+1}=\gamma_{k}A_{s,k}$, where k is a parameter step, which brings about
a rapid convergence of $\gamma$ $\rightarrow$1\. The maximum number of calls
for power spectra computations for any of the runs was less than 600.

In applications where a parameter equivalent to $A_{s}$ is lacking, calls for
computing the associated function are incremented in a two dimensional grid.
We estimate that the number of required computational calls is on the order of
thousands. That number is nonetheless low enough for the analysis to hold its
own when compared with MCMC.

## 2D maximum likelihood estimate.–

The Planck Collaboration advanced the most accurate and precise CMB results
[2], and we use their data in this study. Though Planck observations covered
the entire sky, the data used to derive the CMB power spectrum excluded the
Galactic plane. The data was processed by the Planck Collaboration to provide
the equivalent of a full sky. At high multipoles the likelihood distribution
for a full sky tends to a Gaussian [4]. The Planck data and the computed power
spectra are binned as described in Planck 2018 results V [5]. We select a bin
size of $\Delta\ell$ =30 spanning multipoles 30-2508. That yields 82 full size
bins and a partial last bin. We omit multipoles 2-29 to approach a Gaussian
distribution, and to remove anomalies found in those first few multipoles,
including a power deficit [6]. An additional simplification is the omission of
polarization in the analysis. The approximate Gaussian distribution validates
using the ln likelihood function, $ln$$\mathcal{L}(p_{1},p_{2})$ [7],

$ln\mathcal{L}(p_{1},p_{2})=\sum_{i=2}^{{84}}\frac{(\mathcal{\hat{D}}_{i}^{TT}-\mathcal{D}_{i}^{TT}(p_{1},p_{2}))^{2}}{\hat{\sigma}_{i}^{2}}$
(1)

where,

$\mathcal{D}^{TT}_{i}(p_{1},p_{2})=\frac{\ell_{i}(\ell_{i}+1)\mathcal{C}_{i}^{TT}(p_{1},p_{2})}{2\pi}$
(2)

denotes the binned computed temperature power-spectrum multipole, entries
$\mathcal{\hat{D}}_{i}^{TT}$ and ${\hat{\sigma}_{i}^{2}}$ represent the binned
data points and corresponding variance, respectively, $\ell_{i}$, the weighted
multipole, and $\mathcal{C}_{i}^{TT}(p_{1},p_{2})$, the weighted multipole
moment, both within each bin, and $p_{1}$ and $p_{2}$ the two parameters that
we choose to vary. We use the CLASS code to evaluate temperature power spectra
[8].

Six parameters determine the power spectrum configuration, $h$ (as defined by
the Hubble parameter, $H_{0}=100h$ km $\text{s}^{-1}$ $\text{Mpc}^{-1}$),
$\Omega_{b}h^{2}$, the baryon fraction, $\Omega_{c}h^{2}$, the cold dark
matter fraction, $\tau$, the reionization optical depth, $n_{s}$, the spectral
index, and $A_{s}$, the power spectrum amplitude. Since the evaluation of
$\tau$ depends upon polarization data we adopt the value of 0.0543 found in
Ref.[2], and use that value throughout. We will later find it convenient to
substitute $\theta_{s}$, the angular scale of the first acoustic peak, for
$h$.

We initially choose $h$ as one of the variables. A key to rapid computation is
the selection of the power spectrum amplitude as the second variable. We have
incorporated lensing in our investigation, which negates the linear
relationship between $A_{s}$ and the power spectra. Nevertheless, for small
excursions from the nominal amplitude, linearity is assumed with little error.
That error will be eliminated when an iterative process is implemented to
evaluate the power spectra parameters. We refer to the analysis with $A_{s}$
as quasi-linear. $A_{s}$ is set to a constant value within
$\mathcal{D}_{i}^{TT}(h,A_{s})$ and we assume that
$\mathcal{D}^{TT}_{i}(h,A_{s})=\gamma y_{i}^{TT}(h)$, with all variation of
$A_{s}$ relegated to $\gamma$.

With that assumption, we take partial derivatives of Eq. (1) and setting them
to zero find,

$0=\sum_{i=2}^{84}\frac{(\mathcal{\hat{D}}_{i}^{TT}-\gamma
y_{i}^{TT}(h))}{\hat{\sigma}_{i}^{2}}\gamma\partial y_{i}^{TT}(h)/\partial h$
(3)

and,

$0=\sum_{i=2}^{84}\frac{(\mathcal{\hat{D}^{TT}}_{i}-\gamma
y_{i}^{TT}(h))}{\hat{\sigma}_{i}^{2}}y_{i}^{TT}(h).$ (4)

Figure 1: Plots of $\partial{ln\mathcal{L}(h,\gamma)}/\partial{(h)}$, and
$\partial{ln\mathcal{L}(h,\gamma)}/\partial{(\gamma)}$ versus $H_{0}$ for
three pairs of partial derivatives, each delineated by a specific value of
$\gamma$. The solution is found where a pair crossing lands exactly on the
zero line, $\gamma$=1.00007 in this instance. The dashed line illustrates that
the solution can be derived directly from the pairs that straddle the zero
line. The fixed parameter values are [2], $\Omega_{b}h^{2}$=0.02237,
$\Omega_{c}h^{2}$=0.1200, $n_{s}$=0.9649, and $\tau$=0.0543.

Fig. 1 depicts the results obtained using Planck parameter values [2]. The
Hubble parameter is plotted along the abscissa, while the values of the
partial derivatives of the ln Likelihood are plotted along the ordinate. Each
pair of the partials is defined by a specific value of $\gamma$. Much like
finding the root of a function with a single variable, once the zero line is
bracketed by two pairs, we have the equivalent of the intermediate value
theorem, and the solution on the zero line is readily found. The solution,
with $\gamma=1.00007$ was derived by this “bracketing and bisection” [7]
technique. However, the dashed line in the figure indicates that we could have
derived it directly from the pairs that straddle the zero line. The
intersection of the dashed line with the zero line furnishes the $h$ value,
and interpolation provides the $\gamma$ value of 1.00007. That intersection
and interpolation procedure is particularly useful when neither parameter
variable is linear, in that it minimizes the number of required power spectra
computation calls.

We turn to a problem of estimating how much the galactic foreground
components, cosmic infrared background (CIB), and Poisson point sources (PS),
might influence the measured value of $H_{0}$. The Planck collaboration has
gone to great lengths to subtract out those constituents. Thus, we have no
reference with regard to what trace, if any, remains.

Figure 2:  Plots of $\partial{ln\mathcal{L}(h,\gamma)}/\partial{(h)}$, and
$\partial{ln\mathcal{L}(h,\gamma)}/\partial{(\gamma)}$ versus $H_{0}$ for two
pairs of partial derivatives, one with no foreground addition, the other with
the addition of foreground. With no foreground addition, the solution for
Addison et al. [1], hypothetically placed within the multipole range 30-2508,
is much like that found in the 2-1000 multipole range. The fixed parameter
values, as estimated from Ref. [1] figure 1, are, $\Omega_{b}h^{2}$=0.02252,
$\Omega_{c}h^{2}$=0.115, and $n_{s}$=0.966. As throughout this paper,
$\tau$=0.0543. The foreground components for CIB and PS can be found in Ref.
[9], table 5. The appropriate weighting for each of those components is given
in Ref. [10], figure 16.

We take the following approach. Addison, et al.[1] have evaluated the CMB
power spectrum parameters with a maximum multipole of 1000, in contrast to the
Planck collaboration multipole maximum of 2508. In particular, they find a
value of $H_{0}$=69.7$\pm$1.7 km $\text{s}^{-1}$ $\text{Mpc}^{-1}$, as
compared to the Planck value of 67.36$\pm$ 0.54 km $\text{s}^{-1}$
$\text{Mpc}^{-1}$. If we assume that the Addison, et al.[1] values hold at a
multipole maximum of 2508, we have a hypothetical situation in which there are
no perturbations to the CMB power spectrum between multipole 1000 and
multipole 2508. Although the Planck data does not support this solution, we
have forced it by maintaining the values of the fixed parameters, which were
evaluated at a maximum multipole at 1000. We add the CIB, and PS components
found in that gap to the power spectrum. The values of these components at the
relevant frequencies,143 Ghz, and 217 Ghz, and their weighting can be found in
Planck collaboration publications [9], and [10], respectively. In Fig. 2 we
show the solutions when no foregrounds are added, and when the full complement
of CIB and PS components is added to the power spectra. The addition of the
foregrounds moves the computed value of $H_{0}$ substantially in the direction
of the Planck value.

The Atacama Cosmology Telescope (ACT) collaboration, working at multipoles up
to 4000, derived a value of $H_{0}$=67.9$\pm$1.5 km $\text{s}^{-1}$
$\text{Mpc}^{-1}$, in agreement with the Planck results [11]. As with Planck,
ACT also took extensive measures to subtract foregrounds from their data. In
view of those results, we have no explanation for the discrepancy between the
outcomes at multipole maximum 1000, and at multipole maximum 2508, and higher.

## Iterative solutions.–

Though until now we have focused on varying $h$, any of the other parameters,
i.e., $\Omega_{b}h^{2}$, $\Omega_{c}h^{2}$, and $n_{s}$ can also be paired
with $A_{s}$ through $\mathcal{D}^{TT}_{i}(p,A_{s})=\gamma y_{i}^{TT}(p)$.
Kosowsky, et al.[12] introduced the use of the angular scale of the first
acoustic peak, $\theta_{s}$, as one of the parameters, and subsequently that
was adopted in CMB analysis as a substitute for $h$ [13]. We do the same here.
$\theta_{s}$ is the CLASS code nomenclature, and is similar, but not identical
with that defined by CAMB $\theta_{*}$.

Table 1: A sample from an iteration run that converged after 31 iterations. Shown are the initial parameter values, the results of the first two iterations, the result of the 15th iteration, and the results of the final two iterations. See text for further details. | 100$\theta_{s}$ | $\Omega_{b}h^{2}$ | $\Omega_{c}h^{2}$ | $n_{s}$
---|---|---|---|---
|  $10^{9}A_{s}$/$\gamma$ |  $10^{9}A_{s}$/$\gamma$ |  $10^{9}A_{s}$/$\gamma$ |  $10^{9}A_{s}$/$\gamma$
Initial values | – | 0.03 | 0.1 | 1.2
Iteration #1 | 1.05190 | 0.01758 | 0.0708 | 1.1658
| 1.050275/1.78167 | 1.871243/0.9935 | 1.859080/0.96855 | 1.800611/1.0141
Iteration #2 | 1.04372 | 0.02331 | 0.0649 | 1.1246
| 1.826000/0.9866 | 1.801532/0.99417 | 1.791029/0.99323 | 1.778903/1.0174
| … | … | … | …
Iteration #15 | 1.04225 | 0.02251 | 0.1160 | 0.9744
| 2.074503/0.9998 | 2.074088/0.99995 | 2.073984/1.00165 | 2.077406/1.0012
| … | … | … | …
Iteration #30 | 1.04180 | 0.02227 | 0.1203 | 0.9649
| 2.101274/0.999995 | 2.101263/0.999995 | 2.101253/1.00004 | 2.101337/0.99997
Iteration #31 | 1.04180 | 0.02227 | 0.1203 | 0.9649
| 2.101274/0.999995 | 2,101263/0.99999 | 2.101242/1.00003 | 2.101305/0.99998

Sequentially computing the 2D solution for each parameter, and substituting
the latest value into $\mathcal{D}^{TT}_{i}(p,A_{s})=\gamma y_{i}^{TT}(p)$ is
an iterative process. We have conducted four iterative runs, each with a
different set of initial parameter values. In two of the runs (one with 10
iterations to convergence, the other with 16) we retain the Planck value of
$A_{s}=2.100549\text{x}10^{-9}$, and within a specified range select the other
parameter values at random. For the other two runs ( 19 and 31 iterations), we
introduce a significant modification by updating $A_{s}$ at each parameter
step with the algorithm $A_{s,k+1}$=$\gamma_{k}$$A_{s,k}$, where k is the
parameter step number. That enables us to arbitrarily set an initial value for
$A_{s}$, and add the requirement that $\gamma$ $\rightarrow$1 at convergence.
In Table 1 we provide a sample from one of those runs (31 iterations).
100$\theta_{s}$ is computed first, hence the lack of an initial value. The
initial value for $A_{s}$ was set at half that of the fiducial Planck value
(By contrast, the $A_{s}$ initial value in the second run was set at twice the
fiducial value. Initial values of the other parameters were likewise chosen to
be dissimilar). The four runs yield the following results,
100$\theta_{s}$=1.04179$\pm$0.000026, $\Omega_{b}h^{2}$=0.02227$\pm$0.000017,
$\Omega_{c}h^{2}$=0.1205$\pm$0.00026, $n_{s}$=0.9646$\pm$0.00057, and
$10^{9}A_{s}$=2.101971$\pm$0.0013. We also extract values of $H_{0}$ from each
run, and find, $H_{0}$=67.58$\pm$0.1 km $\text{s}^{-1}$ Mpc. The variance
between runs is tight. Despite our several simplifications, all parameter
values are in close agreement with the Planck results [2].

Figure 3:  Plots of
$\partial{ln\mathcal{L}(\theta_{s},\Omega_{b}h^{2})}/\partial{(\Omega_{b}h^{2})}$,
and
$\partial{ln\mathcal{L}(\theta_{s},\Omega_{b}h^{2})}/\partial{(\theta_{s})}$
versus $\Omega_{b}h^{2}$ for two pairs of partial derivatives. All points for
both pairs require several calls for power spectra computation. The pairs
straddle the zero line, and the intercept connecting the crossing of each
pair, and interpolation, provide the solution, $\Omega_{b}h^{2}$=0.02171, and
$\theta_{s}$=1.03994. The connecting line is omitted in the figure for
clarity. This example was taken after two iterations of a quasi-linear run, to
assess any differences. The results with $\gamma$ as a variable were
$\Omega_{b}h^{2}$=0.02175, and $\theta_{s}=1.03996$, virtually identical.

In the general case, creating a family of ln Likelihood partials cannot be
achieved by simple multiplication of one pair, and each pair requires multiple
calls of the power spectra, or an analogous entity for a non-CMB application.
Maintaining the focus on CMB, the metric by which we judge the effectiveness
of a computational method is the number of calls to generate power spectra.
For the example illustrated in Table 1 fewer than 600 calls were required. The
other three runs each required 400 or fewer calls. Symmetric numerical
differentiation is implemented to evaluate the partial derivatives of the ln
Likelihood. For our quasi-linear approach, an isolated point requires three
power spectra calls. Generally, within a neighborhood, we evaluate equally
spaced points to minimize the number of calls used exclusively to compute
derivatives. When searching for the cross-point of the two partials, we
routinely create a pair of points to determine which direction to further
proceed. That necessitates four power spectra calls. Fig. 3 depicts a two
parameter context in which 100$\theta_{s}$, and $\Omega_{b}h^{2}$ are the two
variables, and we no longer have a quasi-linear condition. The pairs of points
in that figure require eight power spectra calls for each pair, as derivatives
must be generated for both components of the partials. As indicated earlier,
once the zero line is straddled by two pairs, the intercept of the line
connecting the cross-points of each pair with the zero line identifies one
parameter value, and interpolation identifies the other. That is one of two
mitigating factors with regard to the required number of power spectra calls
in the general case. The second, is that after the first few iterations, the
parameter values have migrated to the vicinity of their final values. Thus, in
the example of Table 1, all iterations after the fifth required less than 20
calls, that is 4-6 calls per parameter. Consequently, our estimate of the
number of power spectra calls for deriving the CMB parameters with variables,
other than $A_{s}$, though significantly higher than under quasi-linear
circumstances, in the thousands, rather than hundreds, is not prohibitive.
Other applications require an assessment of their specific circumstances. We
posit that in general, two dimensional analysis compares favorably with the
MCMC method.

## References

  * [1] G. E. Addison, Y. Huang, D. J. Watts, C. L. Bennett, M. Halpern, G. Hinshaw, and J. L. Weiland, Quantifying discordance in the 2015 Planck CMB spectrum, Astrophys. J., 818, 132 (2016)
  * [2] Planck Collaboration, Planck 2018 results VI Cosmological parameters, Astron. Astrophys. 641, A6 (2020)
  * [3] W.K. Hastings, Monte Carlo sampling methods using Markov chains and their applications, Biometrika, 57, 97 (1970)
  * [4] W. J. Percival and M. L. Brown, Likelihood techniques for the combined analysis of CMB temperature and polarization power spectra, Mon. Not. R. Astron. Soc. 372, 1104 (2006)
  * [5] Planck Collaboration, Planck 2018 results V. CMB power spectra and likelihoods, Astron. Astrophys. 641, A5 (2020)
  * [6] Planck collaboration, Planck 2013 results. I Overview of products and scientific results, Astron. Astrophys. 571, A1 (2014)
  * [7] William H. Press, Saul A Teukolsky, William T. Vetterling, and Brian P. Flannery, Numerical Recipes, 3rd Edition (Cambridge University Press, Cambridge, 2007)
  * [8] Diego Blas, Julien Lesgougues, and Thomas Tram, The Cosmic Linear Anisotropy Solving System (CLASS). Part II: Approximation schemes, JCAP, 7, 034 (2011)
  * [9] Planck Collaboration, Planck 2013 results. XVI. Cosmological parameters, Astron. Astrophys. 571, A16 (2014)
  * [10] Planck Collaboration, Planck 2015 results XI. CMB power spectra, likelihoods, and robustness of parameters, Astron. Astrophys. 594, A11 (2016)
  * [11] Steve K. Choi, Matthew Hasselfield, Shuay-Pwu Patty Ho, Brian Koopman, Marius Lungu, Maximilian H. Abitbol, Graeme E. Addison, Peter A. R. Ade, Simone Aiola, David Alonso, et al., The Atacama Cosmology Telescope: a measurement of the Cosmic Microwave Background power spectra at 98 and 150 GHz, JCAP 12, 045 (2020)
  * [12] Arthur Kosowsky, Milos Milosavljevic, and Raul Jimenez, Efficient cosmological parameter estimation from microwave background anisotropies, Phys. Rev D 66, 063007 (2002)
  * [13] L. Verde, H. V. Peiris, D. N. Spergel, M. R. Nolta, C. L. Bennett, M. Halpern, G. Hinshaw, N. Jarosik, A. Kogut, M. Limon, et al., First-year Wilkinson microwave anisotropy probe (WMAP) observations: parameter estimation methodology, Astrophys J Suppl Ser, 148, 195 (2003)