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

· by · in · [ Co-Reading ] ·

本文是宇宙学统计中的非常非常小的一部分。是我刚刚读了一篇文章(arXiv:0911.3105)的一部分后所学到的。如有谬误，欢迎指出。2012-11-30补充了一个 Modern Cosmology 中的实例。此为系列文章第一篇：统计基础。

系列文章目录：

约定

假定我们的宇宙遵从特定的规律，可以设想这样的两种观测者。一种是自然之上的创世者们，比如叫他们宇宙定制公司，他们非常清楚宇宙所遵从的规律，而且也知道宇宙是在什么样的参数下运行；另一种是宇宙之内卑微的生存者，比如我们地球上这些渺小的人类，我们并不知道宇宙是精确的规律是什么，也不知道在某个规律下面宇宙的参数是如何设定的，我们只能实验和观测，更加难堪的是，在宇宙学中，我们基本上只能进行观测，想要尝试去改变一些宇宙的设定来得到某些结果，目前还只能出现在科幻小说之中。

那么对于宇宙定制公司的人们来说，他知道宇宙的精确的理论是 H ，（对于我们来说是个假设而已，因为从我们的角度来看就是个 Hypothesis，简称 H ），而我们生活的宇宙可以用一个模型 \(\theta(\theta_1,\theta_2,…)\) 来描述，其中 \(\theta_i\) 为模型参数，而且他们知道宇宙运行之后的数据 D，（比较幸运的是，对于很多数据 D，我们也是可以得到的，因为我们就限定 D 为我们可以得到的数据）。

基本概念

那么从宇宙定制公司的角度来看：

可以根据 H 来计算 \(\theta\) 模型下的某个数据出现的理论概率，比如现在在一个范围 A 内，共有 M 颗1型超新星，和 N 颗2型超新星。他们可以计算在拿到了 m 颗1型超新星和 n 颗2型超新星之后，随便抽取一颗是1型的概率，因为整个宇宙在任何时空点的配置他们都清楚的很。这种我们统统称之为计算 \(P(D|\theta(\theta_1,\theta_2,..))\)，即 likelihood. 意思是说，我们现在知道一个模型\(\theta\)和参数\(\theta_1,\theta_2,…\)，得到某个结果 D 的概率是多少。

这差不多是宇宙定制公司用的最多的计算了吧。

而对于我们来说，情况变得就复杂多了。我们只有观测数据，并不知道宇宙的理论 H 和具体的模型及参数 \(\theta(\theta_1,\theta_2,…)\)。而我们希望通过这些观测了解到宇宙背后的理论和具体模型及参数。
因此从我们的角度来看却只能计算 Posterior，记作 \(P(\theta|D)\)，意思是说，我们有一堆观测数据 D ，那么一个特定的模型\(\theta\)并且取某些特定的参数 \(\theta(\theta_1,\theta_2,…)\)的概率是多少。

可以看到这样一种记法，| 是符合条件符合，可以直接读作”符合条件”，因为这个条件概率的条件放在后面。而 \(P(x | y)\) 的意思就是，给定 \(y\) 发生 \(x\) 的概率。

为了更明白，我重复 arXiv:0911.3105 中的一个例子。

If you have an urn with N red balls and M blue balls and you draw one ball at the time then probability theory can tell you what are your chances of picking a red ball given that you have already drawn n red and m blue: P(D|H). However this is not what you want to do: you want to make a few drawn from the urn and use probability theory to tell you what is the red vs blue distribution inside the urn: P(H|D). In the Frequentist approach all you can compute is P(D|H).

这里 \(P(D|H)\) 是 likelihood， 而 \(P(H|D)\) 是 posterior. 可以明显的看到，likelihood 是指我们明了系统的详细配置的情况下计算某件事情发生的概率，这个很强大;而 posterior 是我们不知道系统的配置的情况下，可以根据很多的实验来计算出的一个概率。这个例子中提到，如果我们可以做无穷次实验，那么就可以计算 posterior，即使我们是生活在2012年的地球人。==~

Likelihood 的例子

从 Modern Cosmology 中（337页）取一个例子，稍加修改。算是一段插曲。

我们要测量自己的体重，作为科学工作者，大家肯定都知道要多次测量，去平均值，还要给出误差。那么平均值该如何求的呢？误差又该如何求的呢？

