はじめに

X線天文学では、天体から届いた光子を数えることで、その天体の性質を調べます。

たとえば、ある領域で何個の X線イベントが検出されたかを数え、そのカウント数から天体の明るさやスペクトルを調べます。

ここで注意したいのは、source region で観測されたカウントには、天体由来の信号だけでなく背景成分も含まれているということです。

ざっくり書けば、source region で観測される total counts は、

C
=
A
+
B_{\rm src}

のような量です。

ここで、

$C$：source region で観測された total counts
$A$：天体由来の counts
$B_{\rm src}$：source region に含まれる background counts

です。

私たちが知りたいのは、天体由来の信号 $A$ です。

そのため、background region などから source region 内の background を推定し、それを $\hat{B}$ と書くと、素朴には

N_{\rm net}
=
C
-
\hat{B}

として、背景差し引き後の net counts を求めたくなります。

この操作自体はとても自然です。

しかし、カウントデータを統計的に扱うときには、この「引く」という操作には少し注意が必要です。

なぜなら、ポアソン分布に従うカウント同士を引くと、その結果はもう単純なポアソン分布ではなくなるからです。

もちろん、カウントが十分大きい場合には Gaussian 近似で扱えることもあります。

一方で、低カウントでは分布の離散性や非対称性が効いてきます。

そこでこの記事では、

ポアソン分布に従うカウントを引くと、どんな分布になるのか？

を Python で可視化してみます。

最終的には、背景差し引き後の値だけを見ると、なぜ誤差を直感しにくくなるのかを考えます。

ポアソン分布は足すと壊れない

まず、ポアソン分布の足し算を見ます。

X \sim \mathrm{Poisson}(\lambda_1)

Y \sim \mathrm{Poisson}(\lambda_2)

で、$X$ と $Y$ が独立なら、

X+Y \sim \mathrm{Poisson}(\lambda_1+\lambda_2)

になります。

つまり、ポアソン分布に従う独立なカウントを足すと、またポアソン分布になります。

これはかなり扱いやすい性質です。

平均カウントが $\lambda_1$ の成分と、平均カウントが $\lambda_2$ の成分を足すと、平均カウントは単純に

\lambda_1+\lambda_2

になります。

X線観測で source region を見ると、そこには天体由来の counts と background counts が一緒に入っています。

つまり、

C
=
A
+
B_{\rm src}

のような形です。

天体由来の光子も background の光子も、それぞれ独立なポアソン過程として考えられるなら、観測される total counts $C$ もまたポアソン分布として扱えます。

これは、観測されたカウントをそのままポアソン統計の枠組みで扱いやすい理由の一つです。

では、引くとどうなるか

では、引き算はどうでしょうか。

ここでは、source region で観測された total counts を $C$、background から推定した counts を $\hat{B}$ とします。

背景差し引き後の net counts は、

N_{\rm net}
=
C
-
\hat{B}

です。

ここで、簡単のために

C \sim \mathrm{Poisson}(\lambda_C)

\hat{B} \sim \mathrm{Poisson}(\lambda_B)

のように、どちらも独立なポアソン分布に従うと考えます。

このとき、

N_{\rm net}
=
C
-
\hat{B}

はポアソン分布にはなりません。

理由は単純です。

ポアソン分布は、

0,1,2,3,\cdots

のような非負の整数しか取りません。

しかし、差

C-\hat{B}

は負の値を取ることがあります。

たとえば、

C=3,\quad \hat{B}=7

なら、

C-\hat{B}=-4

です。

これは「天体からの信号が負」という意味ではありません。

あくまで、揺らぎを持つ2つのカウントの差が負になることがある、ということです。

つまり、背景を引いた後の net counts は、元の観測カウントと同じ意味でのポアソンデータではありません。

ここが、X線天文学で背景を単純に引いた量だけを見るときに、少し慎重になる理由の一つです。

誤差伝搬として見ると

まずは、Gaussian 近似の気持ちで考えてみます。

独立な2つの量 $C$, $\hat{B}$ に対して、

