宇宙学中的统计方法【乙】

· by · in · [ Co-Reading ] ·

宇宙学中的统计方法笔记系列，第二篇：CMB 中的统计方法。

系列文章目录：

继续前面的 likelihood 的讨论，不过现在要扩展到 CMB 中来。与前面的不同的地方是，之前的例子中，我们的测量是每次一个值的，但是在 CMB 的测量中，我们的数据是一个分布，而且之前的例子的每个数据点我们是有理论可以完全预测的，而在 CMB 中，我们的分布的一个特定的数据点，我们是没有理论来完整预测的，而是只有一个对所有点的统计的预测。

CMB 的数据记录和相关理论

CMB 测量方法和数据记录待补充。flag 1

CMB 的观测中，我们把天空分成很多像素点，针对每个像素点记录温度各项异性。下面的分析中我们使用的温度各项异性，是相对于平均温度的差异。

对于各项异性的理论预言，我们选择较为简单的 Gaussian 分布的情况，即每个像素点上的温度差异分布符合 Gaussian 分布。

单像素点的 likelihood

现在我们先建立数学描述。对于天空上的一个点，我们可以测量温度各项异性，其值为 \(s\)，而预测值（estimator）记做 \(\Delta\)。跟前面的例子一样，我们同样会有一个代表预测结果可靠性的量（ Modern Cosmology 中把这个叫做 variance of the estimator）的量 \(C_N\)，类似前面的例子，对于 \(s\) 和 \(\Delta\) 同样有一个分布，
\[ P[s] = \frac{1}{\sqrt{ 2\pi C_N }} \exp\left( -\frac{(\Delta – s)^2}{2C_N}\right) \]
对于理论预测，如果我们取暴涨模型所给定的 Gaussian 分布，那么就有下面的这个关于 \(S\) 的分布，其中参量有 \(s\) 和 \(C_S\)
\[ P[\theta(\Delta,s)] = \frac{1}{\sqrt{2\pi C_S}}\exp\left( -\frac{s^2}{2C_S}\right) \]

现在，我们关心的两个量是 \(\Delta\) 和 \(C_S\)，因为 \(\Delta\) 是我们要的温度各项异性值，而 \(C_S\) 直接跟我们后面要了解的一个量 \(C_l\) 关联，\(C_l\) 是我们要的最终结果，而 \(s\) 这样的量我们并不关心。所以，现在这种情况下，我们需要求的 likelihood 其实是 \(P[\Delta|C_S]\)，也就是给定一个带有参数 \(C_l\) 的理论，得到数据 \(\Delta\) 的概率。
那么根据概率论（之前的例子中有具体计算步骤）
\[ P[\Delta|C_S] = \sum_s P[\Delta|s] \times P[s|C_S] \]

在 Modern Cosmology 中，Scott Dodelson 并没有对这里的 likelihood 进行很多解释。这里 flag 一下，提醒以后注意回来检查精简这段描述。flag 2

于是 Likelihood function 是概率密度的积分
\[ {\scr L} = \int^\infty_{-\infty} \mathrm ds \frac{1}{\sqrt{2\pi C_S}}\exp\left( -\frac{s^2}{2C_S}\right) \frac{1}{\sqrt{ 2\pi C_N }} \exp\left( -\frac{(\Delta – s)^2}{2C_N}\right) \]

定义 \(C=C_S+C_N\) 是 full covariance matrix，把 \(s\) 积分积掉，

\begin{eqnarray}
{\scr L} &=& \int \mathrm ds \frac{1}{2\pi \sqrt{C_S C_N}} \exp\left( -\frac{s^2 C_N + (\Delta – s)^2C_S}{2C_S C_N} \right) \\
&=& \frac{1}{2\pi\sqrt{C_S C_N}} \int \mathrm ds \exp\left( -\frac{ C(s- C_S \Delta/C)^2 + \Delta^2 (C_S C_N/C) }{2C_N C_S} \right) \\
&=& \frac{1}{\sqrt{2\pi C_S}}\exp\left( – \frac{\Delta^2}{2C} \right) \int \mathrm d(s – C_S\Delta/C)\exp\left( – \frac{C(s – C_S\Delta/C)^2 }{2C_S C_N} \right) \\
&=& \frac{1}{\sqrt{2\pi C}}\exp\left( -\frac{\Delta^2}{2C} \right)
\end{eqnarray}

多像素点情况

待补充。详细说明为何可以使用矩阵形式。flag 3

