​

# 以下是示例代码，用于展示相关计算过程（非完整代码）
# 假设已经有了协方差矩阵 C 和其他相关参数定义
import numpy as np
# 进行 Cholesky 分解
L = np.linalg.cholesky(C)
# 生成随机数（这里简化示意，实际可能更复杂）
delta_p0 = np.random.normal(0, j**2 / dj, (dj,))
# 投影回原始参数空间
delta_p = np.dot(L, delta_p0)

​
没有精确的规则来确定过采样因子 F2,⋯,FN。对于更快的似然函数和/或块中更多的干扰参数，应增加这些因子。如果数量过低，则无法享受快慢参数分解的优势。如果数量过高，则在快速参数的（∑jFjdj−Ns）次迭代中花费的时间可能与在慢速参数的 Ns次迭代中花费的时间相比显著增加，并且参数结果的收敛将被延迟。例如，对于普朗克似然函数的干扰参数，我们通常应用 4 的过采样因子。 ​

自动跳跃因子自适应还将使代码在从非常差的提议密度开始时更强大（例如，当向先前的运行添加许多新的自由参数时，或者当仅在约束性差得多的集合的协方差矩阵可用时研究一组非常有约束性的似然函数时）。为此，我们在 v3.0.0 版本中添加了新选项–superupdate，它与 –update互补：它们通常应结合使用，但由于–update默认已激活，因此只需传递–superupdate 标志。 请注意，过去已经研究了其他更新完整提议密度（而不仅仅是协方差矩阵）的方案。在对单链尝试此方法后，我们将其自由地改编为多链情况，以一种与传统协方差矩阵更新方案（–update）兼容的方式。 总体策略：当使用 –superupdate SU运行时，代码从可以通过命令标志手动设置的跳跃因子开始（例如–f 2.2），但如果未传递任何值，则默认使用 2.4。代码记录最后 $SU$ 个周期（即最后 $SU\times FPM$ 步）的接受率 a.r. 和跳跃参数 $c$。它还跟踪这些最后 $SU$ 个周期的平均 a.r. 和 $\bar{c}$。此信息用于计算开始和停止标准。由于 –superupdate要求–update` 处于活动状态，因此运行可以分为几个“更新序列”，这是两个连续更新之间的步骤集合。第一个“更新序列”只是直到第一次更新的时间。超级更新的基本原理是根据递归关系在每一步调整跳跃因子：

---

---

​

# 示例代码示意（非完整代码）
# 假设相关参数已定义
k_update = 0 # 第一个序列的起始步编号
for k in range(1, total_steps + 1):
 if k >= k_update + SU * FPM:
 # 计算跳跃因子更新（这里简化示意，实际有更多逻辑）
 j_new = j * (target_a.r. / a.r.) ** (1.0 / (k - k_update))

其中 kk 是当前步编号，而 kupdatekupdate 是每个新“更新序列”的第一步编号。对于第一个序列，kupdate=0kupdate=0。此递归关系导致在每个新序列开始时更快地更新，并在一段时间后具有更慢的更新和安全的收敛特性。 ​

​

# 示例代码示意（非完整代码）
# 假设相关参数已定义
N = total_parameters # 自由（慢 + 快）参数的数量
c_n1 = c_n * (np.linalg.det(C_n) / np.linalg.det(C_n1)) ** (1.0 / N)

其中 N 是自由（慢 + 快）参数的数量。当协方差矩阵的演变仅来自一个或几个参数时，这种重新缩放可能不是非常有效，但总的来说，这是可以做出的最佳简单猜测。就跳跃因子而言，这给出： ​

​

# 示例代码示意（非完整代码）
# 假设相关参数已定义
j_n1 = j_n * (np.linalg.det(C_n) / np.linalg.det(C_n1)) ** (1.0 / (2 * N))

请注意，对于协方差矩阵的第一次更新，这种重新缩放背后的逻辑不成立。实际上，如果我们从较差的输入协方差矩阵开始，跳跃因子的第一次重新缩放可能完全不切实际。为了安全起见，在第一次更新时，我们将跳跃因子重置为输入值（通过 --f 提供 [默认：2.4]）。
停止跳跃因子更新：我们调整跳跃参数，直到满足三个条件： ​

​

•  接受率应收敛到 26%，容差为 1 个百分点（在许多情况下，接受率开始较低并增加到最佳值，然后当代码达到 25% 时适应将停止）：

# 示例代码示意（非完整代码）
# 假设相关参数已定义
if abs(a.r. - 0.26) <= 0.01:
 stop_jump_factor_update = True

• 除了接受率标准之外，为了停止跳跃参数的适应，我们还要求它是稳定的：

# 示例代码示意（非完整代码）
# 假设相关参数已定义
if abs(c - c_mean) <= tolerance: # tolerance 为设定的稳定阈值
 stop_jump_factor_update = True

其中 c¯是最后 SU×FPM 步的跳跃参数的平均值。此外，我们不希望允许跳跃参数收敛到任意低的值，因为链陷入局部最小值的风险会增加。因此，我们引入跳跃因子的最小值，对应于初始跳跃因子的 10%。如果需要较小的跳跃参数，建议使用 --f 输入较低的值。

• 最后，我们要求从链的后半部分计算的最大（R−1）对于所有参数低于 0.4：这与停止协方差矩阵更新的条件相同。因此，超级更新机制将仅在最终的“更新序列”中停止其活动，在此期间可以积累大量真正的马尔可夫步骤（使用恒定的提议密度生成）。
  非 MPI 用户友好性：我们在蒙特卡罗 Python 中以与更新机制相同的编码原则实现了超级更新机制。因此，由于链之间的通信通过文件而不是 MPI 命令进行，它也可以在有或没有 MPI 的情况下使用。在运行目录中，文件 jumping_factors.txt 存储了所有已使用的跳跃因子序列，而文件 jumping_factor.txt 仅包含最终的跳跃因子，可在下一次运行中用作输入值。当使用 --restart 命令在同一目录中重新启动链时，这将自动完成。
  在最终结果中仅保留马尔可夫步骤：当 --superupdate 被激活时，代码仍会在每个新“更新周期”开始时在链中写入一些注释行，其中包含有关当前最大（R−1R−1）、jj 和 a.r. 值的信息。此外，当跳跃因子更新停止时，它会写入一行额外的注释。当在 info 模式下分析链时，除非用户传递选项 --keep - non - markovian，否则之前的所有行都将被丢弃，最终的数字和绘图将基于纯粹的马尔可夫链。
  单链运行的替代实现：超级更新机制原则上需要多个链，因为它使用基于 Gelman – Rubin 统计量的收敛测试。对于单链运行，代码将链拆分为三个单独的链以计算 Gelman – Rubin 统计量，这种做法可能比运行多个链不太可靠。此其他机制由标志 --adaptive 而不是 --superupdate 激活。它不使用 Gelman – Rubin 统计量。

​

​

著名的 Fisher 矩阵由有效 χ2相对于模型参数在 χ2χ2 最小值（即似然最大值）处的二阶导数构建而成：

# 示例代码示意（非完整代码）
# 假设相关函数和参数已定义
def fisher_matrix(parameters):
 # 计算二阶导数（这里简化示意，实际计算可能更复杂）
 second_derivatives = calculate_second_derivatives(parameters)
 return second_derivatives

根据最大似然点的定义，Fisher 矩阵必须是正定且可逆的。其逆矩阵是似然在最佳拟合点附近的高斯近似的协方差矩阵。如果二阶导数矩阵不是在该点计算的，则它可能不可逆。

（一）计算 Fisher 矩阵的动机

我们在Python v3.0.0 中使用有限差分法直接从似然和公式（8）计算 Fisher 矩阵，我们将在下一节详细介绍。此计算背后有两个动机： ​

​

•  对于使用模拟数据的参数预测，我们通常使用基准光谱作为观测光谱，而不生成随机实现。因此，最大似然恰好与基准模型重合，用户事先知道该模型。然后，新的 Fisher 方法效果特别好。
• 对于从真实数据中提取参数，我们最多知道最佳拟合点的近似值。然后，人们可能希望，如果真实和近似最佳拟合点之间的距离与有限差分法中使用的步长相比很小，则在后一点计算的近似 Fisher 矩阵仍将是正定且可逆的。然而，我们发现这种情况并不经常发生，因此将这种新方法用于真实数据的可能性仍然有些不确定。为了增加机会，我们在蒙特卡罗 Python 中纳入了一些来自 optimize Python 库的最小化算法，这至少可以帮助在尝试 Fisher 矩阵计算之前更接近真实的最佳拟合点。然而，我们的测试表明，目前实现的最小化算法不是非常稳健，特别是在存在许多干扰参数的情况下。
  总之，新的 Fisher 计算肯定会改善所有 MCMC 预测并且它也可能改善使用真实数据的 MCMC 运行，除非最终得到的似然形状对于找到最小值和/或运行 Fisher 算法来说过于复杂。

2. 替代 MCMC 运行

当知道给定运行的后验应该接近高斯分布，或者当对后验的细节不感兴趣（例如具有一些偏度、峰度、香蕉形状等的非平凡参数相关性）时，很想用简单的 Fisher 矩阵计算来替代整个 MCMC 参数提取运行。然后，逆 Fisher 矩阵将给出一些近似的一维置信区域和二维椭圆轮廓。这对于灵敏度预测尤其直接，因为在那种情况下，最大似然点事先已知。如果最大似然点已知到良好的近似值，例如在使用蒙特卡罗 Python 的新 --minimize 选项运行之后（其成功并非保证），也可以设想用于真实数据。
在文献中，绝大多数参数预测基于 Fisher 矩阵计算。这些通常由特定代码执行，利用在几步解析计算之后，FijFij 可以重新表示为可观测量相对于参数的导数的函数（例如 ∂C‘/∂pi 或 ∂P(k)/∂pi）。相反，蒙特卡罗 Python 执行的 Fisher 矩阵计算是基于直接似然评估，因为我们计算 ∂lnL/∂pi。这两种方法在数学上是等价的，但后者可能提供一些实际优势。实际上，主要用于模拟给定实验的量是似然。跳过导致像 ∂C‘/∂pi或 ∂P(k)/∂pi这样的导数的解析步骤有时可以避免复杂的表达式、引入近似的需要以及进一步出错的风险。
在这两种方法中，都必须使用给定的步长计算一些数值导数。对于纯高斯似然，步长应该无关紧要，只要它不是太小以至于数值误差（来自玻尔兹曼代码或似然代码）开始占主导地位。通常，玻尔兹曼代码经过优化，以便为最具约束性的实验（通常是现在的普朗克实验）提供 χ2effχeff2 的精度约为 δχ2∼O(10−1)，因为实现更好的精度不会改变置信区间的结果，因此将是浪费计算时间。因此，使用使得：

# 示例代码示意（非完整代码）
# 假设相关参数已定义
if step_size < 0.1 / some_factor: # some_factor 根据具体情况确定
 step_size = 0.1 / some_factor

显著小于 0.1 的步长 Δpi是危险的。这大致提供了 ΔpiΔpi 的下限。上限的问题更加微妙，特别是当似然不是高斯分布时：不同的选择可能会返回显著不同的 Fisher 矩阵和置信限。在这个问题上，社区存在不同的方法。一种观点建议使用尽可能小的步长，直到数值噪声起作用，以便尽可能接近二阶导数的数学定义。另一种观点更喜欢选择步长，使得如果最终目标是为参数提供 68%（或 95%）置信限的预测，则 Δχ2≈1（或 4），因为在那种情况下，Fisher 矩阵给出的似然的高斯近似恰好在与参数界限相关的区域内有效。许多 Fisher 代码甚至没有针对 Δχ2 的任何特定给定数量级，并且任意选择步长 Δpi。在 Python v3.0.0 中，我们通过让用户选择 ΔlnL=2Δχ2的目标值来解决这个问题。默认情况下，代码将首先尝试使用 ΔlnL∼0.1 获得可逆的 Fisher 矩阵，并在结果不可逆的情况下迭代增加该值，如下一节所述。然而，用户可以选择 ΔlnL的第一个值，例如可以将其设置为 0.5 或 2（使用标志 --fisher - delta）。

（二）从逆 Fisher 矩阵绘制似然轮廓

拥有 Fisher 矩阵及其逆矩阵后，用户可以轻松编写一个小脚本，绘制对应于参数空间中二维置信水平的椭圆。蒙特卡罗 Python 本身无法做到这一点，它只能在存在一些 MCMC 链时进行绘图。然而，蒙特卡罗 Python 的绘图工具对于比较 MCMC 运行结果与逆 Fisher 矩阵给出的高斯后验近似很有用。为此，可以使用通常的 info 模式分析蒙特卡罗 Python 结果，只需添加一个输入标志：--plot - fisher。然后代码将检查是否在用户尝试绘制的链所在的同一目录中计算并存储了 Fisher 矩阵。如果是这种情况，Fisher 椭圆将绘制在 MCMC 轮廓之上，如图 所示。

    y2 = info.extent[3]
        ax2dsub.fill_between(x95,
                             y1,
                             y2,
                             facecolor=contour_color[0],
                             edgecolor=contour_color[1],
                             linewidth=1,
                             alpha=contour_alpha)
        ax2dsub.fill_between(x68,
                             y1,
                             y2,
                             color=contour_color[1],
                             alpha=contour_alpha)

​
此图显示了基于来自 DESI 调查的模拟 BAO 数据与普朗克数据相结合的灵敏度预测，用于具有大质量中微子和动态暗能量的模型（νwCDM，8 个自由参数）。实际上，Fisher 矩阵在启动链之前计算，其逆矩阵用作输入协方差矩阵。最终的 MCMC 轮廓证明在这种情况下，Fisher 近似非常出色。这个逆 Fisher 矩阵不仅为 MCMC 运行提供了良好的提议密度，还给出了参数界限的出色估计，并且可以替代整个 MCMC 结果。可能会出现这样的情况：首先在最佳拟合点的猜测周围计算 Fisher 矩阵，然后使用它来启动将以真实最佳拟合点为中心的 MCMC 链（在高斯后验的情况下）。
在比较图中，可能会发现 Fisher 椭圆形状正确，但相对于真实最佳拟合点有偏移。为了获得更好的图，蒙特卡罗 Python 用户可以使用输入标志 --center - fisher。这将自动将 Fisher 椭圆中心定位在从 MCMC 链中提取的最大似然点上，而不是使用在 log.param 文件中读取的中心值，即使 Fisher 矩阵实际上是在该点计算的。

性能示例

为了说明超级更新的性能以及使用逆 Fisher 矩阵作为输入协方差矩阵的影响，我们选择了一些数据集和模型，并在使用或不使用这些不同选项的情况下进行了一些拟合。在最困难的情况下进行比较特别有趣：大量自由参数、先验知识少（即对输入协方差矩阵的猜测不佳）等。

（一）先验知识少的预测

我们首先对来自 DESI 调查的模拟 BAO 数据和普朗克数据的组合进行一些 MCMC 预测。我们使用名为 fake_desi_vol 的模拟 DESI 似然函数（在附录 D 中有文档记录）。对于预测，我们不需要使用真实的普朗克数据。相反，我们使用一种大致模拟普朗克卫星近似情况的似然函数，但使用与普朗克最佳拟合模型对应的一些合成数据。这个似然函数称为 fake_planck_realistic。
​
表 1：对于模拟数据和几个模型，使用 Gelman – Rubin 收敛准则比较三种采样选项。
我们使用最小的 6参数 ΛCDMΛCDM 模型以及几个具有多达 12 个自由参数的扩展模型来拟合这些数据集。这些扩展列在表 1 的第一列中，包括大质量中微子、具有恒定状态方程的动态暗能量、额外的相对论自由度、原始谱指数的跑动、空间曲率，以及最后具有 CPL 参数化的动态暗能量。我们在这些不同情况下使用 Metropolis – Hastings 算法和三种不同方法运行蒙特卡罗 Python：

• update [--update]：协方差矩阵的定期更新，
• superupdate [--superupdate]：跳跃因子的额外自适应，
• superupdate + Fisher [首先 --method Fisher；然后 --superupdate]：相同，但从蒙特卡罗 Python 计算的逆 Fisher 矩阵开始。
  对于更新和超级更新运行，提议密度初始化为“普朗克 2015 协方差矩阵”，即作为蒙特卡罗 Python 包公开分发的协方差矩阵，该矩阵源自基于普朗克 2015 似然函数的良好收敛运行的分析，并假设 6 参数 ΛCDM 模型。此外，在这两次运行中，跳跃因子最初设置为 2.4（因此在更新方法中它保持等于此值）。对于 Fisher 矩阵计算，我们将用于生成模拟数据的精确最佳拟合模型传递给代码。对于每个模型和方法，我们使用 8 个链启动代码，其中每个链在 6 个核心上运行，总共使用 48 个核心。在 12 或 48 小时后，我们计算所有参数上的较差收敛准则（R−1），去除每个链的初始 10 – 20%（取决于燃烧阶段的持续时间，但对于给定的模型和实验组合始终相同）。
  这些运行的困难在于对输入协方差矩阵的不良猜测。在六参数运行中，输入协方差矩阵仅来自普朗克数据，而 DESI BAO 数据非常具有约束性。这意味着提议密度最初太宽，需要收缩到 DESI 数据允许的小区域。当添加额外参数时，情况更糟。对于额外参数，代码不依赖于输入协方差矩阵，而是依赖于输入文件中写入的标准偏差（对于这些我们插入普朗克误差条）。因此，提议密度需要学习这些新方向上跳跃的正确数量级以及涉及额外参数的参数相关性。
  我们发现，对于最简单的模型（6 和 8 个参数），所有三种方法都成功地至少获得了 R−1=0.03的收敛，尽管超级更新和超级更新 + Fisher 表现更好，获得的 R−1小达 2 倍。对于更复杂的模型（10、11 和 12 个参数），从 Fisher 矩阵开始并使用超级更新产生了很大的差异，因为只有从 Fisher 矩阵开始的运行在限制我们仅运行 12 小时的情况下设法获得了任何收敛水平。然而，如果我们允许更长的运行时间（48 小时），更新和超级更新方法也设法收敛，这得益于协方差矩阵的定期更新，超级更新和超级更新 + Fisher 运行显示出比单独更新好 2.6 到 3.7 倍的收敛性。图 2 明确显示了为什么跳跃因子自适应和 Fisher 矩阵计算会导致此运行收敛速度显著加快。
• 
  ​
  图 2：表 1 中最后一次运行（10 参数模型，模拟普朗克 + DESI，48 小时）的接受率、跳跃因子和收敛估计量的演变，使用 U=50，SU=20，FPM=10，因此 Nupdate=500。跳跃因子信息在每一步更新，而关于（R−1）的信息每 Nupdate步更新（在链的最后 50% 上计算）。此信息始终可以从代码输出中提取。相反，上图中的局部接受率是为了绘制此图而对链进行后处理计算的，并且通过在大约 500 步上平均来定义（因此这与超级更新算法使用的量 a.r. 略有不同，后者仅在 SU×FPM=200 步上平均）。“Fisher”运行从一开始就基本上捕获了正确的协方差矩阵、跳跃因子和接受率。“超级更新”运行减小其跳跃因子以快速积累许多点并获得良好的协方差矩阵估计；一旦完成，它会增加跳跃因子以避免接受率过大。最后，“更新”运行在进入具有良好协方差矩阵和接受率的有效采样制度之前需要大约 3500 多个步，这在我们的 48 个核上对应于大约 12 小时：因此我们可以说在这个特定示例中，“超级更新”为单次运行节省了大约 600 小时。

（二）当前数据

为了比较我们的新方法在使用当前数据时的效率，我们仅考虑 6 参数 ΛCDM 模型，我们将其拟合到两个数据集：一个仅包含普朗克和 BAO 似然函数的小数据集，以及一个更大的数据集，包括大尺度结构（LSS）似然函数（来自 SDSS 的星系团和来自 CFHTLenS 的弱引力透镜）。更多详细信息在表 2 中给出。我们再次以三种不同方式（更新、超级更新、超级更新 + Fisher）进行拟合，与上一节完全相同（即从代码分发的“普朗克 2015 协方差矩阵”开始更新和超级更新，并从跳跃因子 2.4 开始）。我们使用与预测相同数量的链和核心，并允许链运行 12 或 48 小时。
这些运行说明了从提议密度的不良猜测开始的情况，因为输入协方差矩阵仅考虑普朗克，并且需要收缩到与 BAO 和 LSS 数据兼容的较小区域。与上一节的运行相比，还有其他显著差异。

首先，当我们使用“小”数据集时，我们拥有普朗克高 TT 似然函数的所有干扰参数，这些参数具有强烈的非高斯后验并且彼此相关。这意味着最佳跳跃因子与 2.4 显著不同（实际上更接近 1.9）。这也意味着由于参数数量众多、相对于某些参数的似然函数的非高斯性以及我们在参数空间中仅具有最佳拟合点的不良近似（我们仅从普朗克数据的早期运行的链中提取的近似最佳拟合点计算 Fisher 矩阵），Fisher 矩阵计算很困难。对于扩展数据集，代码实际上无法使用完整的普朗克 TTTEEE + BAO + LSS 数据获得可逆的 Fisher 矩阵，因此我们不得不切换到普朗克 – 精简 TTTEEE 似然函数以摆脱干扰参数。
使用小数据集的运行显示了自动跳跃因子更新的影响：使用超级更新，代码迅速将跳跃因子调整到约 1.9，而使用更新时它保持在 2.4，导致接受率较低。表 2 显示使用超级更新时（R−1）提高了两倍。然而，这次运行也表明，在当前数据和许多非高斯参数的情况下，使用逆 Fisher 矩阵并不总是一个好主意，因为 Fisher 矩阵可能是似然函数的如此差的近似（特别是在非高斯干扰参数的方向上），以至于它实际上是比从先前仅普朗克运行的链推断出的协方差矩阵更差的输入协方差矩阵。因此，初步的 Fisher 计算使性能比超级更新 + “普朗克 2015 协方差矩阵”差三倍。
使用大数据集（但使用普朗克 – 精简 TTTEEE 似然函数）的运行显示出与预测相同的趋势：超级更新和 Fisher 都带来了显著的改进，（R−1）提高了多达两倍。
这些不同的情况使我们得出以下结论，这与我们进行的其他一些未在此处包含的测试相符：

• 使用超级更新基本上总是一个好主意。唯一可以考虑坚持使用更新的情况是最简单的情况，即当新运行涉及与先前运行非常相似的数据集和模型，以至于我们已经对协方差矩阵和最佳跳跃因子有了出色的了解。在那种情况下，更新和超级更新几乎等效，但在最不幸的情况下，超级更新可能会有一个过渡阶段，在此期间跳跃因子会偏离最佳值，然后渐近地回到最佳值，并且会在效率上略有损失。这通常是微不足道的，我们可以安全地建议在所有情况下使用超级更新：然后，根据运行的“难度”，改进范围将从可忽略不计到很大。
• 当存在许多非高斯参数，例如普朗克干扰参数时，Fisher 矩阵计算经常失败，即使它没有失败，逆 Fisher 矩阵通常也是提议密度的不良近似，与从具有相同干扰参数的链推断出的任何输入协方差矩阵相比。例如，这意味着当使用完整的普朗克高 似然函数时，应该使用分布式的“普朗克 2015 协方差矩阵”，或者自己从先前运行中获得的协方差矩阵，而不是 Fisher 选项。在几乎所有其他情况下，我们发现计算并从逆 Fisher 矩阵开始是加速收敛的非常强大的方法。

​

综上所述，蒙特卡罗 Python 通过超级更新策略和 Fisher 矩阵计算的引入，在参数推断方面取得了显著的进展。超级更新在处理不同复杂度模型和数据集时展现出良好的适应性，能够有效提升收敛速度并减少计算资源的消耗。尤其是在面对先验知识匮乏的情况，如在处理模拟数据与实际数据结合的复杂场景中，其优势更为突出。

然而，Fisher 矩阵计算虽然在多数 MCMC 预测场景中表现卓越，但在面对大量非高斯参数以及真实数据且似然形状复杂的情况下，可能会面临挑战，甚至可能导致性能下降。因此，在实际应用中，需要根据具体的数据特征和模型需求，灵活选择是否采用 Fisher 矩阵计算以及如何结合超级更新策略，以达到最佳的参数推断效果。

未来，蒙特卡罗 Python 的发展可以进一步聚焦于提升 Fisher 矩阵计算在复杂情况下的稳定性和准确性，例如优化最小化算法以更好地处理非高斯似然函数和众多干扰参数的情况。同时，对于超级更新策略，可以探索更加智能的跳跃因子调整机制，进一步提高其在各种场景下的效率，从而为研究提供更加强有力的工具支持，助力科学家们更深入地探索宇宙的奥秘，揭示参数背后隐藏的物理规律。

参考文献

[1] B. Audren, J. Lesgourgues, K. Benabed, S. Prunet, Conservative Constraints on Early Cosmology: an illustration of the Monte Python cosmological parameter inference code, JCAP 1302 (2013) 001.

[2] A. Lewis, S. Bridle, Cosmological parameters from CMB and other data: A Monte Carlo approach, Phys. Rev. D66 (2002) 103511. arXiv:astro-ph/0205436, doi:10.1103/PhysRevD.66.103511.

[3] A. Lewis, Efficient sampling of fast and slow cosmological parameters, Phys. Rev. D87 (10) (2013) 103529.

[4] J. Zuntz, M. Paterno, E. Jennings, D. Rudd, A. Manzotti, S. Dodelson, S. Bridle, S. Sehrish, J. Kowalkowski, CosmoSIS: Modular Cosmological Parameter Estimation, Astron. Comput. 12 (2015) 45–59. arXiv:1409.3409, doi:10.1016/j.ascom.2015.05.005.

[5] J. Lesgourgues, The Cosmic Linear Anisotropy Solving System (CLASS) I: Overview. arXiv:1104.2932.

[6] D. Blas, J. Lesgourgues, T. Tram, The Cosmic Linear Anisotropy Solving System (CLASS) II: Approximation schemes, JCAP 1107 (2011) 034.

[7] J. Lesgourgues, The Cosmic Linear Anisotropy Solving System (CLASS) III: Comparision with CAMB for ΛCDM. arXiv:1104.2934.

[8] J. Lesgourgues, T. Tram, The Cosmic Linear Anisotropy Solving System (CLASS) IV: efficient implementation of non-cold relics, JCAP 1109 (2011) 032.

[9] M. Zumalacarregui, E. Bellini, I. Sawicki, J. Lesgourgues, hi´class: Horndeski in the Cosmic Linear Anisotropy Solving System. arXiv:1605.06102.

[10] G. W. Pettinari, C. Fidler, R. Crittenden, K. Koyama, D. Wands, The intrinsic bispectrum of the cosmic microwave background, J. Cosmology Astropart. Phys. 4 (2013) 3. arXiv:1302.0832, doi:10.1088/1475-7516/2013/04/003.

[11] P. Stocker, M. Krämer, J. Lesgourgues, V. Poulin, Exotic energy injection with ExoCLASS: Application to the Higgs portal model and evaporating black holes. arXiv:1801.01871.

[12] B. Bolliet, B. Comis, E. Komatsu, J. F. Mac´ías-Perez, Dark Energy from the Thermal Sunyaev Zeldovich Power Spectrum. arXiv:1712.00788.

[13] E. Di Dio, F. Montanari, J. Lesgourgues, R. Durrer, The CLASSgal code for Relativistic Cosmological Large Scale Structure, JCAP 1311 (2013) 044. arXiv:1307.1459, doi:10.1088/1475-7516/2013/11/044.

[14] Y. Dirian, S. Foffa, M. Kunz, M. Maggiore, V. Pettorino, Non-local gravity and comparison with observational datasets. II. Updated results and Bayesian model comparison with ΛCDM, JCAP 1605 (05) (2016) 068. arXiv:1602.03558, doi:10.1088/1475-7516/2016/05/068.

[15] A. Lewis, A. Challinor, A. Lasenby, Efficient computation of CMB anisotropies in closed FRW models, Astrophys. J. 538 (2000) 473–476. arXiv:astro-ph/9911177, doi:10.1086/309179.

[16] C. Howlett, A. Lewis, A. Hall, A. Challinor, CMB power spectrum parameter degeneracies in the era of precision cosmology, JCAP 1204 (2012) 027. arXiv:1201.3654, doi:10.1088/1475-7516/2012/04/027.

[17] A. Refregier, L. Gamper, A. Amara, L. Heisenberg, PyCosmo: An Integrated Cosmological Boltzmann Solver. arXiv:1708.05177.

[18] J. Dunkley, M. Bucher, P. G. Ferreira, K. Moodley, C. Skordis, Fast and reliable MCMC for cosmological parameter estimation, Mon. Not. Roy. Astron. Soc. 356 (2005) 925–936. arXiv:astro-ph/0405462, doi:10.1111/j.1365-2966.2004.08464.x. ​