MCMC和贝叶斯统计在宇宙微波背景辐射（CMB）中应用

此文为笔者笔记，描述的是MCMC（马尔科夫链蒙特卡洛）方法和贝叶斯统计在CMB中的应用，描述的只是笔者的观点，有待修改（更关键的是，笔者主要描述的是CosmoMC程序相关的东西）。

进入正题

假设待测宇宙学参数只有六个，那么描述的则是六维空间。在参数ini文件中（params_CMB_defults.ini），每个参数都有一个事先估计的范围和最可几的值(路径为：~/programs/CosmoMC/batch2> vim params_CMB_defaults.ini)

这里的参数是在运行数据拟合之前给的一个值，接着将所有这些参数对应的这个center的值放到CAMB中，此时CAMB中即有了初始值（注意这里的权限比CAMB中的权限高，一旦这里的值取定了之后CAMB中params.ini所给定的值则无效），这个时候就可以根据这些给定的宇宙学参数的值画出功率谱，并给出卡方的值（我们希望卡方越小越好，如此误差才会越小），卡方又是与似然，误差相关的，因此可以得到一系列相关的参数（见后文）。接着上图中刚才提到的六个参数又开始改变（加减步长大小），将改变之后的值又放到CAMB中，CAMB再由给定的值输出功率谱和卡方，接着对比这两个卡方的大小，MCMC将调节参数赋值的方向，使得这个方向往卡方小的地方赋值。同样的进行第三次，第四次…….并每次都将本次的卡方和上一次的卡方的值对比以决定下次赋值的方向（即步长往哪边走）；值得注意的是，倘若本次的卡方比上次的小，那么并不能保证下次的卡方一定比本次小，因为这个是不能确定的，每走一步都是试探性的，但总体趋势是往卡方小的方向走。
假设我们只有六个宇宙学参数，这六个宇宙学参数是要通过数据来确定的。那么我们第一次赋值是赋在center这个值上，以omegabh2为例，则是0.0221。这六个值的center值带入到CAMB中，于是相当于给定了初始值initial value, 当这些参数都被赋值了以后，CAMB就能够生成TT谱，BB谱，EE谱，TE谱等。接着通过以下公式计算卡方：

以上卡方只是用了TT谱，同样可以在EE谱中使用，另外，参考龚云贵的《宇宙学基本原理》P44,当我们将先验分布看成高斯型，似然函数表达成如下形式（现代宇宙学P286）：

则卡方还有以下定义：

第二个式子是在参数很多且每个参数的方差不相等的情况下使用
从上面其实可以近似将卡方和似然函数写成如下形式 ：

卡方越大（在方差sigma不变的情况下），则似然（即概率分布）越小，即实验值与理论值相差越大，这是我们不希望的；我们想要的是似然大的（即概率大的），因此需要找到卡方最小时对应的参数。另外，误差的宽度是与似然的宽度成正比的，
（似然）正态分布（高斯分布）图如下：

当方差sigma2越小，则正太分布越“高瘦”，那么置信度越大，即越有可能相信理论和实验符合得很好；因此，实验上现在主要是降低卡方的分母（即方差）的值来得到理想的结果。以下分步骤进行讲解：
1，六个宇宙学参数给定初值（.ini文件）

2，科学上而言，我们只给出卡方最小的时候对应的参数值是不够有说服力的，也不具备科学条件，我们还要给出取这些值时的误差，现在我们就来分析这个误差如何计算
首先我们已经得到了卡方的值，在此强调一句，对于多次测量，卡方的计算是通过协方差矩阵进行的（现代宇宙学P287）（单个就是一维矩阵）。
假设我们将这六个参数取了n=600次，即得到了600个卡方，这600个卡方有一个最小值，我们称为卡方min，那么相当于我们将这些参数在卡方min的周围进行了599次取值

离卡方越近的点越密集，因为实际上如果假设是高斯型分布，那么计算机很容易就找到卡方最小的点，即很快就可以找到这六个参数使得其对应的卡方最小，即跟实验上最接近。值得一提的是，这600次实验可能在400次或者更小时就已经找到卡方最小的点，并且之后的两百次卡方会一直是这个不变的值，叫做马尔科夫链（马氏链）的平稳分布。
这些六百次的测试将得到关于这六个参数的一个分布（想象一下六维空间中每个正交的轴代表一个参数，那么这些参数共同构成一个概率分布），当我们想知道其中一个的概率分布时，如omegabh2,就需要进行五次积分，将其他五个参数积分掉，（即联合概率分布化为边缘概率分布）于是得到如下的图;

也可以只积分掉四个，留下两个，看他们之间的关系，（如ns与omegabh2）则图如下：

这里内圈代表一个sigma,外圈代表2个sigma,sigma对应的表可以在维基百科找到，如下表：

像这样的圈比较大的表示横纵坐标这两个参数相关性不大，比如1个sigma的置信范围内当横坐标取0.0224时，纵坐标可以取0.965-0.978这样一个范围内的值，即变化还是很大的，当是接近线形图的时候才表示相关性很大，如下图：

这个时候横坐标一旦取定，纵坐标可变的范围极小，有时候几乎可以认为是线性相关，这样就可以减少一个维度（图中是正相关）
那么既然我们得到了600个卡方，这些卡方也是可以画出一个分布图的，如下：

这里横坐标是事件标记（每次试验看成一个事件，比如400-600次时的卡方应该是接近同一个值，那么也就是说这样的卡方对应的值在这600百次中出现的概率最大，纵坐标表示的就是卡方出现的概率）所以在峰值最高的时候表示卡方最小（对应的纵坐标事件标记可能是n=600，其附件可能是n=590范围内）
通过画出的图便可表示误差范围

补充：
贝叶斯公式表示如下：

这里x是样本（实验数据），theta是参数（如上面的六个参数），通常情况下，我们是有一个理论参数对应的值，然后多次实验得到很多数据，再用这些数据去拟合得到这个参数；而我们的CMB恰恰相反，因为实验数据只有一组，也就是天图只有一张（即一个确定的实验功率谱），然后通过多次调节理论参数，使得这个参数和实验很好的符合，因此用的是后验分布（等式右边，在有实验数据的前提下去得到参数），等式右边分子第一项是似然函数，第二项是先验分布，分母是一个积分（积分后是一个常数）；贝叶斯统计需要先从经验上主观给定先验，这正是上世纪很多人反对贝叶斯统计的关键所在。贝叶斯统计是在给定先验分布的前提下计算出后验分布的矩（后验均值、后验方差）、后验概率密度函数等，当参数比较多的时候分母将涉及到高维度的积分，使得计算变得异常复杂。MCMC则通过模拟的方法对高维积分进行计算，消除掉积分困难问题。