粒子物理與核物理實驗中的數(shù)據(jù)分析 - 清華大學_12908

1、粒子物理與核物理實驗中的粒子物理與核物理實驗中的數(shù)據(jù)分析數(shù)據(jù)分析楊振偉楊振偉清華大學清華大學第第十十一講：解譜法（開拆法）一講：解譜法（開拆法）2本講要點本講要點n數(shù)學公式數(shù)學公式, ,反應函數(shù)反應函數(shù)( (矩陣矩陣) )n求反應矩陣的逆求反應矩陣的逆n修正因子修正因子n正規(guī)化的解譜法正規(guī)化的解譜法n估計量的方差與偏置估計量的方差與偏置n正規(guī)化參數(shù)的選擇正規(guī)化參數(shù)的選擇n舉例舉例a) Tikhonov 規(guī)則規(guī)則b) MaxEnt 規(guī)則規(guī)則3圖像還原問題圖像還原問題一個常見的問題一個常見的問題：由于實驗儀器的原因而出現(xiàn)由于實驗儀器的原因而出現(xiàn)圖像變形圖像變形，例如，例如如果如果已知探測器響應（可

2、通過探測器模擬已知探測器響應（可通過探測器模擬得到探測器響應的形式），得到探測器響應的形式），探測器響應探測器響應實驗觀測分布實驗觀測分布能否還原出不受實驗儀能否還原出不受實驗儀器影響的器影響的真實分布真實分布？Unfolding(解譜法解譜法)真實分布真實分布4解譜問題的表述解譜問題的表述(1)考慮隨機變量考慮隨機變量y，目標：尋找概率密度函數(shù)目標：尋找概率密度函數(shù)f(y)bin ( )1,.,jjjtotjpf y dyjMp “真實的直方真實的直方圖圖” 1)如果已知參數(shù)化形式如果已知參數(shù)化形式 ，);(yf最大似然法最大似然法);(yf2)不知概率密度函數(shù)的形式，可以構造直方圖不知概率