N_{\rm net}
=
C
-
\hat{B}

とすると、誤差伝搬では

\sigma_{\rm net}^2
=
\sigma_C^2
+
\sigma_{\hat{B}}^2

になります。

引き算なのに、分散は足し算です。

これは、符号がマイナスでも、不確かさとしては両方のばらつきが効くためです。

ポアソン分布では、分散は平均に等しいので、

\sigma_C^2
=
\lambda_C

\sigma_{\hat{B}}^2
=
\lambda_B

です。

したがって、

\sigma_{\rm net}^2
=
\lambda_C
+
\lambda_B

となります。

一方で、平均は

E[N_{\rm net}]
=
\lambda_C
-
\lambda_B

です。

つまり、

E[C-\hat{B}]
=
\lambda_C-\lambda_B

ですが、

\mathrm{Var}(C-\hat{B})
=
\lambda_C+\lambda_B

になります。

ここで大事なのは、

平均は引き算されるが、分散は足し算される

ということです。

背景を引くと、平均値としては小さくなります。

しかし、background の揺らぎまで消えるわけではありません。

ただし、分布そのものは Gaussian とは限らない

ここまでの話は、誤差伝搬として見るとかなり自然です。

カウント数が十分大きければ、

C-\hat{B}
\approx
\mathcal{N}(\lambda_C-\lambda_B,\lambda_C+\lambda_B)

のように Gaussian 近似で考えることもできます。

しかし、ポアソン分布は本来、離散分布です。

さらに、低カウントでは Gaussian のような左右対称な滑らかな分布とは違い、離散性や非対称性が効いてきます。

そのため、低カウントの差を考えるときには、

誤差伝搬では分散はわかるけれど、分布の形そのものはどうなっているのか？

が気になります。

この「2つの独立なポアソン分布の差」として現れる分布が、Skellam分布です。

Skellam 分布

2つの独立なポアソン分布

C \sim \mathrm{Poisson}(\lambda_C)

\hat{B} \sim \mathrm{Poisson}(\lambda_B)

に対して、

N_{\rm net}
=
C-\hat{B}

を考えます。

このとき、

N_{\rm net}
\sim
\mathrm{Skellam}(\lambda_C,\lambda_B)

となります。

Skellam分布は、整数全体

\cdots,-2,-1,0,1,2,\cdots

に確率を持つ離散分布です。

ポアソン分布が非負の整数だけを取るのに対して、Skellam 分布は負の値も自然に取れます。

そのため、背景を引いた後の net counts の分布を考えるときに出てくる自然な分布です。

Python で可視化してみる

では、実際に Skellam 分布を Gaussian 近似と比較してみます。

ここでは、

$\lambda_C$：source region で観測される total counts の平均
$\lambda_B$：background として推定される counts の平均

とします。

背景差し引き後の net counts は、

N_{\rm net}
=
C
-
\hat{B}

です。

実装コードは Google Colab からもお試しいただけます：
https://colab.research.google.com/drive/1NrIWAlsE9GBb_ez_tluRgDb19v4LRTdC?usp=sharing

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import skellam, norm

# ======================================================
# Parameters
# ======================================================
lambda_c = 20   # mean observed total counts in the source region
lambda_b = 12   # mean estimated background counts

# Difference k = C - B
k_min = -20
k_max = 40
k = np.arange(k_min, k_max + 1)

# Skellam distribution
p_skellam = skellam.pmf(k, mu1=lambda_c, mu2=lambda_b)

# Gaussian approximation
mean = lambda_c - lambda_b
sigma = np.sqrt(lambda_c + lambda_b)
p_gauss = norm.pdf(k, loc=mean, scale=sigma)

# Monte Carlo sample
rng = np.random.default_rng(0)
n_sample = 200_000

c = rng.poisson(lambda_c, n_sample)
b = rng.poisson(lambda_b, n_sample)
n_net = c - b

# ======================================================
# Plot
# ======================================================
fig, ax = plt.subplots(figsize=(8, 5))

bins = np.arange(k_min - 0.5, k_max + 1.5, 1)

