陳威威, 王洪光, 張顏榮

?

基于貝葉斯推斷的脈沖星偏振位置角擬合算法

陳威威, 王洪光, 張顏榮

(廣州大學(xué)物理與電子工程學(xué)院, 廣東廣州, 510006)

脈沖星輻射束的偏振位置角與相位的關(guān)系可用旋轉(zhuǎn)矢量模型(RVM)描述, 現(xiàn)有的Levenberg-Marquardt非線性擬合和格點(diǎn)搜尋算法對(duì)RVM的擬合存在效率低、易過(guò)擬合等問(wèn)題。本文發(fā)展了一套基于貝葉斯推斷的馬爾可夫鏈蒙特卡洛算法(MCMC), 并用該方法對(duì)已有文獻(xiàn)中的10顆脈沖星偏振位置角進(jìn)行了擬合。結(jié)果表明, MCMC算法不僅能得到與已有文獻(xiàn)一致的結(jié)果, 而且該算法具有高效、不易過(guò)擬合、能更好地估計(jì)參數(shù)的可信度區(qū)間等優(yōu)點(diǎn)。

脈沖星; 輻射束; 偏振位置角; 旋轉(zhuǎn)矢量模型

脈沖星是一種快速自轉(zhuǎn)的致密天體, 具有超強(qiáng)的磁場(chǎng), 整體上是一個(gè)偶極磁場(chǎng)。在星體表面附近可能存在多級(jí)場(chǎng), 通常認(rèn)為其射電輻射來(lái)源于磁層中偶極場(chǎng)占主導(dǎo)的區(qū)域, 其旋轉(zhuǎn)矢量模型如圖1所示。偶極場(chǎng)磁軸與自轉(zhuǎn)軸之間的夾角稱為磁傾角(), 觀測(cè)者視線與自轉(zhuǎn)軸之間的夾角稱為視線角(), 而視線與磁軸之間的最小夾角通常被稱為撞擊角(), 這3個(gè)角度只有2個(gè)角獨(dú)立, 且滿足方程=+。脈沖星的輻射來(lái)源于其磁層, 在上述與輻射幾何相關(guān)的3個(gè)角度不同的情況下, 觀測(cè)者看到的現(xiàn)象會(huì)有差異, 磁層中的動(dòng)力學(xué)過(guò)程也與磁傾角密切相關(guān), 因此, 如何從觀測(cè)中確定磁傾角和撞擊角(或視線角)是脈沖星研究中的重要課題之一。1969年Radhakrishnan和Cooke[1]提出了旋轉(zhuǎn)矢量模型(Rotation Vector Model, RVM), 給出了磁傾角、撞擊角和線偏振位置角的幾何關(guān)系:

其中0和0為磁軸與自轉(zhuǎn)軸所在磁力線平面對(duì)應(yīng)的相位和線偏振位置角。線偏振位置角可以從偏振觀測(cè)數(shù)據(jù)中的Stokes參量計(jì)算出來(lái), 因此通過(guò)使用RVM, 可以擬合出、、0以及0。相比其他確定和的方法, RVM擬合被認(rèn)為是最簡(jiǎn)潔可靠的方法。目前用這類方法限定了和的脈沖星不到90顆, 較為系統(tǒng)性的有Everett等[2]、Keith等[3]、Rookyard等[4]的工作。目前主要使用2種方法來(lái)擬合和, 一種是格點(diǎn)計(jì)算, 將定義域劃分為一定數(shù)量的格點(diǎn), 通過(guò)計(jì)算各格點(diǎn)參數(shù)值下的2來(lái)尋找最優(yōu)參數(shù)值; 另一種是介于高斯—牛頓法與梯度算法之間的Levenberg-Marquardt(簡(jiǎn)稱LM)[5]非線性擬合方法, 該方法借助模型對(duì)函數(shù)進(jìn)行求導(dǎo), 尋找2的極小值點(diǎn)。

圖1 磁傾角α與撞擊角β