3、密度函數(shù)的形式，可以構造直方圖可以構造直方圖可以構造直方圖(M個區(qū)間個區(qū)間)，概率密度函數(shù)離散化，概率密度函數(shù)離散化第第j個區(qū)間的概率為個區(qū)間的概率為為為 j( (或或pj) )構造估計量構造估計量( (參數(shù)的數(shù)目參數(shù)的數(shù)目= =區(qū)間的數(shù)目區(qū)間的數(shù)目M)離散化的離散化的p.d.f.p.d.f.5解譜問題的表述解譜問題的表述(2)問題問題：y 的測量不可能沒有誤差的測量不可能沒有誤差1）第）第i個區(qū)間的個區(qū)間的y在第在第j個區(qū)間測量到個區(qū)間測量到 （測量分辨率）（測量分辨率）2）第）第i個區(qū)間的某些事例沒有測量到個區(qū)間的某些事例沒有測量到 （測量效率）（測量效率）3）某些區(qū)間的事例完全測量不到）

4、某些區(qū)間的事例完全測量不到 （接收度）（接收度）后果后果：f(y) 模糊化，峰展寬，模糊化，峰展寬， 甚至分布形狀完全變形。甚至分布形狀完全變形。測量測量后果嚴重時需要考慮解譜法還原真實分布后果嚴重時需要考慮解譜法還原真實分布6響應矩陣響應矩陣測量誤差的影響可用積分方程表示測量誤差的影響可用積分方程表示: :y真值真值；觀測值觀測值:x( )(|)( )meastruefxR x y fy dy 離散化離散化1,1,.,MiijjjRiN 觀測直方圖觀測直方圖( (期待值期待值) )真實直方圖真實直方圖響應矩陣響應矩陣 Rij=P(觀測值在第觀測值在第 i 區(qū)區(qū)|真實值在第真實值在第 j 區(qū)區(qū)

5、)( | )R x y：響應函數(shù)，表示真值為：響應函數(shù)，表示真值為y觀測值為觀測值為x的概率的概率對對y求積分即可得到測量值為求積分即可得到測量值為x的概率的概率7響應矩陣響應矩陣數(shù)據(jù)：數(shù)據(jù)：1(,.,)NiinnnE n 這里這里Rij=P(觀測值在第觀測值在第 i 區(qū)區(qū)|真實值在第真實值在第 j 區(qū)區(qū))注意：注意： 是常數(shù)，是常數(shù)， 會受會受到統(tǒng)計漲落的影響。到統(tǒng)計漲落的影響。,n真實直方圖真實直方圖離散化的離散化的p.d.f.p.d.f. 觀測數(shù)據(jù)觀測數(shù)據(jù)和數(shù)據(jù)的期待值和數(shù)據(jù)的期待值8效率、本底效率、本底 有時，事例可能會不被探測到有時，事例可能會不被探測到: : 效率效率 11( |

6、)(|()NNijiijRPijPj 觀觀測測值值在在第第區(qū)區(qū) 真真實實值值在在第第區(qū)區(qū)觀觀測測值值在在全全范范圍圍 真真實實值值在在第第區(qū)區(qū)) )效效率率 i 是在觀測直方圖上預期的本底數(shù)目。是在觀測直方圖上預期的本底數(shù)目。 有時，無真實事例發(fā)生，也有事例被觀測到：有時，無真實事例發(fā)生，也有事例被觀測到：本底本底 1MiijjijR 真實直方圖真實直方圖第第 j 區(qū)探測效率區(qū)探測效率真值真值效率效率測量值測量值9各關鍵量匯總各關鍵量匯總“真實真實”直方圖直方圖: :概率概率: :觀測直方圖的期待值觀測直方圖的期待值: : 觀測直方圖觀測直方圖: :響應矩陣響應矩陣: :效率效率: :11(,

7、.,),MMtotjj 1(,.,)/Mtotppp 1(,.,)N 1(,.,)Nnnn ( | )ijRPij 觀觀測測值值在在第第區(qū)區(qū) 真真實實值值在在第第區(qū)區(qū)1NjijiR 預期的本底預期的本底: :1(,.,)N E nR 或關聯(lián)矩陣或關聯(lián)矩陣(各區(qū)間不獨立各區(qū)間不獨立)Vij=covni,nj一般通過一般通過構造構造 logL 或或 2尋求尋求 的估計量，這需要相關的的估計量，這需要相關的概率理論，例如概率理論，例如: :泊松分布泊松分布(各區(qū)間獨立各區(qū)間獨立)(;)!iiniiiiP nen M 個區(qū)間個區(qū)間N 個區(qū)間個區(qū)間10為什么要用解譜法？為什么要用解譜法？ 但是，不將實驗

8、數(shù)據(jù)進行解譜處理，結果發(fā)表后，有但是，不將實驗數(shù)據(jù)進行解譜處理，結果發(fā)表后，有關反應矩陣的知識將不在保留。而且，解譜后的分布可以直關反應矩陣的知識將不在保留。而且，解譜后的分布可以直接與各種理論的預言相比較，也可以與別的實驗經(jīng)過解譜以接與各種理論的預言相比較，也可以與別的實驗經(jīng)過解譜以后的分布相比較。后的分布相比較。 一般而言，我們并不需要解譜法，例如當比較現(xiàn)有理論一般而言，我們并不需要解譜法，例如當比較現(xiàn)有理論的預期值時，最好是將探測器相應疊加到理論中去。即在預的預期值時，最好是將探測器相應疊加到理論中去。即在預期值中包含探測器效應并與未修正的原始數(shù)據(jù)期值中包含探測器效應并與未修正的原始數(shù)據(jù)

9、 相比較。相比較。n 通常解譜后的結果更為有用，否則當反應矩陣不可恢復通常解譜后的結果更為有用，否則當反應矩陣不可恢復時，即使對結果又有新的理論解釋，也很難進行理論檢驗。時，即使對結果又有新的理論解釋，也很難進行理論檢驗。在粒子物理研究中，解譜法常用的領域為：在粒子物理研究中，解譜法常用的領域為： 強子強子結構函數(shù)結構函數(shù) 的譜函數(shù)的譜函數(shù)( (也就是強子不變質量譜也就是強子不變質量譜) ) 強子事例形狀分布強子事例形狀分布 粒子多重數(shù)分布粒子多重數(shù)分布11為什么用解譜法(舉例)中微子與鐵原子核相互作用，產(chǎn)生電子，測量電子的能譜。實驗上只可能測量電子在鐵中的射程。能量與平均射程的關系因子修正法