ax.hist(
    n_net,
    bins=bins,
    density=True,
    alpha=0.35,
    label=r"Monte Carlo: $C-\hat{B}$"
)

ax.plot(
    k,
    p_skellam,
    marker="o",
    linestyle="-",
    label="Skellam distribution"
)

ax.plot(
    k,
    p_gauss,
    linestyle="--",
    label="Gaussian approximation"
)

ax.axvline(
    mean,
    linestyle=":",
    label=r"mean = $\lambda_C - \lambda_B$"
)

ax.set_xlabel(r"$k = C-\hat{B}$")
ax.set_ylabel("Probability")
ax.set_title(
    rf"Background-subtracted counts: "
    rf"$\lambda_C={lambda_c}$, $\lambda_B={lambda_b}$"
)
ax.legend()
ax.grid(alpha=0.3)

plt.show()

ここでは、

\lambda_C=20,\quad \lambda_B=12

としました。

このとき、背景差し引き後の平均は

\lambda_C-\lambda_B=8

です。

一方で、標準偏差は

\sqrt{\lambda_C+\lambda_B}
=
\sqrt{32}

です。

背景を引いたので、分布の中心は $8$ にあります。

しかし、分布の広がりは $\lambda_C+\lambda_B$ で決まります。

つまり、背景を引いた後の値は $8$ 付近に見えますが、その揺らぎには background のカウントも効いています。

いくつかの場合を比較する

次に、いくつかの $\lambda_C,\lambda_B$ について見てみます。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import skellam, norm

examples = [
    (5, 5),
    (20, 12),
    (50, 50),
    (5, 20),
]

fig, axes = plt.subplots(2, 2, figsize=(10, 8))

rng = np.random.default_rng(0)

for ax, (lam_c, lam_b) in zip(axes.ravel(), examples):
    mean = lam_c - lam_b
    sigma = np.sqrt(lam_c + lam_b)

    k_min = int(mean - 6 * sigma) - 3
    k_max = int(mean + 6 * sigma) + 3
    k = np.arange(k_min, k_max + 1)

    p_skellam = skellam.pmf(k, mu1=lam_c, mu2=lam_b)
    p_gauss = norm.pdf(k, loc=mean, scale=sigma)

    c = rng.poisson(lam_c, 100_000)
    b = rng.poisson(lam_b, 100_000)
    n_net = c - b

    bins = np.arange(k_min - 0.5, k_max + 1.5, 1)

    ax.hist(
        n_net,
        bins=bins,
        density=True,
        alpha=0.3,
        label="MC"
    )

    ax.plot(k, p_skellam, "o-", markersize=3, label="Skellam")
    ax.plot(k, p_gauss, "--", label="Gaussian")
    ax.axvline(mean, linestyle=":")

    ax.set_title(
        rf"$\lambda_C={lam_c},\ \lambda_B={lam_b}$"
    )
    ax.set_xlabel(r"$C-\hat{B}$")
    ax.set_ylabel("Probability")
    ax.grid(alpha=0.3)

axes[0, 0].legend()

fig.suptitle("Background-subtracted counts and Skellam distributions", fontsize=14)
fig.tight_layout()
plt.show()

この図を見ると、

$\lambda_C=\lambda_B$ なら中心は $0$
$\lambda_C>\lambda_B$ なら正の方向にずれる
$\lambda_C<\lambda_B$ なら負の方向にずれる
カウント数がある程度大きいと Gaussian 近似とかなり近い

ことがわかります。

背景が強い場合、背景差し引き後の値が負になることも自然に起こります。

これは「天体からの信号が負」という意味ではありません。

あくまで、source region で観測された total counts と background estimate の差が、統計的な揺らぎによって負になることがある、という意味です。

低カウントではどうなるか

次に、低カウントの場合を見ます。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import skellam, norm

examples = [
    (0.2, 0.2),
    (1.0, 1.0),
    (1.0, 5.0),
]