Rookyard等人[4]使用格點(diǎn)搜索的辦法擬合RVM, 得到了26顆脈沖星的和限定結(jié)果。這種方法面臨的問(wèn)題有: (1) 在多維情況下計(jì)算量巨大, 耗時(shí)長(zhǎng); (2) 由于沒(méi)有明顯的優(yōu)化路徑, 極大地消耗計(jì)算資源; (3) 精度受限制。

在Everett等[2]和Keith等[3]的工作中, 他們使用LM算法對(duì)磁傾角與撞擊角進(jìn)行擬合, 分別給出10顆和5顆脈沖星的磁傾角和撞擊角的最優(yōu)擬合值。在其他文獻(xiàn)中有不少也沿用這種方法[6]。LM以及其他類似算法的特點(diǎn)是優(yōu)化路徑朝著極值點(diǎn)推進(jìn), 沿著2減少的方向取參數(shù)值, 直至找到局部極值點(diǎn)。在RVM的擬合中, 使用此方法將會(huì)面臨幾個(gè)問(wèn)題: (1) 優(yōu)化路徑的選擇過(guò)度依賴2的大小, 容易導(dǎo)致過(guò)擬合; (2) 優(yōu)化路徑以及最后的極值點(diǎn)容易受初始值的影響, 最終停在局域極值點(diǎn); (3) 這種優(yōu)化方式更關(guān)注極值點(diǎn), 對(duì)置信區(qū)間的分析能力較弱; (4) 在每換一套參數(shù)值后, 都要重新對(duì)似然函數(shù)進(jìn)行求導(dǎo), 當(dāng)參數(shù)維度較高時(shí), 計(jì)算量將急劇增大, 耗時(shí)長(zhǎng)。Everett等[2]指出, 使用LM算法的主要困難在于和存在強(qiáng)相關(guān), 在偏振位置角曲線相位較窄或數(shù)據(jù)對(duì)RVM有一定程度偏離的情況下, 擬合程序有時(shí)會(huì)給出趨向于180°或者0°及趨向于0°的結(jié)果, 這就是所謂的過(guò)擬合現(xiàn)象(盡管這仍符合RVM的數(shù)學(xué)關(guān)系)。發(fā)生這種情況的時(shí)候,2在(,)趨向于(180°, 0°)或(0°, 0°)的角落緩慢地減少。這種現(xiàn)象在0與0作為自由參量的情況下尤為嚴(yán)重, 擬合程序有時(shí)鎖定在一組(0,0), 使得無(wú)法在降低χ2值的同時(shí)避免出現(xiàn)過(guò)擬合。為了減輕該問(wèn)題的影響, Everett等[2]加入了Powell[7]擬合算法, 先改進(jìn)初始值的選擇, 再把改進(jìn)了的初始值輸送到LM算法中進(jìn)行擬合。他們首先使用其他研究者給出的結(jié)果作為初始值, 再手動(dòng)微調(diào)0與0來(lái)擬合數(shù)據(jù), 最后使用擬合程序?qū)?、?和0一起擬合。過(guò)擬合問(wèn)題的存在不僅增加了擬合的繁瑣程度, 不利于大樣本的自動(dòng)化擬合, 而且也限制了此方法所能處理的脈沖星樣本。

鑒于以上2種方法存在的問(wèn)題與不足, 有必要發(fā)展新的擬合和參數(shù)的方法。馬爾科夫鏈蒙特卡洛(MCMC)方法是一種基于貝葉斯推斷, 能適用于多參數(shù)擬合的有效方法, 在宇宙學(xué)等許多天文領(lǐng)域都有應(yīng)用, 但在限定脈沖星輻射幾何的研究中還未見(jiàn)應(yīng)用。本文將使用該方法對(duì)參數(shù)和進(jìn)行擬合, 并與已有文獻(xiàn)結(jié)果進(jìn)行對(duì)比分析。

1 基于貝葉斯推斷的MCMC擬合方法