10、：能譜(圓圈)與真值(直方圖)解譜法12響應矩陣的逆響應矩陣的逆若數(shù)據(jù)是泊松分布若數(shù)據(jù)是泊松分布最大似然法的估計量為最大似然法的估計量為假設假設 的逆存在的逆存在: :R 1()R (;)!iiniiiiP nen 1log ( )(log)NiiiiLn 1()nRn 則有則有若若的非對角元太大，的非對角元太大，即區(qū)間寬度比分辨率要即區(qū)間寬度比分辨率要小時，會導致上式有很小時，會導致上式有很大的方差，以及在相鄰大的方差，以及在相鄰區(qū)間產(chǎn)生很強的負關聯(lián)。區(qū)間產(chǎn)生很強的負關聯(lián)。n？13錯誤的原因錯誤的原因(I)考慮一個簡單的例子分辨率很差完美探測器，其中:0:110,111121,2121Abb

11、bxxxbxA意義。的，但結果并沒有物理是嚴格的，是無偏有效的方法，盡管理論上接求反應矩陣獲得真值所以，通常情況下，直湮沒。信息完全被非物理振蕩放大后，結果中有用的這個隨機數(shù)被實際上是一個隨機數(shù)。時，考慮到統(tǒng)計漲落，尤其當xbbbb/12121項起決定作用很小時，第2112112111121111121111121212111bbbbbAxA解決辦法是進行平滑處理，消除無意義的統(tǒng)計漲落。但平滑會帶來偏向性，需要在漲落與偏向性之間找到平衡。14錯誤的原因錯誤的原因(II)考慮周期分布函數(shù)20),sincos(2)(10 xxbxaaxf2220122221()( )exp( )()22cossi

12、n(cossin)2exp,exp22yxg yf x dxg(y)xxg(y)xxab 如果每一點的測量都按照偏差為 的高斯分布彌散，那么測量到的分布為卷積將展開為和的函數(shù)：其中將完全占據(jù)統(tǒng)治地位。大后中的統(tǒng)計誤差被指數(shù)放和時尤其是對于一般函數(shù)，放大。這個統(tǒng)計誤差被都有一定的統(tǒng)計誤差，和但，和的系數(shù)布，則可嚴格求出真實分和的系數(shù)如果可以準確知道, 0,)2/(exp)(22babaxfg(y)實際上，非物理的劇烈漲落主要是有高頻部分造成的。平滑處理，主要是消除或壓低高頻部分無意義的漲落。15錯誤的原因錯誤的原因(III)假設假設 真的有精細結構真的有精細結構R作用后作用后( (理論上理論上)

13、 )得到得到R大部分精細結構被抹平，大部分精細結構被抹平，但會隱藏精細結構的信息。但會隱藏精細結構的信息。 但我們沒有但我們沒有 只有只有 （存在漲落）（存在漲落） n R-1“認為認為”這是隱藏的精細結構信息，不遺余力地恢復這是隱藏的精細結構信息，不遺余力地恢復這種精細結構，從而造成劇烈的振蕩效應這種精細結構，從而造成劇烈的振蕩效應。R-1 作用回作用回 應完全恢復精細結構：應完全恢復精細結構： 1R16重新研究最大似然法的解重新研究最大似然法的解(1)無偏無偏! !計算估計量的方差計算估計量的方差1( )ERE n 11,1111cov,() () cov,() ()Nijijikjlkl

14、k lNikjkkkURRnnRR ,iklklknnn coco假假設設 是是獨獨立立的的泊泊松松v v變變量量時時， 估計量的平均值估計量的平均值17重新研究最大似然法的解重新研究最大似然法的解(2)利用利用 RCF 邊界做無偏估計量邊界做無偏估計量211log()NikilklikliR RLUE 倒數(shù)后得到倒數(shù)后得到111() ()NijikjkkiURR 為了減小方差，必須引入一些偏置量為了減小方差，必須引入一些偏置量策略：接受小的偏置量策略：接受小的偏置量( (系統(tǒng)誤差系統(tǒng)誤差) )以換取大幅以換取大幅 減小方差減小方差( (統(tǒng)計誤差統(tǒng)計誤差) )。1()R 1log ( )(lo

15、g)NiiiiLn 與剛才得到的結果一致。與剛才得到的結果一致。ML估計量在所有無偏估計中給出的方差最小。估計量在所有無偏估計中給出的方差最小。但這個方差太大了！但這個方差太大了！18簡單方法：修正因子法簡單方法：修正因子法對對 做相同的分區(qū)，并取做相同的分區(qū)，并取, ()iiiiC n ()MCiiMCiC 修修正正因因子子 與與 用蒙特卡羅模擬得到（無本底）。用蒙特卡羅模擬得到（無本底）。MCi MCi 2cov,cov,ijijiijUCn n 通常通常 ，因此方差不會被放大。，因此方差不會被放大。(1)iCO ( (相當于相當于R-1-1取為對角矩陣取為對角矩陣) ) ( (R-1)

16、)ii= =Ci19修正因子法中的偏置問題修正因子法中的偏置問題iiibE,MCsigsigiiiiiiiiiMCsigiib = =其其中中。注意：注意：該偏置量存在把該偏置量存在把 拉向拉向 的傾向，造的傾向，造 成模型檢驗的困難。成模型檢驗的困難。 MC 1 1）如果分區(qū)寬度如果分區(qū)寬度 幾倍的分辨率，結果不會太壞。幾倍的分辨率，結果不會太壞。 2 2）實際應用中，中。該方法常用于事例形狀變量的分布研究）實際應用中，中。該方法常用于事例形狀變量的分布研究除非模擬采用的模型無誤，使得除非模擬采用的模型無誤，使得 ，否則，否則上式不為零，需要考慮對應的系統(tǒng)誤差。上式不為零，需要考慮對應的系統(tǒng)

17、誤差。修正因子法的偏置修正因子法的偏置MCii 20例子例子：脈沖形狀的還原：脈沖形狀的還原將理論將理論( (真實真實) )的直方圖除以受實驗儀器影響的直方圖除以受實驗儀器影響 的直方圖得到修正因子的直方圖得到修正因子將觀測直方圖乘以修正因子直方圖得到理論將觀測直方圖乘以修正因子直方圖得到理論 ( (真實真實) )的直方圖的直方圖21正規(guī)化的解譜法正規(guī)化的解譜法考慮考慮“合理的合理的”估計量，對選定的估計量，對選定的 滿足滿足Llogmaxlog ( )loglogLLL ( )log ( )( )LS( )S 正則化函數(shù)正則化函數(shù)( (光滑性的量度光滑性的量度) ) 正則化參量正則化參量(

18、(其選擇與給定的其選擇與給定的 對應對應) )Lloglog Ln 描描述述了了數(shù)數(shù)據(jù)據(jù)與與期期待待值值之之間間的的 “距距離離”。 估計量滿足該不等式并且最光滑，等價于估計量滿足該不等式并且最光滑，等價于將下式求最大值將下式求最大值 22正規(guī)化的解譜法（續(xù)一）正規(guī)化的解譜法（續(xù)一）另外，要求解譜后對總事例數(shù)的估計為無偏的另外，要求解譜后對總事例數(shù)的估計為無偏的tot1,NNiijjii jRn 1( , )log ( )( )NtotiiLSn : :拉格朗日乘子拉格朗日乘子在約束情況下將下式求最大值在約束情況下將下式求最大值, R 因因所所以以是是的的函函數(shù)數(shù)1/0Nitotin 顯然，顯

19、然，23正規(guī)化的解譜法正規(guī)化的解譜法（續(xù)二）（續(xù)二）0 給出最光滑的解給出最光滑的解( (與數(shù)據(jù)無關與數(shù)據(jù)無關) )給出最大似然解給出最大似然解( (方差太大方差太大) )( )S 需要需要正規(guī)化函數(shù)正規(guī)化函數(shù) 以及選取以及選取 值的方案值的方案。所得到的估計量的好壞由它們的偏置和方差所得到的估計量的好壞由它們的偏置和方差來判斷。來判斷。a) Tikhonov 規(guī)則規(guī)則b) MaxEnt 規(guī)則規(guī)則1( , )log ( )( )NtotiiLSn 24Tikhonov 規(guī)則規(guī)則取光滑度等于第取光滑度等于第 k 階導數(shù)平方的均值階導數(shù)平方的均值，有，有2( )( ),1,2.ktruetruek