我们现在可以假定一个理论，说每次测量的体重有两部分组成，一个是体重的精确值 \(theta_0\)，另一个是噪音（noise）。如果我们假定测量是个 Gaussian distribution，
\[ P[D|\theta(\theta_0,\sigma)] \equiv \scr L(D,\theta(\theta_0,\sigma)) = \frac {1}{ \sqrt{2\pi\sigma^2}} exp\left({-\frac{(\theta – \theta_0)^2}{2\sigma^2}}\right) \]
那么精确值 \(\theta_0\) 是分布的极值点，而每次测量的噪音的大小可以用分布中的 \(\sigma\) 来表征。这样我们这个理论中有两个参数：\(\theta_0\)， \(\sigma\)。

上面的 \(P[D|\theta(\theta_0,\sigma)]\) 是指的一次测量中，给定一个理论 \(\theta(\theta_0,\sigma)\)，我们测得一个数据 \(D\) 的概率。但是实际上我们可能进行 N 次独立测量，这样的话，likelihood 应该是
\[ \scr L(D,\theta(\theta_0,\sigma)) = \frac {1}{ (2\pi\sigma^2)^{N/2}} exp\left({-\frac{\sum^{N}_{i=1}(\theta_i – \theta_0)^2}{2\sigma^2}}\right) \]

但是，我们做实验的人并不想要上面所求的 \( P[D|\theta(\theta_0,\sigma)]\)，相反，我们要的是 \(P[\theta(\theta_0,\sigma)|D]\)。

概率论中，我们有
\[ P[A\cap B] = P[A|B]P[B] = P[B|A]P[A] \]
这样我们就可以得到 likelihood 跟我们要求的量之间的关系。
\[ P[\theta(\theta_0,\sigma)|D] = \frac{ P[D|\theta(\theta_0,\sigma)P[\theta(\theta_0,\sigma)]] }{ P[D] } \]

式中分母 \(P[D]\) 并不会影响我们的概率极值在 \(\theta_0 ~ \sigma\) 相空间的位置，所以我们暂时把这一项丢掉，剩下的分子上的 prior，即 \(P[\theta(\theta_0,\sigma)]\) 项，我们可以假设是均匀分布的。这样我们就可以简化到下式，
\[ P[\theta(\theta_0,\sigma)] \sim \scr L \]

这样我们就把要求解的量 \(P[\theta(\theta_0,\sigma)]\) 与 likelihood 联系起来了，而 likelihood 我们是可以写出表达式的。

那么下面我们就是要找到参数空间的最佳拟合点，也就是参数空间中概率最大的地方。分别对两个参数求导，\(\frac{\partial \scr L}{\partial {\theta_0}}\)，\(\frac{\partial \scr L}{\partial \sigma}\)。

高斯分布的情况下，情况就简单了。
\[ \frac{\partial \scr L}{\partial {\theta_0}} = \frac{\sum^N_{i=1} (\theta_i – \theta_0) }{ \sigma^2 (2\pi\sigma^2)^{N/2} } \exp \left( -\frac{\sum^N_{i=1} (\theta_i – \theta_0)}{2\sigma^2} \right)\]
令上式为零，得到
\[ \sum^N_{i=1}(\theta_i – \theta_0) = 0 \]
也就是说，\(\theta_0\) 是 \(\theta_i\) 的平均值，记做 \(\bar\theta\)。

这段的意思是说，如果我们做了一堆实验，那么实验的理论值（按照我们一开始提出的理论）应该是所有值的代数平均 \(\bar\theta\)。

对于另外一个参数 \(\sigma\)，

\[ \frac{\partial \scr L}{\partial \sigma} = \scr L \left( -\frac{N}{2\sigma^2} + \frac{\sum^N_{i=1} (\theta_i – \theta_0)^2 }{2\sigma^4} \right) \]

上式为零，得到
\[ \sigma^2 = \frac{\sum^N_{i=1} (\theta_i – \theta_0) }{ N } \]
也就是说， \(\sigma^2\) 是测量数据统计的方差。

这只是找到了参数空间的最佳拟合，这样的两个数据有时候并不够用，我们需要进一步写出更加方便的表达测量误差的方法，一个比较好的方法是看分布的半高宽，但是我们的表达式里面半高宽并不容易看出来，于是我们可以把e指数里面的求和式写开来，
\[ \exp\left( – \frac{(\theta_N – \frac 1 2 \theta_0)}{\sigma^{2/N}} – \cdots \right) = \exp\left( – \frac{(\theta_N – \theta_0)}{(2\sigma)^{2/N}}\right) \cdots \exp\left( – \frac{(\theta_1 – \theta_0)}{(2\sigma)^{2/N}}\right) \]

这样我们可以把每个因子的e指数里面的二次项系数合起来，\(N/\sigma^2\)，可以认为整个e指数的宽度是 \( \frac{\sigma}{N^{1/2}} \)。当然这只是一个大致的说法，详细的定义可以看概率论的书（？）或者看 Modern Cosmology P339。
这个结果很合理，分布的越散，测量的可信度就越小，测量次数越多，测量的可信度就越大。