MCMC是一類基于貝葉斯推斷的、通過(guò)建立以目標(biāo)后驗(yàn)分布為平穩(wěn)分布的馬爾科夫鏈并對(duì)其進(jìn)行采樣的算法。有了后驗(yàn)分布, 模型參數(shù)的最可能值及其可靠取值范圍即可從中獲得。這類算法的思路是: 為了得到未知的模型參數(shù)在其取值空間中的概率分布(后驗(yàn)分布), 可以先猜測(cè)一個(gè)初始的分布(先驗(yàn)分布), 通過(guò)和觀測(cè)數(shù)據(jù)的比對(duì)來(lái)修正猜測(cè), 其用到了貝葉斯推斷的原理。修正的過(guò)程是不斷地從已知的分布(最初猜測(cè)的和后來(lái)調(diào)整的)中采樣得到當(dāng)前參數(shù)值, 再按照馬爾科夫鏈平穩(wěn)分布的要求產(chǎn)生新參數(shù)值的過(guò)程。理論上, 被更新的參數(shù)值的分布應(yīng)當(dāng)不斷逼近目標(biāo)后驗(yàn)分布, 當(dāng)馬爾可夫鏈?zhǔn)諗繒r(shí), 大量的參數(shù)值達(dá)到穩(wěn)定的統(tǒng)計(jì)分布, 便近似得到后驗(yàn)分布, 擬合過(guò)程結(jié)束。MCMC有多種對(duì)參數(shù)值進(jìn)行采樣和更新的策略, 本文采用了Metropolis-Hasting算法。貝葉斯推斷、馬爾可夫鏈和采樣策略是MCMC算法的核心。

1.1 貝葉斯推斷

貝葉斯推斷是一種利用觀測(cè)數(shù)據(jù)來(lái)修正先驗(yàn)假設(shè)的統(tǒng)計(jì)推斷。它的理論基礎(chǔ)是貝葉斯定理,

(|) =(|)()/()μ(|)(), (2)

式中:為觀測(cè)數(shù)據(jù);是對(duì)某1組未知參量的假設(shè)值, 本文中其為、、0、0四個(gè)自由參量;()為先驗(yàn)概率, 是在分析觀測(cè)數(shù)據(jù)前認(rèn)為參數(shù)值出現(xiàn)的概率;(|)為似然概率, 是模型參數(shù)值為時(shí), 恰好能得到觀測(cè)數(shù)據(jù)的可能性;(|)為后驗(yàn)概率, 表示在觀測(cè)數(shù)據(jù)時(shí), 模型參數(shù)值恰好為的概率;()為邊緣概率, 是所有可能的參數(shù)取值下觀測(cè)到數(shù)據(jù)的概率, 是一個(gè)常數(shù)。

后驗(yàn)概率可以理解為借助觀測(cè)數(shù)據(jù), 對(duì)先驗(yàn)概率進(jìn)行修正后的概率。本文所做擬合的目的是希望得到模型參數(shù)的分布, 在實(shí)際使用過(guò)程中, 要先假設(shè)參數(shù)服從某種形式的概率分布函數(shù)(下文稱為先驗(yàn)分布), 然后從中隨機(jī)抽取參數(shù)值, 代入模型, 再計(jì)算似然概率。即可得到后驗(yàn)概率, 其與真正的后驗(yàn)概率存在一個(gè)常系數(shù)的差別。

為了得到每個(gè)參數(shù)的后驗(yàn)分布, 以及其中的最大后驗(yàn)值和參數(shù)可靠區(qū)間, 一種做法是在參數(shù)定義域內(nèi)計(jì)算所有后驗(yàn)值。在模型復(fù)雜、參數(shù)維度高的場(chǎng)合, 這種做法將需要大量的計(jì)算資源, 因?yàn)槊總€(gè)參數(shù)的高概率分布區(qū)域可能落在相對(duì)比定義域小得多的范圍內(nèi)。一般的蒙特卡洛算法由于沒(méi)有一個(gè)高效的搜索策略, 將會(huì)把大量的計(jì)算資源浪費(fèi)在參數(shù)分布概率低的區(qū)域。為了解決這個(gè)問(wèn)題, 人們提出對(duì)后驗(yàn)分布進(jìn)行采樣, 通過(guò)大量樣本對(duì)后驗(yàn)分布進(jìn)行近似。本文選擇了MCMC算法, 它為多維分布的采樣提供了一個(gè)高效率的方案。

