正規分布について考えるため、以下の図で記号を定義する。
出典: PDG
信頼区間(confidence interval) を定義したとき、有意度(σ)と、信頼区間から外れる確率(両側の場合はα、片側の場合はα/2)は以下のように計算できる。
from scipy.stats import norm for sigma in [1,2,3,4,5,6]: print("δ: {:.1f}σ --> α : {:.2e}".format(sigma,(1-norm.cdf(sigma))*2)) print("δ: {:.1f}σ --> α/2: {:.2e}".format(sigma,1-norm.cdf(sigma))) print() for p_value in [0.3,0.2,0.1,0.01,0.001,1e-4,1e-5]: print("α : {:.1e} --> δ: {:.2f}σ".format(p_value,norm.ppf(1-p_value/2))) print("α/2: {:.1e} --> δ: {:.2f}σ".format(p_value,norm.ppf(1-p_value)))
cdfはCumulative Distribution Function、つまり累積分布関数の頭文字、ppfはPercent Point Functionの頭文字でcdf(累積分布関数)の逆関数である。
出力
δ: 1.0σ --> α : 3.17e-01 δ: 1.0σ --> α/2: 1.59e-01 δ: 2.0σ --> α : 4.55e-02 δ: 2.0σ --> α/2: 2.28e-02 δ: 3.0σ --> α : 2.70e-03 δ: 3.0σ --> α/2: 1.35e-03 δ: 4.0σ --> α : 6.33e-05 δ: 4.0σ --> α/2: 3.17e-05 δ: 5.0σ --> α : 5.73e-07 δ: 5.0σ --> α/2: 2.87e-07 δ: 6.0σ --> α : 1.97e-09 δ: 6.0σ --> α/2: 9.87e-10 α : 3.0e-01 --> δ: 1.04σ α/2: 3.0e-01 --> δ: 0.52σ α : 2.0e-01 --> δ: 1.28σ α/2: 2.0e-01 --> δ: 0.84σ α : 1.0e-01 --> δ: 1.64σ α/2: 1.0e-01 --> δ: 1.28σ α : 1.0e-02 --> δ: 2.58σ α/2: 1.0e-02 --> δ: 2.33σ α : 1.0e-03 --> δ: 3.29σ α/2: 1.0e-03 --> δ: 3.09σ α : 1.0e-04 --> δ: 3.89σ α/2: 1.0e-04 --> δ: 3.72σ α : 1.0e-05 --> δ: 4.42σ α/2: 1.0e-05 --> δ: 4.26σ
例えば、ある観測結果が背景事象で説明できる確率が1×10-5のとき、それが正規分布の片側のみに現れる背景事象の場合は、正規分布に焼き直すと観測結果は4.26σの有意度であるといえる。
次に、ポアソン分布を考える。
Wikipediaより
例えば、設定した条件内に入ってくる背景事象の期待値が1の場合、観測イベントが0イベントになる確率は0.368、1イベントの確率は0.368、2イベントの確率は0.184となる。これを計算するには、
from scipy.stats import poisson for expected in [1,4,10]: for k in [0,1,2,3,4,5,10,15,20]: print("λ={} k={} P={:.3f}".format(expected, k, poisson.pmf(k, expected)))
とすればよい。出力は、
λ=1 k=0 P=0.368 λ=1 k=1 P=0.368 λ=1 k=2 P=0.184 λ=1 k=3 P=0.061 λ=1 k=4 P=0.015 λ=1 k=5 P=0.003 λ=1 k=10 P=0.000 λ=1 k=15 P=0.000 λ=1 k=20 P=0.000 λ=4 k=0 P=0.018 λ=4 k=1 P=0.073 λ=4 k=2 P=0.147 λ=4 k=3 P=0.195 λ=4 k=4 P=0.195 λ=4 k=5 P=0.156 λ=4 k=10 P=0.005 λ=4 k=15 P=0.000 λ=4 k=20 P=0.000 λ=10 k=0 P=0.000 λ=10 k=1 P=0.000 λ=10 k=2 P=0.002 λ=10 k=3 P=0.008 λ=10 k=4 P=0.019 λ=10 k=5 P=0.038 λ=10 k=10 P=0.125 λ=10 k=15 P=0.035 λ=10 k=20 P=0.002
背景事象の期待値が1イベントで、実験で観測されたイベントが10イベントの場合、この10イベントが全て背景事象で説明できる確率はどう求めればよいか?
これは、期待値が1イベントで、観測が0~9イベントまでの累積分布を求めれば良い。
from scipy.stats import poisson bg = 1 events = 10 print("p-value: {:.2e}".format(1-poisson.cdf(events-1, bg))) print("p-value: {:.2e}".format(poisson.sf(events-1, bg)))
この結果はどちらも 1.11e-07となるだろう。このp値が小さすぎる場合は、poisson.sf関数の方が正確な結果になる。
このp値を正規分布における有意度に焼き直すには、背景事象がポアソン分布の片側だけに出てくることに注意すると、norm.ppf(1-poisson.sf(events-1, bg))より5.2σが導出できるだろう。