到此我们把理论值和误差都表达出来了，也就是说，我们的整个 likehood 的可以由两个两个参数，分别是 estimator \(\bar\theta\) 和 variance of estimator \(\frac{\sigma}{N^{1/2}}\)，来确定。那么，现在我们就可以放弃之前带有复杂求和式的表达形式，而重新定为一个更加简单的形式：
\[ \scr L = \frac{1}{\sqrt {2\pi C_N}} \exp\left( -\frac{(\theta – \bar\theta)}{2C_N} \right) \]
其中用到 \(\sigma^2 = \frac{\sum^N_{i=1} (\theta_i – \bar\theta)}{N}\)。而 \(C_N = \frac{\sigma^2}{N}\)。

\(\chi^2\) fit / Chisquare fit

现在宇宙学积累了很多的数据，而且现在的很多数据也都是很精确的。同时我们有很多的理论模型，各式各样的，很科幻的都有。那么我们如何把数据和理论结合起来呢？这就需要用到拟合，fit.

一个常用的技巧是 least square 拟合。这里我们用到的是 chisquare 拟合。

列出我们现在所拥有的:

  • 观测数据。这里我们使用1a超新星的例子。我们有这样一堆数据：超新星名字，超新星红移 \(zo
    \)，超新星亮度/distance modulus \(do\)，超新星的distance modulus error \(e\).

  • 理论模型。这里使用最简单的 LCDM 作为例子。理论上可以计算红移对应的距离 \(dt(z,H0,\Omega m0,\Omega d0,\Omega k0) = c\cdot (1+z)\int \frac{1}{H0\sqrt{\Omega m0 (1+zz)^3 + \Omega d0 + \Omega k0(1+zz)^2}}\mathrm d zz\)
    实际上 \(dt\) 并没有那么多独立参量。比如我们可以这样建立各个参量之间的关系：
    \begin{eqnarray}
    \Omega m0 + \Omega d0 +\Omega k0 &=& 1 \\
    \end{eqnarray}
    如果是一个平坦的宇宙，那么我们还可以把 \(\Omega k0\) 设成 0 . 为了适应不同的情况，我们就抽象的写成 \(dt(z,x1,x2,…)\). \(z\)是变量。

那么我们可以定义这样一个量

\begin{equation}
\chi^2(x1,x2) = \sum_{l} \frac{(do[i] – dt(z[i],x1,x2,…))^2}{\sigma_i^2}
\end{equation}

其中 \(\sigma_i\) 是第 i 个数据点的测量误差。

一旦这样定义，我们可以看到，如果我们的理论和数据完全符合，那么 \(\chi^2(x1,x2,…)=0\). 但是，如果这样，很有可能发生了不正常的事情，实际观测和实验不太可能做到这样。但是很明显，如果 \(\chi^2\) 太小，那么也不太正常，因为这很可能意味着对观测值的误差估计太大了等等。

那么多大的时候才比较好呢？如果观测有 n 个数据点，模型有 m 个参数，定义一个量 N = m – n. 当 \(\chi^2/N \sim 1\) 时，是比较正常的，因为这正好意味着数据点的 error bar 差不多都与理论曲线相交。

这样我们可以去遍历所有的数据点或者 steep 方法或者 MCMC 方法得到 \(\chi^2\) 的最小值 \(\chi^2_{min}\)，此时对于的参数为 \(x1_m,x2_m,…\). 这组参数就是我们这个模型在这组数据点下面的 best fit.

然后我们可以去绘制置信度的 contour . 即绘制在 \(\chi^2_{min}\)，周围差别 \(\Delta\chi^2\) 的地方的曲线。对于不同的 N 值的 best fit 和不同的置信度，有不同的 \(\Delta\chi^2\) 值。

这些详细的技巧比如 covariance matrix 等等会在后续的 post 中出现。现在不再赘述。

我有一个关于 Chisquare fit 的 git 项目。这个所用到的理论基本上是中学生就可以理解的。但是到2012年6月25日为止我依然感觉结果有些问题，虽然方法是正确的。

另外，如果你喜欢在豆瓣九点等这种地方读这篇文章，可以安装 math anywhere 插件，这样公式就可以显示了：
https://chrome.google.com/webstore/detail/gebhifiddmaaeecbaiemfpejghjdjmhc

长程相关:

  1. 宇宙学中的统计方法【乙】
  2. 宇宙学中的统计方法【丙】
  3. 宇宙学中的统计方法【丁】
  4. Hubble膨胀，Galilean-invariant，non-equilibrium
Tags: , ,