1.2 馬爾可夫鏈

對(duì)后驗(yàn)分布進(jìn)行采樣意味著不斷地更新模型的參數(shù)值, 這可以看成狀態(tài)的變化過(guò)程。一個(gè)隨機(jī)過(guò)程中, 如果下一個(gè)狀態(tài)的出現(xiàn)概率只取決于當(dāng)前狀態(tài), 而與之前任一狀態(tài)都無(wú)關(guān), 即P(X+1=|1=1,2=2,×××,X=x) =P(X+1=|X=x), 那么這個(gè)隨機(jī)過(guò)程就叫做馬爾可夫鏈。

從當(dāng)前狀態(tài)過(guò)渡到下一個(gè)狀態(tài)的概率可以由轉(zhuǎn)移矩陣描述,, 其中任一元素p代表狀態(tài)從變?yōu)榈陌l(fā)生概率, 它其實(shí)是一個(gè)條件概率, 即在給定狀態(tài)的前提下,狀態(tài)出現(xiàn)的概率。因此, 轉(zhuǎn)移概率矩陣描述了在狀態(tài)空間中, 所有可能的狀態(tài)變化的概率, 也包括一個(gè)狀態(tài)之后出現(xiàn)相同狀態(tài)的概率, 如11、22等。

馬爾可夫鏈有一個(gè)特性, 只要滿足每一個(gè)狀態(tài)的出現(xiàn)概率都不小于0, 以及狀態(tài)轉(zhuǎn)換不存在周期性這2個(gè)條件, 無(wú)論這條鏈開始于哪一個(gè)狀態(tài)點(diǎn)(0), 持續(xù)使用相同的轉(zhuǎn)移矩陣, 經(jīng)過(guò)一定次數(shù)()的狀態(tài)轉(zhuǎn)換后, 它的各種狀態(tài)的分布將會(huì)收斂到一個(gè)穩(wěn)定的分布——馬爾可夫鏈的平穩(wěn)分布(用表示), 且0×?。在滿足Kolmogorov準(zhǔn)則的情況下, 當(dāng)馬爾可夫鏈?zhǔn)諗康狡椒€(wěn)分布后, 狀態(tài)的轉(zhuǎn)移將滿足“細(xì)致平衡(Detailed balance)”條件[8], 即=。其中:π和π為處在和狀態(tài)的平穩(wěn)分布;為從狀態(tài)轉(zhuǎn)換到狀態(tài)的馬爾科夫鏈轉(zhuǎn)移矩陣。

1.3 Metropolis-Hasting算法