20、d fyS fydykdy 其其中中通常取通常取 k=2，使得使得 S 約等于曲率平方的平均約等于曲率平方的平均值。對直方圖而言，也就是值。對直方圖而言，也就是22121( )(2)MiiiiS 注意注意：2 階階導數(shù)對導數(shù)對直方圖的直方圖的第一和最后的區(qū)第一和最后的區(qū) 間沒有很好的定義間沒有很好的定義。Sov. Math.5(1963)103525Tikhonov 規(guī)則（續(xù)）規(guī)則（續(xù)）如果在如果在 下，下，采用采用Tikhonov(k=2)規(guī)則規(guī)則21log2L 2( , )log ( )( )( )( )2LSS 是是 i 的二次項。的二次項。對對 求偏微分求偏微分，給出線性方，給出線性方

21、程，得到程，得到 估計值與方差。估計值與方差。在高能物理界現(xiàn)有好幾個現(xiàn)成的程序在高能物理界現(xiàn)有好幾個現(xiàn)成的程序：RUN，Blobel, SVD, Hcker,26最大熵最大熵(MaxEnt)規(guī)則規(guī)則另一種表征光滑度的方法基于熵。對于一組概另一種表征光滑度的方法基于熵。對于一組概率而言，熵定義為率而言，熵定義為1logMiiiHpp 所有所有 相等意味著熵最大相等意味著熵最大( (最光滑最光滑) )ip有一個有一個 ，其它為零，則意味著熵最小。，其它為零，則意味著熵最小。1ip Ann. Rev. Astron. Astrophys.24 (1986)127totiip27最大熵最大熵(MaxE

22、nt)規(guī)則（續(xù)）規(guī)則（續(xù)）用熵作為正規(guī)化函數(shù)，用熵作為正規(guī)化函數(shù)，1()()logl(ogMMiiitottottotSH 填填入入個個區(qū)區(qū)間間中中各各種種可可能能的的總總數(shù)數(shù)) )有時侯，根據(jù)貝葉斯統(tǒng)計有時侯，根據(jù)貝葉斯統(tǒng)計這里，我們仍然采用經(jīng)典近似這里，我們仍然采用經(jīng)典近似: : 估計量的好壞估計量的好壞由偏置，方差來判斷。由偏置，方差來判斷。注意注意：熵并不取決于區(qū)間的順序。：熵并不取決于區(qū)間的順序。( )S的先驗概率密度函數(shù)的先驗概率密度函數(shù)(?)(?)28 的方差與偏置的方差與偏置一般來說，決定一般來說，決定 的方程是非線性的。當?shù)玫姆匠淌欠蔷€性的。當?shù)玫秸?guī)函數(shù)后，將到正規(guī)函數(shù)后，