这样我们就得到了天空中一个像素点的情况，可以推广到 N 个像素点的情况，
\begin{equation}
{\scr L} = \frac{1}{(2\pi)^{N/2}(\mathrm {det} \bf C)^{1/2}} \exp\left( – \frac{1}{2} \bf \Delta \bf C^{-1} \bf\Delta \right) \label{eqn-multipixellikelihood}
\end{equation}

上式（方程(\ref{eqn-multipixellikelihood})）中 \(\Delta\) 是长度为 \( N \)测量数据矢量，而 \(C\) 是 \(N\times N\) 的 full corariance matrix。

Covariance Matrix

上面出现的三个 covariance matrix 分别是 \(C_S\)，\(C_N\)，\(C\)。其中 \(C_S\) 是出现在 \(s\) 的分布中的，也就是说是理论上 \(s\) 的分布宽度，称之为 signal covariance matrix，与测量中的噪声无关。从 \(C_S\) 到 \(C\) 其实是个很复杂的过程，严格的说，\(C\) 中需要考虑到任意的两个 parameter 直接的关联。

详细说明如何关联。flag 4

Covariance Matrix 的求解

如果我们之考虑 covariance matrix 的对角项，
\[ C_{S,ii} \equiv \langle s_i s_i \rangle \]
其中元素
\[ s_i = \int \mathrm d \hat n\Theta(\hat n) B_i(\hat n) \]
其中 \(\Theta(\hat n)\) 是温度背景，\(B_i\) 是 beam pattern，也就是代表信号的性质的量。对于每个像素点都要把所有的方向的光子都考虑进来。
把 \(s_i\) 带入 \(C_{S,ii}\) 可以得到 covariance matrix

\[ \frac{C_{S,ii}}{T^2} = \int \mathrm d\hat n \int \mathrm d\hat n’ B_i(\hat n)B_i(\hat n’) \sum_{lm} Y_{lm}(\hat n)\sum_{l’m’} Y_{l’m’}^*(\hat n’)\langle a_{lm}a_{l’m’}^* \rangle \]

其中用到了
\[ \Theta(\hat n) = \sum_{i=1}^\infty \sum_{m=-l}^{l} a_{lm} T Y_{lm}(\hat n) \]

*Modern Cosmology* 中写 \(\Theta\) 时不含 \(T\)，但是书中 11.35 式却含有 \(T\)，所以我把 \(T\) 添加到了 \(\Theta\) 的分解中。 flag 5

把下式带入 covariance matrix 的表达式
\[ \langle a_{lm} \rangle; \langle a_{lm} a_{l’m’}^* \rangle = \delta_{ll’}\delta_{mm’}C_l \]
把 \(l’,m’\) 的求和求出来
\[ \frac{C_{S,ii}}{T^2} = \int \mathrm d\hat n \int \mathrm d\hat n’ B_i(\hat n) B_i(\hat n’) \sum_{l} C_l \sum_m Y_{lm}(\hat n) Y_{lm}^*(\hat n’) \]

我们知道球谐函数有如下性质
\[ \sum_m Y_{lm}(\hat n)Y_{lm}^*(\hat n’) = (2l+1) P_l (\hat n\cdot \hat n’)/4\pi \]
定义
\[ W_{l,ii} = \int \mathrm d\hat n\int \mathrm d\hat n’ B_i(\hat n)B_i(\hat n’) P_l(\hat n\cdot \hat n’) \]
其中的 \(\hat n \cdot \hat n’\) 可以写成 \(\cos(x)\)，其中 \(x\) 是 \(\hat n\) 和 \(\hat n’\) 之间所夹的弧长。这样我们可以在一个二维平面上定义两个矢量，\(\vec x\) 和 \(\vec x’\)，使得 \( \cos(x) = \cos(|\vec x – \vec x’|) \)，即 \(\hat n\cdot \hat n’ = \cos(\vec x – \vec x’)\)
于是
\[ W_{l,ii} = \int \mathrm d^2x\int \mathrm d^2 x’ B_i(x) B_i(x’) P_l(\cos(|\vec x – \vec x’|) ) \]

为什么变成二次导数了？我求的是一次导数。回来看看。下面还是继续按照 *Modern Cosmology* 里面所讲的来写。 flag 6

所以 covariance matrix 变成
\[ \frac{C_{S,ii}}{T^2} = \sum_l \frac{2l+1}{4\pi} C_l W_{l,ii} \]