在模型參數(shù)的擬合中, 雖然可以通過(guò)式(2)計(jì)算后驗(yàn), 但無(wú)法直接對(duì)整個(gè)后驗(yàn)分布進(jìn)行采樣。這個(gè)問(wèn)題可以通過(guò)把參數(shù)的后驗(yàn)分布設(shè)為馬爾可夫鏈的平穩(wěn)分布來(lái)解決。這樣, 對(duì)后驗(yàn)的采樣就轉(zhuǎn)化成對(duì)收斂后的馬爾科夫鏈的采樣, 因此, 關(guān)鍵是如何使馬爾可夫鏈?zhǔn)諗康狡椒€(wěn)分布。Metropolis-Hasting算法正是利用了馬爾可夫鏈的細(xì)致平衡特性, 通過(guò)尋找合適的轉(zhuǎn)移矩陣來(lái)實(shí)現(xiàn)這一目的[9–10]。盡管一開始轉(zhuǎn)移矩陣是未知的, 但是可以先假設(shè)一個(gè), 再通過(guò)細(xì)致平衡進(jìn)行修改。一個(gè)可行的方案是將假設(shè)參數(shù)的先驗(yàn)分布(·|)作為初始的轉(zhuǎn)移矩陣。顯然, 開始這個(gè)轉(zhuǎn)移矩陣并不滿足細(xì)致平衡條件, 例如可能出現(xiàn)(+1|) >+1(|+1), 該式可以解讀為從狀態(tài)轉(zhuǎn)移到狀態(tài)+ 1的概率比反過(guò)來(lái)轉(zhuǎn)移的概率大。為了滿足細(xì)致平衡, 必須引入一個(gè)因子(+1|)來(lái)減少?gòu)臓顟B(tài)到狀態(tài)+ 1的轉(zhuǎn)換概率, 從而實(shí)現(xiàn)平衡, 即(+1|)(+1|) =+1(|+1)。因此可以得到(+1|) =+1(|+1)/(+1|), 其中(t+1|θ)稱為接受概率, 用于判斷是否將+ 1狀態(tài)接受為一個(gè)新的狀態(tài), 還是要繼續(xù)保留狀態(tài)。在實(shí)際操作中, 若(t+1|θ) > 1, 就接受+ 1狀態(tài)作為新的當(dāng)前態(tài)。若(θ+1|θ) < 1, 則從均勻分布(0, 1)中產(chǎn)生一個(gè)隨機(jī)值, 如果(θ+1|θ) >, 則接受+ 1狀態(tài), 否則保持當(dāng)前狀態(tài)。+ 1狀態(tài)被接受后, 參數(shù)的先驗(yàn)分布也隨之被修改。重復(fù)執(zhí)行上述過(guò)程, 經(jīng)過(guò)大量迭代, 產(chǎn)生參數(shù)的馬爾可夫鏈最終收斂到平穩(wěn)分布。此時(shí)對(duì)平穩(wěn)分布進(jìn)行采樣相當(dāng)于對(duì)后驗(yàn)分布進(jìn)行采樣。

假設(shè)被檢測(cè)的后驗(yàn)分布是中值為5、標(biāo)準(zhǔn)差為1.5的高斯分布, 采用的初始先驗(yàn)是中值為50、標(biāo)準(zhǔn)差為0.5的高斯分布。經(jīng)過(guò)2 000次迭代后, 先驗(yàn)被修正成中值為5.09, 標(biāo)準(zhǔn)差為1.54的高斯分布。圖2給出了MCMC的參數(shù)軌跡歷史, 從圖2可以看出在第1000次隨機(jī)過(guò)程之前, 先驗(yàn)一直在被修正。在1 000次之后, 參數(shù)值的變化已基本收斂, 馬爾科夫鏈?zhǔn)諗康狡椒€(wěn)分布, 即所求的后驗(yàn)分布, 那么收斂后的采樣點(diǎn)即可用于后續(xù)的統(tǒng)計(jì)分析。

圖2 馬爾可夫鏈的參數(shù)序列

1.4 擬合步驟

給定一顆脈沖星, 本文使用RVM對(duì)該星若干脈沖相位的偏振位置角數(shù)據(jù)ψ、(= 1, 2, 3,×××,)進(jìn)行擬合, 步驟如下。

(1) 為4個(gè)模型參數(shù)、、0和0指定初始的先驗(yàn)分布, 本文取高斯分布, 然后根據(jù)先驗(yàn)分布隨機(jī)生成一組參數(shù)值θ(,,0,0)作為當(dāng)前參數(shù)值。

(2) 把參數(shù)值θ代入似然概率表達(dá)式計(jì)算后驗(yàn)概率(θ|)。本文采用的似然概率表達(dá)式為, 其中mi是在相位上RVM模型預(yù)言的線偏振位置角, 可通過(guò)求解式(1)得到。

(3) 利用當(dāng)前先驗(yàn)分布(·|θ), 從當(dāng)前參數(shù)值θ生成新的參數(shù)值θ+1。

(4) 把參數(shù)值θ+1代入似然概率表達(dá)式, 計(jì)算后驗(yàn)概率(θ+1|)。