for lam_c, lam_b in examples:
    mean = lam_c - lam_b
    sigma = np.sqrt(lam_c + lam_b)

    k_min = int(np.floor(mean - 6 * sigma)) - 2
    k_max = int(np.ceil(mean + 6 * sigma)) + 2
    k = np.arange(k_min, k_max + 1)

    p_skellam = skellam.pmf(k, mu1=lam_c, mu2=lam_b)
    p_gauss = norm.pdf(k, loc=mean, scale=sigma)

    diff = p_skellam - p_gauss
    max_abs_diff = np.max(np.abs(diff))

    fig, (ax1, ax2) = plt.subplots(
        2,
        1,
        figsize=(7, 5.5),
        sharex=True,
        gridspec_kw={"height_ratios": [3, 1], "hspace": 0.0},
    )

    # ======================================================
    # Upper panel: distribution
    # ======================================================
    ax1.vlines(k, 0, p_skellam, linewidth=2, label="Skellam")
    ax1.plot(k, p_skellam, "o", markersize=5)
    ax1.plot(k, p_gauss, "--", label="Gaussian approx.")
    ax1.axvline(mean, linestyle=":", alpha=0.8)

    ax1.set_ylabel("Probability")
    ax1.set_title(
        rf"$\lambda_C={lam_c},\ \lambda_B={lam_b}$"
        + rf", mean={mean:.1f}, $\sigma={sigma:.2f}$"
    )
    ax1.grid(alpha=0.3)
    ax1.legend()

    # Hide x tick labels on upper panel
    ax1.tick_params(labelbottom=False)

    # ======================================================
    # Lower panel: difference
    # ======================================================
    ax2.axhline(0, linestyle=":", alpha=0.8)
    ax2.vlines(k, 0, diff, linewidth=2)
    ax2.plot(k, diff, "o", markersize=4)

    ax2.set_xlabel(r"$C-\hat{B}$")
    ax2.set_ylabel(r"$\Delta p$")
    ax2.grid(alpha=0.3)

    ax2.text(
        0.98,
        0.88,
        rf"$p_{{\rm Skellam}} - p_{{\rm Gaussian}}$"
        + "\n"
        + rf"max $|\Delta p|$={max_abs_diff:.3g}",
        transform=ax2.transAxes,
        ha="right",
        va="top",
        fontsize=9,
    )

    fig.suptitle(
        "Low-count background subtraction",
        fontsize=13,
        y=0.98,
    )

    fig.tight_layout()
    plt.show()

このように、低カウントでは、

Skellam 分布が Gaussian 近似からずれやすい
同じ net counts でも、もとの observed counts と background counts によって分布の広がりが変わる
背景差し引き後の平均値だけを見ると、その不確かさの違いを見落としやすい

という性質が見えてきます。以下でもう少し詳しく説明します。

次に、この最後の点についてもう少し詳しく見ていきます。

背景を引いた量だけを見ると、誤差を直感しにくい

ここまでで、ポアソン分布に従うカウントを引くと Skellam 分布になることを見ました。

ただ、この記事で一番強調したいのは、分布の名前そのものではありません。

重要なのは、

背景を引くと平均値は小さくなるが、不確かさは background の揺らぎも背負ったまま残る

という点です。

背景を引くこと自体は、とても自然な処理です。

私たちが知りたいのは天体由来の信号なので、

N_{\rm net}
=
C
-
\hat{B}

のように、source region で観測された total counts $C$ から、推定した background counts $\hat{B}$ を引きたくなります。

この操作によって、平均値としては background の分が引かれます。

しかし、background の揺らぎまで消えるわけではありません。

誤差伝搬で考えると、

\sigma_{\rm net}^2
=
\sigma_C^2
+
\sigma_{\hat{B}}^2

です。

ポアソン統計では、ざっくり

\sigma_{\rm net}
\sim
\sqrt{C+\hat{B}}

のようになります。

つまり、背景を引いた後の値 $N_{\rm net}$ は小さく見えても、その不確かさは観測された total counts と background estimate の両方に依存しています。

ここが少し直感とずれやすいところです。

たとえば、差し引き後の値が同じ

N_{\rm net}
=
C-\hat{B}
=
2