23、將 在在 附近展開附近展開)(nobsn)(n1222222( )(), ,1,.,1 1,.,1,0, 1,1, , 1,.,1,.,1, 1,1,.,.obsobsijijiijijjnA B nni jMAiM jMiMjMiM jNnBiMjNn ，G. Cowan, StatisticalData Analysis, OxfordUniversity Press(1998) 為為非非正正規(guī)規(guī)的的似似然然函函數(shù)數(shù)29 的方差與偏置的方差與偏置( (續(xù)續(xù)) )利用誤差傳遞得到協(xié)方差利用誤差傳遞得到協(xié)方差cov,ijijU 1,TUCVCCA B 其其中中以及對偏置的估計量，以及對偏置的估

24、計量，,iiibE11()(),NNiiijjjjjjjjbCnnn 此處此處 而且通常情況下而且通常情況下. Rncov,ijijVn n 30正規(guī)化參數(shù)正規(guī)化參數(shù) 的選取的選取 決定了決定了置于數(shù)據(jù)的權重大小以便能與光滑度置于數(shù)據(jù)的權重大小以便能與光滑度相比較，相比較， = 0 給出最大的光滑估計值，并與數(shù)據(jù)無給出最大的光滑估計值，并與數(shù)據(jù)無關。因此雖然方差為零，但有明顯的偏置。而取大關。因此雖然方差為零，但有明顯的偏置。而取大的的 ，則回到高度振蕩無偏的最大似然解。，則回到高度振蕩無偏的最大似然解。 為了在偏置與方差之間達到最大平衡為了在偏置與方差之間達到最大平衡: :選擇選擇 使均值誤

25、差的平方（使均值誤差的平方（MSE）最小最小 2211, 11MMiiiiiiiiiUbMSEUbMSEMM 或或帶帶權權重重的的或要求偏置不大于它自身的估計方差或要求偏置不大于它自身的估計方差 ?？梢哉铱梢哉业降?的值使得的值使得221cov ,.MibijijiiibMWb bW 這這里里iiW31正規(guī)化參數(shù)選擇的重要性與數(shù)據(jù)無關，光滑，完全MC:0無物理意義劇烈振蕩，完全數(shù)據(jù)，:32迭代法解譜(Bayes方法)NjjkkjkijiiNjjkkiiNjjiiijnpRpRnpkjPpijPnjiPjiPR1111)|()|(1)|(1Bayes)|(區(qū)間真值在區(qū)間測量值在區(qū)間真值在區(qū)間測量

26、值在區(qū)間測量值在區(qū)間真值在定理，可以得到估計量利用區(qū)間真值在區(qū)間測量值在果趨于劇烈振蕩。當?shù)螖?shù)很大時，結。數(shù)，之后用于實驗數(shù)據(jù)用測試數(shù)據(jù)確定迭代次迭代。，之后用選取初始totiii/p,N),(ip1 應用貝葉斯定理33RooUnfold: SVD方法解譜MC:( 22bin)2 2()Rx MCbxbhxbA可以通過獲得：真值模擬測量結果作散點分布設 和 都分 個實際上獲得了一個的數(shù)組 或矩陣1( , ) ,( , )TijiRAj iRAAj i即1111122121122222(1,1) :,(1,2) :,(2,1) :,(2,2):,AxbRAxbRAxbRAxbR真值為測量值為

27、 的事例數(shù)真值為測量值為 的事例數(shù)真值為測量值為 的事例數(shù)真值為測量值為 的事例數(shù)解譜需要先得到響應矩陣解譜需要先得到響應矩陣R，并進行求逆。，并進行求逆。對于性質不太好的矩陣，對于性質不太好的矩陣，SVDSVD方法更可靠一些。方法更可靠一些。34Singular Value Decomposition(SVD)TjjiiTTTTTTUVSAjiSSIVVI,UUVUmnAAeigdiagSmmAAVnnAAUUSVAAmn11,)(。顯然如果滿足對角矩陣矩陣的本征矢構成的矩陣的本征矢構成的其中可以分解成矩陣任意的列。和，基分別是矩陣和其系數(shù)分別構成向量的線性組合，都分解成正交歸一向量和的結果

28、是把，定義VUzdxbbAbUVSdVSVzxdSzdSzbUdxVzbUxSVbxUSVbAxTTTTTTSVD1111斷。的元素分布可以找到截需要截斷，從物理的，高頻部分等價地說，低頻部分是準正態(tài)分布。中相應的元素的服從標向量頻部分，相應的，對應于傅立葉分解中高的情況維例子中參見振蕩的根源部分是結果劇烈實際上，奇異值很小的的奇異值。矩陣對角元素為聯(lián)系起來它們通過對角矩陣，變?yōu)?，未知量變?yōu)橹?，已知量ddASzxdb),02(,SVD35RooUnfold的下載、編譯和使用nRooUnfold中提供了SVD以及Bayes的解譜法。下載和編譯使用方法見網(wǎng)頁：http:/hepunx.rl.ac.uk/adye/software/unfold/RooUnfold.html在training服務器上，已經(jīng)編譯好放到ROOT的庫文件目錄中，不需要另行下載編譯。使用時只要在ROOT腳本文件中加上 gSystem-Load(“l(fā)ibRooUnfold”);36RooUnfoldSvd正規(guī)化參數(shù)的選擇nRooUnfoldSvd中正規(guī)化函數(shù)的選擇實際上是Tikhonov(k=2)規(guī)則。n使用RooUnfoldSvd后自動生成UnfoHist.root文件，其中包含d的直方圖，從d的直方圖可以讀出正規(guī)化參數(shù)，即kterm=k,其中k之后各個bin的值滿足均值為0方差為1的高斯