(6) 若> 1, 則接受+1為當(dāng)前參數(shù)值, 即成為新的, 更新先驗(yàn)分布; 若< 1, 則從均勻分布(0, 1)產(chǎn)生一個(gè)隨機(jī)值, 如果>, 則接受θ+1為當(dāng)前參數(shù)值, 否則拒絕θ+1并保持θ。

(7) 重復(fù)步驟(3)～(6), 直至參數(shù)值變化收斂。

(8) 繼續(xù)產(chǎn)生充分多的參數(shù)樣本, 得到參數(shù)的統(tǒng)計(jì)分布, 確定其最高后驗(yàn)概率值對(duì)應(yīng)的參數(shù)值, 以此作為參數(shù)的最可能取值, 并根據(jù)分布確定可信區(qū)間。

2 樣本數(shù)據(jù)及結(jié)果

2.1 數(shù)據(jù)

為了與已有文獻(xiàn)的方法進(jìn)行對(duì)比, 檢驗(yàn)MCMC方法的有效性, 本文選擇了與文獻(xiàn)[2]中相同的樣本進(jìn)行擬合, 共有10顆具有高質(zhì)量觀測(cè)數(shù)據(jù)的脈沖星。數(shù)據(jù)來(lái)自European Pulsar Network (EPN)脈沖星輪廓數(shù)據(jù)庫(kù), 與文獻(xiàn)[2]一樣, 所選的1418 MHz的數(shù)據(jù)均來(lái)源于Weisberg等[11]的觀測(cè)。Everett等[2]指出, 由于觀測(cè)噪聲的影響, 由Stokes參量直接計(jì)算線偏振強(qiáng)度和線偏振位置角會(huì)帶來(lái)偏差, 因此, 本文按照文獻(xiàn)[2]的方法對(duì)線偏振位置角進(jìn)行預(yù)處理, 修正偏差。

與文獻(xiàn)[2]略有不同的是, 本文只保留了脈沖窗口中的線偏振位置角數(shù)據(jù), 在窗口之外輻射極弱地方的線偏振位置角誤差很大, 本文沒(méi)有采用。與文獻(xiàn)[2]相同的是, 針對(duì)少數(shù)脈沖星在一些相位段上的線偏振位置角完全偏離RVM模型曲線的情況, 在擬合時(shí)刪除了這些數(shù)據(jù), 這些星包括PSR J0826 + 2637、PSR J0953 + 0755、PSR J1543 + 0929和PSR J1932 + 1059。

2.2 參數(shù)擬合值和可信范圍的確定

經(jīng)過(guò)MCMC方法擬合后, RVM模型的每個(gè)參數(shù)都會(huì)有一系列采樣值。去掉后驗(yàn)分布穩(wěn)定前的采樣, 對(duì)剩下的采樣值做統(tǒng)計(jì)直方圖并進(jìn)行平滑, 可以得到參數(shù)的后驗(yàn)分布。將后驗(yàn)分布中出現(xiàn)概率最大處所對(duì)應(yīng)的參數(shù)值作為最可能的參數(shù)取值, 在其左右某一范圍內(nèi)的參數(shù)區(qū)間作為某一可信度水平上的取值區(qū)間, 本文選用95%可信區(qū)間。以PSR J0528 + 2200為例, MCMC擬合得到的4個(gè)參數(shù)的后驗(yàn)分布如圖3所示。其中統(tǒng)計(jì)直方圖是大約使用了4 × 105個(gè)采樣值得到的, 最可能的參數(shù)擬合值用概率密度最高處的豎線表示。對(duì)參數(shù)值的大小進(jìn)行升序排序, 2.5%與97.5%兩個(gè)分位點(diǎn)之間的取值區(qū)域作為該參數(shù)的95%可信區(qū)間, 對(duì)應(yīng)于圖3中左右兩側(cè)豎線之間的區(qū)間。圖3中包絡(luò)線為經(jīng)過(guò)平滑得到的概率密度分布曲線。這2種方法給出結(jié)果的差別非常小。在PSR J0528 + 2200的擬合中, 從產(chǎn)生樣本到輸出所有被接受的樣本一共耗時(shí)8.5 s。其它樣本中, 最少耗時(shí)4.9 s, 最多耗時(shí)31.9 s, 平均耗時(shí)13.4 s。