だったとしても、

C=12,\quad \hat{B}=10

の場合と、

C=102,\quad \hat{B}=100

の場合では、背負っている揺らぎはかなり違います。

前者では、

\sigma_{\rm net}
\sim
\sqrt{12+10}

ですが、後者では、

\sigma_{\rm net}
\sim
\sqrt{102+100}

になります。

どちらも net counts は $2$ です。

しかし、後者の方が background レベルが大きいため、不確かさも大きくなります。

そのため、背景差し引き後のスペクトルや画像だけを見ると、

値としては同じくらいに見えるけれど、どれくらい不確かなのかは background レベルに強く依存している

という情報が見えにくくなります。

これは、可視化を見る側にとって少しミスリーディングになることがあります。

背景を引いたスペクトルは、見た目には「天体信号だけを残したもの」のように見えます。

しかし実際には、観測された total counts と background estimate の両方の揺らぎを背負った量です。

したがって、背景差し引き後の値だけを見るのではなく、もとの total counts や background counts がどれくらいあったのかも一緒に意識しておくと、誤差の大きさを理解しやすくなります。

この意味で、X線天文学では、背景を引いたスペクトルだけで考えるのではなく、観測カウントと background をそれぞれの統計的な量として扱うことが重要になります。

たとえば実際のスペクトル解析では、background-subtracted spectrum を確認用に見ることはあります。一方で、統計的なフィットでは、background を単純に差し引いた量をポアソンデータとして扱うよりも、source spectrum と background spectrum をそれぞれのカウントデータとして likelihood に入れて扱うことが多いです。

このあたりの考え方については、XSPEC manual の Appendix B: Statistics in XSPEC が参考になります。

注意：実際の background scaling について

ここまででは、説明を簡単にするために、

N_{\rm net}
=
C-\hat{B}

として、source region で観測された total counts $C$ から background estimate $\hat{B}$ をそのまま引く形で書きました。

ただし、実際のスペクトル解析では、background region の面積や exposure の違いを考慮して、background counts を source region 相当にスケールして扱います。

たとえば、background region で測ったカウントを $B_{\rm raw}$、スケール係数を $\alpha$ とすると、

\hat{B}
=
\alpha B_{\rm raw}

のように書けます。

このとき、

N_{\rm net}
=
C
-
\alpha B_{\rm raw}

です。

実際、XSPEC では spectrum file の EXPOSURE や BACKSCAL keyword などを用いて、source spectrum と background spectrum の正規化が行われます（参考: On the BACKSCAL keyword - HEASARC）。BACKSCAL は、spectrum を抽出した detector area に対応する量として使われます。

この場合の分散は、

\mathrm{Var}(N_{\rm net})
\sim
C
+
\alpha^2 B_{\rm raw}

のように、background 側の揺らぎもスケール係数つきで効いてきます。

この記事では、まず「ポアソン分布の差はポアソン分布ではない」という直感を見るために、スケール係数を省いた単純な場合を扱いました。

まとめ

この記事では、X線天文学の背景差し引きを出発点にして、ポアソン分布に従うカウント同士を引くと何が起こるのかを見てみました。

source region で観測される total counts には、天体由来の信号だけでなく background も含まれています。そのため、background estimate を引いて net counts を考えるのは自然な処理です。

しかし、ポアソン分布に従うカウント同士を引くと、その結果はもうポアソン分布ではなく、Skellam 分布になります。特に重要なのは、背景を引くことで平均値は小さくなっても、background の揺らぎまで消えるわけではないという点です。

つまり、背景差し引き後の値だけを見ると、どれくらいの不確かさを背負っているのかが見えにくくなることがあります。同じ net counts でも、もとの observed counts と background counts が大きければ、不確かさも大きくなります。

この意味で、X線天文学では、背景を引いた量だけを見るのではなく、もとの observed counts と background counts を意識しながら、ポアソン統計の枠組みで扱うことが大切になります。

この記事が、背景差し引きとポアソン統計の関係を直感的に理解する一助になれば幸いです。

参考