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則通過模擬的方法對高維積分進行計算,消除掉積分困難問題。
總結
以上是生活随笔為你收集整理的MCMC和贝叶斯统计在宇宙微波背景辐射(CMB)中应用的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: Android 福彩3D体彩排列(源码+
- 下一篇: 刻录系统盘实战学习