圖3 PSR J0528 + 2200的RVM模型參數(shù)后驗(yàn)概率密度分布

2.3 擬合結(jié)果

本文使用MCMC算法成功地?cái)M合了10顆脈沖星, 得到了模型參數(shù)的最可能取值和95%可信區(qū)間, 擬合結(jié)果見(jiàn)表1, 用4個(gè)參數(shù)的最可能值得到的擬合曲線見(jiàn)圖4～7。為了方便對(duì)比, 文獻(xiàn)[2]的擬合結(jié)果也分別在表1和圖4～7中給出。由于文獻(xiàn)[2]沒(méi)有給出0的擬合值, 因此表1未列出該參數(shù)。圖4中: 圖4(a, b, c)為本文MCMC算法的擬合結(jié)果, 圖4(d, e, f)為L(zhǎng)M算法擬合結(jié)果[2]; 圖4(a, d)為脈沖偏振輪廓, 曲線1為總強(qiáng)度, 曲線2為線偏振強(qiáng)度, 曲線3為圓偏振強(qiáng)度; 圖4(b, e)為線偏振位置角數(shù)據(jù)和RVM模型的最佳擬合曲線; 圖4(c, f)為擬合殘差, 其計(jì)算方法為: 觀測(cè)值減去擬合值, 再除以觀測(cè)誤差。圖5、6、7中各分圖說(shuō)明與圖4相同。

表1 2種方法擬合結(jié)果的比較/(°)

續(xù)表1

J1841+0912B1839+0986.1±11.42.30±0.040.34±0.02 J1917+1353B1915+1373.0±19.45.4±0.56.8±0.1 J1918+1444B1916+14118.0±33.4-1.00±0.330.38±0.01 J1932+1059B1929+1035.97±0.9525.55±0.87-19.6±0.2

圖4 PSR J0304 + 1932的RVM模型擬合曲線

圖5 PSR J0528 + 2200的RVM模型擬合曲線

圖6 PSR J0826 + 2637的RVM模型擬合曲線

圖6(a, b, c)在相位(175°, 178.6°)范圍內(nèi)的線偏振位置角不滿足RVM, 未用于擬合, 與文獻(xiàn)[2]的處理相同。

圖7 PSR J0953 + 0755的RVM模型擬合曲線

圖7(a, b, c)在相位(175.8°, 181.0°)范圍內(nèi)的線偏振位置角不滿足RVM, 未用于擬合, 與文獻(xiàn)[2]的處理相同。

從表1和圖4～7可以看出, 本文方法得到的和與Everett等[2]的結(jié)果十分接近, 尤其是PSR J0304 + 1932、PSR J0528 + 2200、PSR J0826 + 2637、PSR J0953 + 0755、PSR J1917 + 1353和PSR J1918 + 1444這6顆星。LM方法僅給出了對(duì)稱的不確定度, 而MCMC方法則能夠給出非對(duì)稱可信區(qū)間。

3 結(jié)論

本文首次將MCMC算法應(yīng)用于脈沖星輻射幾何限定的研究, 通過(guò)使用旋轉(zhuǎn)矢量模型對(duì)脈沖星射電波段線偏振位置角的觀測(cè)數(shù)據(jù)進(jìn)行擬合, 得到了10顆脈沖星的磁傾角和碰撞角。與已有文獻(xiàn)基于LM算法的擬合結(jié)果相比, 本文的結(jié)果與相關(guān)文獻(xiàn)的結(jié)果高度一致, 證明了該方法的有效性。與LM算法和格點(diǎn)計(jì)算擬合方法相比, MCMC算法具有如下優(yōu)點(diǎn): (1) 對(duì)多維參數(shù)的擬合效率高; (2) 對(duì)參數(shù)的初始值不敏感, 能夠較好地抑制過(guò)擬合; (3) 能夠給出更為可靠的參數(shù)不確定度估計(jì); (4) 相比格點(diǎn)計(jì)算, 參數(shù)擬合精度高。在磁傾角和碰撞角的擬合工作中, 效率和過(guò)擬合是2個(gè)主要的問(wèn)題。本文研究顯示MCMC方法在這2方面, 特別是效率方面有突出的優(yōu)勢(shì), 為后續(xù)開展大樣本的擬合工作提供了有力的工具。