关于上面提到的 window function \(W_{l,ii}\)， 还可以继续化简。首先，我们把 Legendre polynomials 表达出来，不过这里要用一个近似，就是在 l 很大的时候， Legendre polynomials 近似变成 Bessel function.
\[ P_{l} (\cos(|\vec x – \vec x’|)) \rightarrow J_0 (l|\vec x – \vec x’|) = \frac{1}{2\pi} \int_0^{2\pi} \mathrm d\phi \exp{-il |\vec x – \vec x’| \cos\phi} \]
其中 \(\phi\) 的意思是 \(l\) 和 \(\vec x – \vec x’\) 之间的额夹角。\(l\) 的方向可以选做某条坐标轴的方向比较方便。
这样我们就可以把 window function 简化成
\[ W_{l,ii} = \frac{1}{2\pi} \int_0^{2\pi}\mathrm d\phi |\tilde B_i(\vec l)|^2 \]
其中
\[ \tilde B_i(\vec l) \equiv \int \mathrm d^2x B_i(\vec x)e^{-i \vec l \vec x} \]
是 \(B_i\) 的 Fourier 变换。

这样我们就可以把 covariance matrix 的对角项表示出来了，选择特定的 window function，就可以求解。

但是这只是对角项，实际上非对角项也是存在的，严格的计算需要把非对角项也包括进来。

Window Function 的例子

Gaussian Beam

对于 Guassian beam 的情况，beam pattern 写作
\[ B_i(\vec x) = \frac{1}{2\pi \sigma} \exp\left( -\frac{(\vec x -\vec x_i)^2}{2\sigma^2} \right) \]
先令 \(\vec x_i = 0\)，方便 Fourier 变换。
Fourier 变换得到
\[ \tilde B_i(\vec l) =\frac{1}{2\pi \sigma^2} \int \mathrm d^2x e^{-i\vec l \cdot\vec x} \exp\left( -\frac{x^2}{2\sigma^2} \right) = e^{-l^2\sigma^2/2} \]
因为 Fourier 变换之后并不涉及 \(\vec l\) 方向，所以 window function 可以立刻得出来
\[ W_{l,ii} = e^{-l^2\sigma^2} \]

这个函数的行为我们都很熟悉，如下所示，

也就是说对于很大的 \(l\)，window function 会消失。换句话说，因为很大 \(l\) 正好对应很小的角分辨率，所以就是说，对于那些角分辨率小于 beam 的宽度的时候，这部分对于 covariance matrix 的贡献就自动被 window function 过滤掉了。

下面是一个例子。左上图是实空间里面的 beam 的 pattern，右上图表示 beam 做 Fourier 展开后在 Fourier 空间中的 pattern，下部图表示 window function。

Differencing a Gaussian Beam

Beam function 也可以选择如下形式
\[ B(x,y) = \delta(y)[\delta(x-x_0) – \delta(x + x_0)] \]
这样的 beam function 意思是？

做 Fourier 变换，得到
\[ \tilde B (\vec l) = 2i\sin(l_x x_0) \]

带入 window function 得到
\begin{eqnarray}
W_l &=& \frac{1}{2\pi} 4 \int _0^{2\pi} \mathrm d\phi \sin^2(lx_0\cos\phi) \\
&=& \frac{1}{\pi}\int_0^{2\pi} \mathrm d\phi (1-\cos(2lx_0\cos\phi)) \\
&=& 2(1-P_l[\cos(2x_0)])
\end{eqnarray}
这里又用了前面用过的 Bessel function 和 Legendre function 之间的近似。

那么如果 beam 的宽度不是无穷小，而是有限大的，那么前面的 \(\delta\) 函数就变化成积分就好了，也就是说
\[ B(x,y) = \frac{ 2\pi \sigma^2 } \int \mathrm dx’\mathrm dy’\exp\left( -\frac{(x-x’)^2 + (y -y’)^2}{2\sigma} \right) \times \delta(y’) [\delta(x’ – x_0) – \delta(x’ + x_0)] \]

那么相应的 window function 变换就是
\[ W_l = e^{-l^2\sigma^2} (1 – P_l[\cos(2x_0)]) \]

下面是一个例子。左上图是实空间里面的 beam 的 pattern，右上图表示 beam 做 Fourier 展开后在 Fourier 空间中的 pattern，下部图表示 window function。

下图是一个更加紧致的 beam 的情况。

本文共有 6 个 flag。请仔细检阅并消除所有 flag。

长程相关:

  1. 宇宙学中的统计方法【甲】
  2. 宇宙学中的统计方法【丙】
  3. 宇宙学中的统计方法【丁】
  4. 我们的宇宙可能由闭宇宙演化为开宇宙么？[第二节]
Tags: , ,