致謝:感謝中國(guó)科學(xué)院高能所夏俊卿及北京大學(xué)李柯伽給予本工作的討論。

[1] Radhakrishnan V, Cooke D J. Magnetic poles and the polarization structure of pulsar radiation [J]. ApL, 1969, 3: 225–229.

[2] Everett J E, Weisberg J, M. Emission beam geometry of selected pulsars derived from average pulse polarization data [J]. ApJ, 2001, 553: 341–357.

[3] Keith M J, Johnston S, Bailes M, et al. The high time resolution universe pulsar survey-IV. discovery and polarimetry of millisecond pulsars [J]. MNRAS, 2012, 419: 1 752–1 765.

[4] Rookyard S C, Weltevrede P, Johnston S. Constraints on viewing geometries from radio observations of γ-ray-loud pulsars using a novel method [J]. MNRAS, 2015, 446: 3 367–3 388.

[5] Marquardt D W. An algorithm for Least-Squares estimation of nonlinear parameters [J]. Journal of the Society for Industrial and Applied Mathematics, 1963, 11(2): 431–441.

[6] Wang H G, Pi F P, Zheng X P, et al. A Fan Beam Model for radio pulsars. I. observational evidence [J]. ApJ, 2014, 789: 73–102.

[7] Powell M. J. D. An efficient method for finding the minimum of a function of several variables without calculating derivatives [J]. Computer Journal, 1964, 7(2): 155–162.

[8] Boltzmann L. Lectures on gas theory [M]. Berkeley: U of California Press, 1964.

[9] Hastings W K. Monte Carlo sampling methods using Markov chains and their applications [J]. Biometrika, 1970, 57(1): 97–109.

[10] Chib S, Greenberg E. Understanding the Metropolis-Hastings Algorithm [J]. The American Statistician, 1995, 49(4): 327–335.

[11] Weisberg J M, Cordes J M, Lundgren S C, et al. Arecibo 1 418 MHz polarimetry of 98 pulsars: full stokes profiles and morphological classifications [J]. ApJS, 1999, 121: 171–217

(責(zé)任編校: 江河)

A bayesian approach to fitting the polarization position angle of pulsars

Chen Weiwei, Wang Hongguang, Zhang Yanrong

(School of Physics and Electronic Engineering, Guangzhou University, Guangzhou 510006, China)

For pulsar emission beams, a widely used relationship between the polarization position angle (PPA) and the pulse phase is given by the Rotation Vector Model (RVM). Current methods of fitting the PPA with RVM are on the basis of Levenberg-Marquardt nonlinear fitting or grid searching algorithms, however, they suffer from inefficiency and overfitting problems. A Markov Chain Monte Carlo (MCMC) algorithm on the basis of Bayesian inference is developed and applied to the fitting of PPA for 10 pulsars, which have been studied in literature with old methods. It is shown that the MCMC algorithm can not only obtain consistent results as those in literatures, but also have advantages in efficiency, the abilities of overcoming over-fitting and estimating the uncertainties of model parameters.

pulsar; emission beam; polarization position angle; Rotating Vector Model

http://www.cnki.net/kcms/detail/43.1420.N.20160526.1115.002.html

10.3969/j.issn.1672–6146.2016.04.007

P 145.6

1672–6146(2016)04–0027–08

陳威威, anforone@msn.com;王洪光, hgwang@gzhu.edu.cn。

2016–05–20

國(guó)家自然科學(xué)基金(NSFC 11178001; 11573008)。