Pythonで複数ピークを自動検出する方法と原理（find_peaksの仕組み・plateau対策・内部実装）

はじめに

Pythonで一次元データから複数ピークを自動検出するニーズは非常に多くあります。

• スペクトル解析で輝線・吸収線を検出したい
• FFT後の周波数スペクトルの共振ピークを拾いたい
• 時系列データの「山」を自動的に検出したい
• ノイズの中から意味あるイベントだけを抽出したい

しかし、

• scipy.signal.find_peaks の パラメータの意味がよくわからない
• plateau（同値ピーク）やノイズで 正しくピークが検出されない
• argrelmax / argrelextrema と find_peaks の 違いが理解できない

といった悩みを持つことも多いでしょう。

本記事では、

✅ Pythonでピーク検出をする際の 基本原理
✅ 差分ベースピーク検出の 落とし穴
✅ SciPy の find_peaks / argrelmax の 内部実装と違い
✅ plateau（同値ピーク）の扱い方
✅ 実用的な パラメータ調整方法と用途別設定例

を、図・コード・FAQ付きで解説します。

python を用いて複数個の peak を自動検出について方法だけでなく、
その原理も含めて解説してみます。

まず、python で、複数個の peak を自動検出する方法の代表格は、

• find_peaks 同値ピークが存在する場合でも正しくピーク検出できる。
• argrelmax 単独の極大値を返す。(同値のピークがある場合は取り逃がす)
• argrelmin 単独の極小値を探す。(同値のピークがある場合は取り逃がす。)
• argrelextrema 任意の判定条件を用いて、極値が判定できる。

の４つの関数で、この中から目的に合致するものをチョイスすることが多いです。

これらの使い方を解説をした記事は色々と出てますので、そちらを参考ください。

ここでは、そもそもどういうパラメータがあって、どういう原理を使ってピークを捉えるのがよいか、について少し説明をしつつ、scipy の中ではどう実装されているかについて解説したいと思います。

複数個のピーク検出の原理について

複数個のピークとは何か？

まずは、ピークは何かを定義しないと、話が噛み合いませんので、この記事ではピークの定義から考えます。実用例から考えてみますと、ピークを検出する必要があるのは、一次元の光のスペクトルの中の輝線や吸収線を見つけたい場合、FFTした周波数信号の共振ピークを見つけたい場合、時系列データの山や谷を捉えたい場合など、でしょうか。この記事ではそのようなピークの検出を想定しています。

これを、コンピュータが解釈できる言葉にすると、ある局所的な範囲において、データが単調に増加し、単調に減少する、その節目のことをピークと呼んでいると思います。数学的には、微係数がゼロ、つまり極大値、極小値を求めるアルゴリズムがあればよいことになります。もしデータが連続的であれば、数学的には必ず微係数がゼロ点がどこかにあるはずですが、現実のコンピュータの世界は離散的にしか数字を扱えないので、ゼロではなく、微係数の符号が変わる点を探すことになります。

つまり、ピークとは、微係数の符号が変わる点を探せばよい、というわけで、もしそれだけで十分な状況であれば、np.diff で微分して符号を調べるだけでよいです。数学的には、連続関数の微分係数が 0 になる点＝極値ですが、デジタルデータは離散値なので 差分の符号変化を見てピークを検出します。

# 微分のイメージ：np.diff で差分をとる
sdiff = np.diff(data)

• 差分が正 → 増加
• 差分が負 → 減少
• 符号が正→負に変わる点 → 極大値

ピークの特徴量

現実的にはもう少し細かくピークを弁別したい状況が多いと思います。ピークが何個もある場合や、見誤りやすいノイズが含まれてる場合など、そう言う場合には、ピークの特徴量を使います。

ピークを特徴付ける量は、

• ピークの高さ
• ピークの位置
• ピークの幅、非対称度
• 一つのピークと前後のピークの間隔

があります。これらを用いて、ピークの選別をします。find_peaks には、あらかじめ用意された選別条件が書かれていて、それについては後述します。ピークの検出さえ正しくできれば、カット条件は自分でコーディングするというスタンスもありかと思います。

ピーク判定後の処理ではなくて、前処理として、低周波や高周波成分を落としたり、サンプリングレートを変える、フィルタの種類を変えるなど、特徴量は前処理にも強く依存しますので、自分のみたい "ピーク" のスケールは環境に合わせたカスタマイズが必要になることが多いかと思います。

差分の符号だけ使った簡単なピーク検出

まずは、もっとも基本の型となる、差分を用いたピーク検出法から紹介します。

に、sin と 対数正規分布の場合の例を紹介しました。

最後の極大値を取得するフィルター生成が肝で、これ図示すると次のような感じです。

差分をとって、１つずらして、それと掛け算で符号を見ると、極値のところだけ符号が変わる。さらに、極大値と極小値を弁別するために、差分の正負で場合分けをする。

sin の場合

import numpy as np
import matplotlib.pyplot as plt

# sin 関数の場合
smin, smax, slen = 0, 15, 100
x = np.linspace(smin, smax, slen)
s = np.sin(x)
plt.plot(x,s,"o")
plt.show()

sdiff = np.diff(s) # 差分を取る。(要素数が１つ減る)

# 極大値のみを取得 : (差分 x 1つシフトした差分が負) & (差分が正)、という条件
sdiff_sign = ((sdiff[:-1] * sdiff[1:]) < 0) & (sdiff[:-1] > 0)

plt.plot(x, s, label="sin curve")
plt.plot(x[1:-1][sdiff_sign], s[1:-1][sdiff_sign],"o",label="peak")
plt.legend()
plt.show()

対数正規分布の場合

を参考に、対数正規分布に従う乱数を生成しています。

rng = np.random.default_rng()
mu, sigma, nsample, nbins =3., 0.1, 500, 100
s = rng.lognormal(mu, sigma, int(nsample))

# 対数正規分布で生成した乱数のヒストグラムを表示
count, bins, ignored = plt.hist(s, bins = nbins, density=True, align='mid')

# 生成に用いた対数正規分布のグラフをプロットする
x = np.linspace(min(bins), max(bins), 10000)
pdf = (np.exp(-(np.log(x) - mu)**2 / (2 * sigma**2)) / (x * sigma * np.sqrt(2 * np.pi)))
plt.plot(x, pdf, linewidth=2, color='r')
plt.yscale("log")
plt.axis('tight')
plt.show() # ヒストグラムで、対数正規分布の生成を確認

x = np.arange(len(s))
sdiff = np.diff(s)
sdiff_sign = ((sdiff[:-1] * sdiff[1:]) < 0) & (sdiff[:-1] > 0)
plt.plot(x, s, label="log-normal curve")
plt.plot(x[1:-1][sdiff_sign], s[1:-1][sdiff_sign],"o",label="peak")
plt.legend()
plt.show() # peak判定の結果を確認

生成した対数正規分布のヒストグラムと、極大値の判定結果です。

簡単なピークサーチの落とし穴

実は、この差分を使うピークサーチは簡単なようで、１つ落とし穴がある。実は、同じ値が連続した時に判定されないのである。それを図説したのがこちら。

差分は、全く同じ値が連続してしまうと、0 になってしまう。0 は、どんな有限な値との掛け算でも 0 になってしまうので、判定の時に漏れてしまう。0 を正しく扱うには、区別がつかない代物なので、特別扱いしないといけない。

極大値に同値がある場合の差分によるピーク検出の例

実際に、sin波の一部をちょんぎった場合の例がこちら。

smin, smax, slen = 0, 30, 100
x = np.linspace(smin, smax, slen)
s = np.sin(x)
s[(s>0.8) & (x<10)] = 0.8
sdiff = np.diff(s)
sdiff_sign = ((sdiff[:-1] * sdiff[1:]) < 0) & (sdiff[:-1] > 0)
plt.plot(x, s, ".--", label="chopped sine")
plt.plot(x[1:-1][sdiff_sign], s[1:-1][sdiff_sign],"o",label="peak")
plt.legend()
plt.show()

scipy の実装例

argrelmax などの差分タイプ

argrelmax の御本尊にあたる部分は、_boolrelextrema という関数を用いている。

それとほぼ同等なコードがこちらに書きくだしたものです。

# argrelmax の御本尊の関数を使う
def _boolrelextrema_local(data, comparator, axis=0, order=1, mode='clip', debug=True):
# using def _boolrelextrema
# https://github.com/scipy/scipy/blob/c1ed5ece8ffbf05356a22a8106affcd11bd3aee0/scipy/signal/_peak_finding.py#L22
  datalen = data.shape[axis]
  locs = np.arange(0, datalen)
  results = np.ones(data.shape, dtype=bool)
  main = data.take(locs, axis=axis, mode=mode)

  for shift in range(1, order + 1):
    if debug: print("shift = ", shift)
    plus = data.take(locs + shift, axis=axis, mode=mode)
    minus = data.take(locs - shift, axis=axis, mode=mode)
    results &= comparator(main, plus)
    results &= comparator(main, minus)
    if debug: 
      print("plus = ", plus)
      print("minus = ", minus)
      print("res1 = ", results)
      print("res2 = ", results)
    if ~results.any(): # 全部 False になったら更新はないので、ループから抜ける。
      return results
  return results    

• データを同じサイズの results という True/False を詰める配列を用意して、numpy.take で配列をシフトし、その時に mode=clipなので、はみ出たアクセスは先頭と最後の値を詰めるだけ。
• comparetor としては、numpy.greater などの比較演算子を入れることで、True/Falseを更新し、results &= comparator(main,plus) などで、新しい比較結果との AND で更新していく。
• 最後まで True が連続し続けたら True が残り、一度でも False がでると False になるので、results 全てが False になると for loop が終了する。

という感じで、配列をグルグル回して、True/Falseの配列に AND 演算を重ねて更新していく、というスタンスです。

# argrelmax と _boolrelextrema でピーク判定してみる。
smin, smax, slen = 0, 30, 100
x = np.linspace(smin, smax, slen)
data = np.sin(x)
data[(data>0.7) & (x<10)] = 0.7
comparator=np.greater
axis=0
order=5
mode='clip'
debug=False

id1 = argrelmax(data, axis=axis, order=order, mode=mode)
id2 = argrelextrema(data, comparator, axis=axis, order=order, mode=mode)
id3 = _boolrelextrema_local(data, comparator, axis=axis, order=order, mode=mode, debug=debug)

plt.figure(figsize=(5,3))
plt.plot(x, data, ".--", label="chopped sine")
plt.plot(x[id1], data[id1],"o", label="peak using argrelmax")
plt.plot(x[id2], data[id2],"x", label="peak using argrelextrema")
plt.plot(x[id3], data[id3],",", label="peak using boolrelextrema")


plt.legend(loc="lower right")
plt.show()

scipy の argrelmax、argrelextrema と、上で自前で定義した _boolrelextrema_local を用いてピーク判定したものはこちらです。

結果は同じで、どれもピークが同値の場合は判定されない。

find_peaks を用いた場合

scipy の find_peaks を用いると、plateau_size=1 を指定すると、プラトー (同値のピーク) の サイズと、左端と右端の配列を返してくれるので、この配列を用いれば、事後処理でプラトーの条件によりイベント弁別ができる。

# argrelmax と _boolrelextrema でピーク判定してみる。
smin, smax, slen = 0, 30, 100
x = np.linspace(smin, smax, slen)
data = np.sin(x)
data[(data>0.7) & (x<10)] = 0.7

# from scipy.signal import argrelmax
peaks_plateau1, properties_plateau1 = find_peaks(data, plateau_size=1)
peaks_plateau2, properties_plateau2 = find_peaks(data, plateau_size=2)
peaks_plateau10, properties_plateau10 = find_peaks(data, plateau_size=10)

plt.figure(figsize=(5,3))
plt.plot(x, data, ".--", label="chopped sine")
plt.plot(x[peaks_plateau1], data[peaks_plateau1],"x",label="peak using find_peaks plateau_size=1")
plt.legend(loc="lower right")
print(properties_plateau1)
plt.show()

plt.figure(figsize=(5,3))
plt.plot(x, data, ".--", label="chopped sine")
plt.plot(x[peaks_plateau2], data[peaks_plateau2],"x",label="peak using find_peaks plateau_size=2")
plt.legend(loc="lower right")
print(properties_plateau2)
plt.show()

plt.figure(figsize=(5,3))
plt.plot(x, data, ".--", label="chopped sine")
plt.plot(x[peaks_plateau10], data[peaks_plateau10],"x",label="peak using find_peaks plateau_size=10")
plt.legend(loc="lower right")
print(properties_plateau10)
plt.show()

プラートが偶数個のデータ点の場合

プラトーが偶数個ある場合はどうなるでしょうか？　それを実際に見てみましょう。

# argrelmax と _boolrelextrema でピーク判定してみる。
data = np.array([0,1,1,1,0,0,1,1,1,1,0,0,0,1,0,0,1,1,1,1,1,0,0])
x = np.arange(len(data))

# from scipy.signal import argrelmax
peaks_plateau1, properties_plateau1 = find_peaks(data, plateau_size=1)
peaks_plateau2, properties_plateau2 = find_peaks(data, plateau_size=2)
peaks_plateau10, properties_plateau10 = find_peaks(data, plateau_size=10)

plt.figure(figsize=(5,3))
plt.plot(x, data, "o-", label="chopped sine")
plt.plot(x[peaks_plateau1], data[peaks_plateau1],"x",label="peak using find_peaks plateau_size=1")
plt.legend(loc="lower right")
print(properties_plateau1)
plt.show()

plt.figure(figsize=(5,3))
plt.plot(x, data, "o-", label="chopped sine")
plt.plot(x[peaks_plateau2], data[peaks_plateau2],"x",label="peak using find_peaks plateau_size=2")
plt.legend(loc="lower right")
print(properties_plateau2)
plt.show()

plt.figure(figsize=(5,3))
plt.plot(x, data, "o-", label="chopped sine")
plt.plot(x[peaks_plateau10], data[peaks_plateau10],"x",label="peak using find_peaks plateau_size=10")
plt.legend(loc="lower right")
print(properties_plateau10)
plt.show()

これを実行した結果が、

となりますが、偶数個の場合は、前の方のデータ点をピークとして選んでいることがわかります。

find_peaks の中身の詳細

find_peaks の _peak_finding.py はインターフェースだけで、実体は、高速化のために、

の中で cython で実装された _local_maxima_1d が御本尊である。cython については、

などを参照ください。

def _local_maxima_1d(const np.float64_t[::1] x not None): の中で、for loop でグルグ回して、極大を探して、極大値が同じなら、ピークの左端(left_edges)と右端(right_edges)を計算する、ということをやってます。計算のメイン部分はこちらで、

_peak_finding_utils.pyx

def _local_maxima_1d(const np.float64_t[::1] x not None):
    """
    Find local maxima in a 1D array.
    This function finds all local maxima in a 1D array and returns the indices
    for their edges and midpoints (rounded down for even plateau sizes).
    """ 
    #(省略)
    with nogil:
        i = 1  # Pointer to current sample, first one can't be maxima
        i_max = x.shape[0] - 1  # Last sample can't be maxima
        while i < i_max:
            # Test if previous sample is smaller
            if x[i - 1] < x[i]:
                i_ahead = i + 1  # Index to look ahead of current sample

                # Find next sample that is unequal to x[i]
                while i_ahead < i_max and x[i_ahead] == x[i]:
                    i_ahead += 1

                # Maxima is found if next unequal sample is smaller than x[i]
                if x[i_ahead] < x[i]:
                    left_edges[m] = i
                    right_edges[m] = i_ahead - 1
                    midpoints[m] = (left_edges[m] + right_edges[m]) // 2
                    m += 1
                    # Skip samples that can't be maximum
                    i = i_ahead
            i += 1

その中で、ピークの値(midpionts)は、

midpoints[m] = (left_edges[m] + right_edges[m]) // 2

で、 left_edges と right_edges の中間を //2 で求めていますが、奇数の場合は //2 は切り捨てなので、e.g., 7//2 = 3 となります。

find_peaks のパラメータの調整

より細かなピークの判定においては、find_peaks のオプションのパラメータについても調整が必要となります。

# 対数正規分布の生成
# https://www.headboost.jp/generator-lognormal/ を参考に、対数正規分布に従う乱数を生成する。
rng = np.random.default_rng()
mu, sigma, nsample =1., 0.5, 1e5
signal = rng.lognormal(mu, sigma, int(nsample))

time = np.arange(len(signal))
plt.plot(time, signal, label="log-normal curve")
plt.legend()
plt.show()

distance=20
peaktime_list_d20, plateau_list_d20 = find_peaks(signal, plateau_size=1, distance = distance) 
distance=40
peaktime_list_d40, plateau_list_d40 = find_peaks(signal, plateau_size=1, distance = distance) 
print("plateau_sizes = ", set(plateau_list["plateau_sizes"])) # check plateau_sizes to know the number of the flat peaks 

diff_peaktime_list_d20 = np.diff(peaktime_list_d20)
diff_peaktime_list_d40 = np.diff(peaktime_list_d40)
w=1
print(plt.hist(diff_peaktime_list_d20, bins=np.arange(1, max(diff_peaktime_list_d40) + w, w), label = "distance = 20", alpha=0.8))
print(plt.hist(diff_peaktime_list_d40, bins=np.arange(1, max(diff_peaktime_list_d40) + w, w), label = "distance = 40", alpha=0.8))
plt.legend()
plt.yscale("log")
plt.xlabel("distance between two adjacent peaks")

この結果は、

のように、discance の違い、ここでは、20, 40 の違いにより、ピークの時間差の分布が異なることがわかる。

ピークの分布についても、

w=0.4
print(plt.hist(signal[peaktime_list_d20], bins=np.arange(1, max(signal[peaktime_list_d40]) + w, w), label = "distance = 20", alpha=0.6))
print(plt.hist(signal[peaktime_list_d40], bins=np.arange(1, max(signal[peaktime_list_d40]) + w, w), label = "distance = 40", alpha=0.6))
plt.legend()
plt.yscale("log")
plt.xlabel("peaks of signals")

のように、cut 条件によって、ピークの分布が異なる。この場合は、discance 20 だと小さいピークをたくさん引っ掛けていることがわかる。

この他にも、

• height
  • ピークの高さの絶対値の判定基準
• threshold
  • ピークの高さの判定条件だが、height と異なり、隣のピークに比べてどのくらい高いか、という基準。
• distance
  • ２つの連続するピークからの距離
• prominence
  • wlen でピークの周囲の探索幅を指定し、その幅内で、ピークの左右の最小値を探す max(left_min, right_min)
  • prominences は、ピークの左右の最小値の大きい方との差で定義される
    • prominences[peak_nr] = x[peak] - max(left_min, right_min)
  • https://github.com/scipy/scipy/blob/v1.10.1/scipy/signal/_peak_finding_utils.pyx
• width
  • ピークの高さを、peak -　prominences x rel_heigth で定義する。
    • height = width_heights[p] = x[peak] - prominences[p] * rel_height
  • この高さまで、信号が下がる幅を計算する。
  • ある window 範囲内で幅が計算できなければ、直線で補完して幅を定義する。
  • https://github.com/scipy/scipy/blob/v1.10.1/scipy/signal/_peak_finding_utils.pyx
• wlen
  • peaks prominences を計算するための幅
• rel_height
  • peaks width を計算するための相対的な高さ
• plateau_size
  • 連続したピークの数

などの条件があるので、問題に応じてこれらのオプションを適切に使うか、自分でオフラインで処理をするとよいです。

【応用編】 ピークの集積方法

応用例として、ピークを集積する方法を紹介します。 find_peaks でピークを検出して、そのピークを stacking して平均的なプロファイルを取得する方法です。

対数正規分布の時系列データの生成

サンプルとして、対数正規分布の時系列データを用意します。

# https://www.headboost.jp/generator-lognormal/ を参考に、対数正規分布に従う乱数を生成する。
rng = np.random.default_rng()
mu, sigma, nsample =3., 2, 8192
signal = rng.lognormal(mu, sigma, int(nsample))

time = np.arange(len(signal))
plt.plot(time, signal, label="log-normal curve")
plt.legend()
plt.show()

など、お好みのパラメータと長さで用意します。

ピークの検出

peak の 判定を行い、ピークを ID できるように peaktime_list というピークの場所が詰まった配列を find_peaks で取得します。

同時ピークの扱いに備えて、plateau_size=1 というオプションで find_peaks を実行しておきます。こうすることで、set(plateau_list["plateau_sizes"]) で、"plateau_sizes" が 1 のみであれば、すべて peak が一点のみだとわかります。 それ以外であれば、複数の同時ピークが存在しうるので、問題に応じて取捨選択を考えましょう。

d = 20
peaktime_list, plateau_list = find_peaks(signal, plateau_size=1, distance=20) 
print("plateau_sizes = ", set(plateau_list["plateau_sizes"])) # check plateau_sizes to know the number of the flat peaks 

ここでは、distance = 20 とした。このようなパラメータに依存することには注意して解析しましょう。

ピーク前後のデータを集積する

ピークの前後 +/- order 個のデータを切り出して、足しこみます。
np.take(signal, np.arange(peaktime-order,peaktime+order+1)) で、ピークの前後を切り出して、それを append して追加していきます。time の方は、offset (time[peaktime]のこと)を引いておくことで、peaktime を offset にした切り出しになります。

times_list, peaks_list,  = [], []
order = 20 # change as you like
debug = False
for i, peaktime in enumerate(peaktime_list):

  # check out of range 
  if ((peaktime - order) < 0) or ((peaktime + (order + 1)) > len(peaktime_list)): # skip when the range is out of input the array
    continue

  # get signal array
  stmp = np.take(signal, np.arange(peaktime-order,peaktime+order+1))  
  peaks_list.append(stmp)

  # get time array from the peaktime 
  ttmp = np.take(time, np.arange(peaktime-order,peaktime+order+1)) - time[peaktime] 
  times_list.append(ttmp)

  # plot if needed to check
  if debug: 
    plt.plot(times_list[-1], peaks_list[-1],"o")
    plt.legend()
    plt.show()    

times_list = np.array(times_list)
peaks_list = np.array(peaks_list)

np.arange(peaktime - order, peaktime + order + 1) で、 order が偶数の場合は、奇数個のサンプルの取得になる。 e.g.,

In [7]: np.arange(100-2,100+2+1)
Out[7]: array([ 98, 99, 100, 101, 102])

となります。

集積したピークのプロットと平均プロファイルの作成

最後に、np.mean で切り出したピーク近傍の配列を平均化します。

# plot every profile to check (<100)
for i, (onet, onep) in enumerate(zip(times_list, peaks_list)):
  if i < 100:
    plt.plot(onet, onep + i, ".-") # plot every profile with offset i 
plt.show()

このように、100個分の profile を表示した例ですが、変なデータ混じっていないか、データをよく見ておきましょう。
OKそうであれば、平均化をしてみましょう。

# 時間の平均 (全部同じなので平均は不要なはず)
folded_time = np.mean(times_list, axis=0)
# ピークの平均化
folded_peak = np.mean(peaks_list, axis=0)

# プロファイルのプロット (ブラックホール連星の解析では、shot profile と呼ばれるもの。)
plt.plot(folded_time, folded_peak, "o-")
plt.show()

これで生成された図が、

で、これが対数正規分布の時系列データのピークの頭を揃えて、足し込んだプロファイルになります。BH連星のデータ解析では shot profile などと呼ばれるものです。

ショット解析の参考文献

手前味噌ですが、、BH連星 Cyg X-1 へのショット解析を使った例は、こういうのがあります。

さて、ここまででは、find_peaks で検出したピークの時刻をそろえ、その前後のデータを平均しました。

このような解析は、ブラックホール連星などの時間変動解析では、一般にショット解析、ショット重ね合わせ解析、あるいは superposed shot analysis と呼ばれます。

ただし、ここで一つ、非常に重要な注意があります。

ショット解析でピークを重ね合わせたからといって、観測された時系列が、互いに独立した同一形状の「ショット」の足し合わせで生成されていると証明されたわけではありません。

この点を理解するには、ショット解析が数学的に何を計算しているのかを考える必要があります。

検出されたピーク時刻を

t_1,~t_2,~\ldots,~t_N

とします。それぞれのピークを時刻ゼロにそろえ、ピークからの相対時刻を $\tau$ とすると、平均ショットプロファイルは概念的には

\bar{x}_{\rm shot}(\tau)
=
\frac{1}{N}
\sum_{i=1}^{N}
x(t_i+\tau)

です。

局所的な平均値 $b_i$ を差し引く場合は、

\bar{x}_{\rm shot}(\tau)
=
\frac{1}{N}
\sum_{i=1}^{N}
\left[
x(t_i+\tau)-b_i
\right]

となります。

これは何を意味するのでしょうか？

ショット解析が直接測っているのは、

「指定したピーク検出条件を満たす時刻を選んだとき、その前後のデータが平均的にどのような形をしているか」

です。

統計的には、ピークが検出されたという条件のもとでの条件付き平均と考えるのが適切です。

\bar{x}_{\rm shot}(\tau)
\sim
E
\left[
x(t+\tau)
\mid
t~\text{がピーク検出条件を満たす}
\right]

したがって、得られた平均プロファイルは、元の天体が持つ時間変動だけでなく、

• ピークの定義
• 検出に用いたエネルギーバンド
• 時間ビン幅
• height
• distance
• prominence
• 局所平均の定義
• 平滑化の方法
• ピーク前後の切り出し幅

にも依存します。

ショット解析では、解析条件そのものが、得られる平均波形の一部を決めていることに注意が必要です。

---

ショット解析と「ショットノイズモデル」は同じではない

ここは用語が混乱しやすいところです。

ショット解析

観測データからピークを選び、その周囲を重ね合わせる解析手法です。

ショットノイズモデル

時系列を、何らかの基本波形のランダムな重ね合わせとして表す物理モデルまたは確率モデルです。概念的には、

x(t)
=
b(t)
+
\sum_i
A_i
s_i(t-t_i)
+
n(t)

のように考えます。

ここで、

• $b(t)$：ゆっくり変化するベースライン
• $A_i$：各ショットの振幅
• $s_i(t)$：各ショットの波形
• $t_i$：発生時刻
• $n(t)$：観測ノイズ

です。

ショット解析は、このようなモデルが正しい場合には、基本波形 $s_i(t)$ の平均的な性質を調べる強力な方法になります。

しかし、実際の時系列が、

• 伝播する降着率揺らぎ
• 複数の時間スケールを持つ連続的な赤色雑音
• 非線形な変動
• 複数成分の重ね合わせ

によって生じていても、局所最大値を選んで平均すれば、何らかの「山型の平均プロファイル」は得られます。

したがって、

平均ショットプロファイルが得られたことと、物理的に独立したショットが存在することは、同じではありません。

ショット解析は、物理モデルを仮定しなくても使える解析手法です。その一方で、得られたプロファイルを特定の物理現象と結びつけるには、追加の検証が必要です。

---

なぜピークを重ねると、意味のある情報が得られるのか？

一つ一つのピークは、光子統計や周囲の変動に埋もれているため、詳細な形状を調べることが難しい場合があります。

そこで、似た条件を満たすピークを $N$ 個重ね合わせます。

各イベントに共通する成分が存在し、イベントごとのランダムな揺らぎが互いに独立であれば、ランダム成分はおおよそ

\frac{1}{\sqrt{N}}

に比例して小さくなります。

このため、一つのイベントでは見えない、

• 立ち上がりと立ち下がりの非対称性
• エネルギーごとのピーク時刻の違い
• 硬度比の時間変化
• スペクトル形状の時間発展
• 反射成分や吸収成分の応答

などを抽出できる可能性があります。

特に、ピーク時刻を一つの高統計エネルギーバンドで決め、同じ時刻を別のエネルギーバンドに適用すれば、低統計のエネルギーバンドでも平均的な時間変化を調べられます。

これは、ショット解析の大きな利点です。

---

Cyg X-1における古典的なショット解析

Cyg X-1のショット解析における代表的な研究の一つが、Negoro et al. (1994) です。

Negoro et al. は、Gingaで観測したCyg X-1のlow stateのライトカーブから多数のピークを選び、ピーク時刻をそろえて重ね合わせました。

その結果、平均的な強度変化には、

• 数秒程度にわたる、ほぼ時間対称な立ち上がりと立ち下がり
• 約 $0.1~{\rm s}$ と約 $1~{\rm s}$ の二つの時間スケール
• ピークへ向かって徐々に軟化し、ピーク通過後に急激に硬化するスペクトル変化

が見いだされました。

ここで興味深いのは、強度プロファイルがほぼ時間対称である一方、スペクトル変化は時間非対称だったことです。

単に「明るくなって暗くなる」というだけでなく、立ち上がりと立ち下がりでは、放射領域の物理状態が異なる可能性を示しています。

• H. Negoro, S. Miyamoto, and S. Kitamoto,
  “Structure of X-Ray Shots of Cygnus X-1 in Its Low State,”
  The Astrophysical Journal Letters, 423, L127, 1994.
  NASA ADS

---

Yamada et al. (2013) による硬X線帯域への展開

Yamada et al. (2013) では、Suzakuで観測したlow/hard stateのCyg X-1にショット解析を適用し、HXD-PINとHXD-GSOを用いて、約 $10$–$200~{\rm keV}$ の広いエネルギー帯域でショットプロファイルを調べました。

この解析でも、

• 強度変化はおおむね時間対称
• 硬度変化は時間非対称

という特徴が確認されました。

さらに、ショットの位相ごとにスペクトルを作成し、熱的Comptonizationモデルを用いて解析した結果、

• 電子温度
• 光学的厚さ
• Compton $y$ パラメータ

がショット位相に伴って変化することが示されました。

特に、Compton $y$ パラメータと電子温度はピークへ向かって徐々に低下し、光学的厚さは増加する傾向を示しました。一方、ピーク通過後には、これらの量が短時間で平均状態へ戻るという非対称な変化が得られました。

この例は、ショット解析が単に平均的なライトカーブを作るだけでなく、

秒以下の時間スケールで、放射領域の温度や光学的厚さがどのように変化するかを調べる「位相分解分光」へ拡張できる

ことを示しています。

• S. Yamada et al.,
  “Rapid Spectral Changes of Cygnus X-1 in the Low/Hard State with Suzaku,”
  The Astrophysical Journal Letters, 767, L34, 2013.
  NASA ADS
  arXiv

---

ショット解析は古い手法なのか？

ショット解析は古くから使われてきた手法ですが、過去の手法というわけではありません。

高時間分解能と広いエネルギー帯域を持つ観測装置が利用可能になったことで、近年でも、

• 軟X線と硬X線の同時ショット解析
• ショット位相分解スペクトル
• 降着円盤とComptonization成分の時間変化
• 複数の観測装置をまたいだ同時解析

に使われています。

例えば、Bhargava et al. (2022) は、AstroSatとNICERの同時観測を用いて、Cyg X-1の約 $0.1$–$80~{\rm keV}$ にわたるショットを解析しました。軟X線側まで含めてショット位相分解分光を行い、ショットに伴う降着円盤成分の変化を議論しています。

• Y. Bhargava et al.,
  “Probing the shot behaviour in Cygnus X-1 using simultaneous AstroSat–NICER observation,”
  Monthly Notices of the Royal Astronomical Society, 512, 6067, 2022.
  arXiv

したがって、ショット解析は、Fourier解析やパワースペクトル解析に置き換えられた古い方法ではありません。

むしろ、

• Fourier解析は「どの時間スケールに、どれだけ変動があるか」
• ショット解析は「特徴的なピークの前後で、強度やスペクトルがどう変化するか」

を見る方法であり、両者は相補的です。

---

最も重要な注意：ピークを選べば、中心には必ず山ができる

ピーク時刻を検出し、その時刻をゼロにそろえて平均しているため、$\tau=0$ には原理的に高い値が集まります。

つまり、

平均プロファイルの中央にピークが存在すること自体は、解析によってほぼ保証されています。

白色雑音や赤色雑音だけで作った時系列であっても、局所最大値を選んで重ねれば、中央に山型の構造が現れます。

このため、議論すべきなのは、

• 中央に山があるか
• 平均プロファイルが滑らかか

だけではありません。

本当に調べるべきなのは、

• 立ち上がりと立ち下がりは対称か
• 特徴的な時間スケールがあるか
• エネルギーによって形状が変わるか
• 硬度比やスペクトルが非対称に変化するか
• 同じパワースペクトルを持つ模擬データでも再現されるか
• 検出条件を変えても結果が維持されるか

です。

---

ピーク選択による「勝者の呪い」

観測値には、天体自身の変動だけでなく、Poisson揺らぎや測定ノイズも含まれます。

ピーク検出では、その中から特に大きな値を選びます。そのため、選ばれたピークには、

• 本当に信号が大きかった
• 偶然、正方向の統計揺らぎが加わった

という二つの効果が混ざります。

同じエネルギーバンドを使ってピークを検出し、そのまま同じデータを重ね合わせると、中心ビンの振幅が統計揺らぎによって上方に偏ることがあります。

これは一般に、選択バイアス、あるいは winner's curse に近い問題です。

対策としては、

• 高統計の参照バンドでピークを検出し、独立した別バンドで形状を測る
• 独立した検出器でピークを検出する
• 偶数・奇数イベントなどにデータを分ける
• Poisson雑音を含む模擬ライトカーブで同じ解析を行う

といった方法が考えられます。

ただし、別のエネルギーバンドでピークを検出すると、そのバンドで明るいイベントだけを選択することになります。

したがって、完全に「無偏差」になるわけではなく、何に対して条件付けした平均なのかを明示する必要があります。

---

distance は単なる重複除去ではない

find_peaks の distance は、隣接するピークの間に必要な最小サンプル数を指定します。

例えば、時間ビン幅が

\Delta t
=
0.1~{\rm s}

で、

distance = 20

とすると、ピーク間隔は最低でも

20
\times
0.1~{\rm s}
=
2~{\rm s}

必要です。

ただし、近接したピークが複数ある場合、find_peaks は小さいピークから除外して、条件を満たすようにします。

したがって、distance を大きくすると、

• 独立したピークだけを選べる
• 同じイベント内の細かい揺らぎを除外できる

という利点がある一方で、

• 小振幅イベントが選ばれにくくなる
• ピークが密集する状態を除外してしまう
• 長時間スケールの広いピークを優先する可能性がある
• ショットの発生間隔分布を人工的に変える

という影響があります。

distance は単なる描画上の調整値ではなく、解析対象とするイベントの定義そのものです。

---

局所平均の定義にも注意する

ショットの振幅を、

A_i
=
x(t_i)-b_i

と定義するとき、局所平均 $b_i$ をどのように求めるかによって結果が変わります。

例えば、ピーク前後の区間全体から平均を求めると、ピーク自身が平均値に含まれます。その結果、

• ピーク振幅が小さく評価される
• 幅の広いショットほど強く差し引かれる
• 長い裾野がベースラインとして消される

可能性があります。

逆に、ピークから十分離れた区間だけで平均値を求めると、低周波変動を正しく追跡できない可能性があります。

そのため、

• 移動平均
• 移動中央値
• ピークから離れた左右の区間
• 低次多項式
• ローパスフィルタで推定したトレンド

など、複数の方法を試し、結果が変わらないか確認する必要があります。

特に移動中央値は、鋭いピークの影響を受けにくいため、局所的なベースライン推定に有効な場合があります。

最初の段階では、平均を引かずに解析する、というのが first choise です。

---

ピーク高さで規格化するかどうか

各ショットをそのまま平均すると、振幅の大きいイベントが平均値に強く影響します。

一方、各ショットをピーク高さで規格化してから平均すると、大小すべてのイベントがほぼ同じ重みになります。

この二つは、異なる物理量を測っています。

規格化しない場合

\bar{x}(\tau)
=
\frac{1}{N}
\sum_i
x_i(\tau)

明るいイベントを含めた、平均的な光子強度の変化を表します。

ピーク振幅で規格化する場合

\bar{x}_{\rm norm}(\tau)
=
\frac{1}{N}
\sum_i
\frac{x_i(\tau)-b_i}{x_i(0)-b_i}

振幅によらない、平均的な相対形状を調べやすくなります。

ただし、ピーク振幅が小さいイベントでは、分母に含まれるノイズが大きくなります。また、振幅と時間幅が相関している場合、規格化によってその情報が失われます。

したがって可能であれば、

• 規格化なし
• 局所平均で規格化
• ピーク振幅で規格化
• 振幅ごとにグループ分け

を比較することが望ましいです。

---

重なり合ったショットをどう考えるか？

現実のライトカーブでは、一つのピークの裾に別のピークが重なっていることがあります。

このとき、重ね合わせ後の平均プロファイルは、単一イベントの形だけではなく、

• 周囲に存在する別イベント
• ピークどうしの時間的な相関
• 低周波のベースライン変動

も含みます。

distance を大きくすれば近接ピークを減らせますが、今度は「孤立したイベントだけを選んだ平均」になります。

それが解析目的に合っているかは別問題です。

したがって、

独立した一つのショット形状を知りたいのか、それともピークの周囲に存在する平均的な変動構造を知りたいのか

を区別する必要があります。

後者であれば、周囲のピークを完全に除く必要はありません。しかし、その場合に得られるものを「単一ショットの形」と呼ぶのは慎重であるべきです。

---

時間ビン幅は物理的な解析条件である

時間ビン幅 $\Delta t$ を粗くすると、

• 細いピークが平均化される
• ピーク振幅が低下する
• ピーク時刻が量子化される
• 短時間の非対称性が見えなくなる

という影響があります。

一方、時間ビン幅を細かくしすぎると、一ビンあたりの光子数が減り、Poisson揺らぎによる偽ピークが増えます。

したがって、最適な時間ビン幅は、

• 観測装置の時間分解能
• カウントレート
• 調べたい時間スケール
• dead time
• telemetry saturation
• pile-up

などによって決まります。

時間ビン幅を変えたときに、平均プロファイルの幅や非対称性がどの程度変わるかを確認することが重要です。

---

平滑化してからピークを検出してもよいか？

ノイズが多い場合、ピーク検出前に平滑化したくなります。

平滑化は偽ピークを減らす一方で、

• ピーク時刻をずらす
• ピーク幅を広げる
• 近接ピークを一つにまとめる
• 非対称な形状を対称に近づける

可能性があります。

一つの方法は、

1. 平滑化したライトカーブでピーク時刻だけを検出する
2. プロファイルの切り出しと平均には、平滑化していない元データを使う

ことです。

それでもピーク時刻には平滑化の影響が残るため、平滑化スケールを変えた場合の再現性を確認する必要があります。

---

ショット解析の不確かさをどう評価するか？

平均プロファイルだけを描いても、どの構造が有意なのかは判断できません。

単純には、各相対時刻 $\tau$ におけるイベント間の標準偏差を用いて、

\sigma_{\bar{x}}(\tau)
=
\frac{
s(\tau)
}{
\sqrt{N}
}

と標準誤差を推定できます。

しかし、これは各ショットが互いに独立であると仮定しています。

ピークが時間的に相関している場合や、同じ長い変動の中から複数のピークを選んでいる場合には、実効的な独立イベント数は $N$ より小さくなります。

より慎重には、

• ピーク単位のbootstrap
• 観測区間単位のbootstrap
• 連続した時間ブロックを再標本化するblock bootstrap
• 同じパワースペクトルと光度分布を持つ模擬時系列との比較

などを使います。

エネルギーごとのショットプロファイルや位相分解スペクトルを比較する場合は、時間ビンやエネルギービンを多数比較するため、多重比較の影響にも注意が必要です。

---

最低限行いたいヌルテスト

ショット解析の結果を物理的に解釈する前に、少なくとも次の確認を行うことをおすすめします。

1. ランダム時刻での重ね合わせ

   実際のピーク数と同じ数だけランダムな時刻を選び、同じ処理を行います。

2. 模擬ライトカーブでの解析

   観測データと似たパワースペクトル、平均カウントレート、Poisson雑音を持つ模擬データに、同じピーク検出を適用します。

3. 検出パラメータ依存性

   height、distance、prominence、時間ビン幅を変えても、主要な結論が維持されるか確認します。

4. データ分割

   観測前半と後半、別の観測区間、異なる検出器などに分けても同じ特徴が得られるか確認します。

5. 検出バンド依存性

   軟X線で検出した場合と硬X線で検出した場合を比較します。

6. ピーク時刻のずらしテスト

   検出時刻を意図的に数ビンずらしたとき、中心付近の特徴がどのように変わるか確認します。

これらのテストで残る特徴ほど、単なるピーク選択効果ではなく、天体固有の時間変動を反映している可能性が高くなります。

---

ショット解析で言えること・言えないこと

比較的直接的に言えること

• 選択条件を満たすピークの平均的な時間構造
• 立ち上がりと立ち下がりの非対称性
• エネルギー依存性
• ピーク前後の平均スペクトル変化
• ピークと他の観測量との平均的な時間関係

解析だけでは直ちに言えないこと

• 各ピークが独立した物理イベントである
• すべてのピークが同じ発生機構を持つ
• 平均プロファイルが一つの素過程の応答関数である
• 得られた時定数が、一意に特定の物理時間を表す
• 中央のピーク自体が統計的に有意な新構造である

平均ショットプロファイルの時定数を、

• 光子の拡散時間
• Compton冷却時間
• 降着流の伝播時間
• 力学時間
• 熱時間
• 粘性時間

などと対応させる場合には、その物理モデルが他の観測量とも整合するかを検証する必要があります。

---

Appendix A：find_peaks パラメータ早見表

パラメータ	何を指定するか	どのようなピークを選ぶか	主な注意
height	ピーク値の絶対的な範囲	一定値より高いピーク	ベースラインが変動すると選択基準も実質的に変わる
threshold	ピークと直近の左右サンプルとの縦方向の差	隣接点より急に突き出したピーク	隣の「ピーク」との差ではない。時間ビン幅やノイズに敏感
distance	ピーク間の最小距離	近接ピークを除いたピーク	単位はサンプル数。近接時は小さいピークから除外される
prominence	周囲の谷を基準にしたピークの目立ち方	局所ベースラインから明確に突出したピーク	wlen により評価する周囲の範囲が変わる
width	指定した相対高さにおけるピーク幅	特定の幅を持つピーク	デフォルトでも、必ずしも通常の意味のFWHMではない
wlen	prominenceやwidthを計算する探索範囲	局所的または大域的に目立つピーク	短すぎると広い構造を正しく評価できない
rel_height	widthを測る高さ	ピーク幅の評価位置	ピークの絶対高さではなく、prominenceを基準にする
plateau_size	平坦な頂上のサンプル数	指定した長さのplateauを持つピーク	ADC量子化だけでなく、飽和やクリッピングでも生じる

時間ビン幅を $\Delta t$ とすると、distance や width を物理時間へ変換するには、

t
=
N_{\rm sample}
\Delta t

を使います。

例えば、

distance = 20

であっても、$\Delta t=0.01~{\rm s}$ なら $0.2~{\rm s}$、$\Delta t=1~{\rm s}$ なら $20~{\rm s}$ です。

コードに設定値を書くときは、

dt = 0.1  # s
minimum_peak_separation = 2.0  # s

distance = int(np.ceil(minimum_peak_separation / dt))

のように、物理時間からサンプル数へ変換すると意味が明確になります。

---

Appendix B：よくある質問（FAQ）

Q1. argrelmax でplateauが検出できないのはなぜ？

argrelmax は、対象点が左右の比較対象よりも厳密に大きいことを要求します。

plateauでは、隣接点と値が等しいため、

x_i
>
x_{i+1}

が成立しません。

したがって、平坦な極大値は検出されません。

単に「差分を使っているから」というより、内部で np.greater による厳密な不等号比較を使っていることが本質です。

---

Q2. find_peaks は、plateau_size を指定しないと平坦ピークを検出できない？

いいえ。

find_peaks は、plateau_size を指定しなくても平坦ピークを局所最大値として検出します。

平坦部の中央のインデックスがピーク位置として返されます。平坦部のサンプル数が偶数の場合は、中央二点のうち小さい側のインデックス、すなわち切り捨て側が返されます。

plateau_size は、

• plateauの長さによってピークを選別する
• plateau_sizes
• left_edges
• right_edges

を取得する

ためのパラメータです。

plateauの情報だけを取得し、plateauの長さでは除外したくない場合は、

peaks, properties = find_peaks(
    data,
    plateau_size=(None, None),
)

と指定できます。

---

Q3. threshold と prominence は何が違う？

threshold は、ピークと直近の左右サンプルとの差を見ます。

そのため、鋭く尖ったピークの検出に敏感ですが、時間ビン幅や高周波ノイズの影響を受けやすい量です。

一方、prominence は、ピークの周囲にある谷を基準として、そのピークが周辺構造からどれだけ目立つかを評価します。

概念的には、

• threshold：すぐ隣と比べて、どれだけ高いか
• prominence：周囲の山谷全体と比べて、どれだけ独立して目立つか

という違いです。

---

Q4. ノイズが多いときは、どのパラメータから調整すればよい？

普遍的に正しい順番はありません。

まず、

• どの時間スケールの現象を検出したいか
• 最低どの程度の振幅が必要か
• ベースラインが一定か
• ピークが鋭いのか広いのか

を考える必要があります。

一般には、

1. 物理的な最小時間間隔から distance
2. 局所的な目立ち方から prominence
3. 必要な絶対強度から height
4. 必要なら width

を検討すると理解しやすいでしょう。

ただし、パラメータを一組だけ採用するのではなく、妥当な範囲で変化させて、結論が維持されるか確認することが重要です。

---

Q5. width はFWHMと同じ？

常に同じではありません。

find_peaks の width は、ピークのprominenceを基準にした高さで測定されます。

デフォルトの

rel_height = 0.5

は、prominenceの半分に相当する高さで幅を測ります。

ベースラインがゼロで、単独のピークが存在する理想的な場合にはFWHMに近くなることがありますが、変動するベースラインや近接ピークがある場合には、通常の意味でのFWHMと一致しません。

---

Q6. 吸収線のような下向きピークを検出するには？

データの符号を反転します。

absorption_peaks, properties = find_peaks(
    -data,
    prominence=prominence_min,
)

ただし、連続光のトレンドがあるスペクトルでは、先に連続光モデルや局所ベースラインを評価する必要があります。

---

Q7. 平均ショットプロファイルがきれいなら、同一形状のショットが存在する？

必ずしもそうではありません。

形状が異なる多数のイベントでも、ピークをそろえて平均すれば滑らかな平均形状が得られます。

また、連続的な赤色雑音から局所最大値を選んだ場合にも、滑らかな山型の平均プロファイルが得られることがあります。

個々のプロファイルの分布、振幅と幅の相関、模擬データとの比較も調べる必要があります。

---

Q8. 平均値と中央値のどちらで重ねるべき？

平均値は、総光子数や平均的な強度変化を議論する場合に自然ですが、極端に明るい少数のイベントに影響されます。

中央値は外れ値に強い一方、光子数の加算という意味は失われ、統計誤差の扱いも平均値とは異なります。

可能であれば、

mean_profile = np.mean(peaks_list, axis=0)
median_profile = np.median(peaks_list, axis=0)

の両方を比較するとよいでしょう。

両者が大きく異なる場合、少数の強いイベントが平均プロファイルを支配している可能性があります。

---

Q9. ショットの誤差は、各時間ビンのPoisson誤差を足せばよい？

光子統計だけを考えるならPoisson誤差は重要ですが、実際にはイベントごとの形状の違いも不確かさに含まれます。

したがって、

• 光子統計による誤差
• ショット間のばらつき
• ピーク時刻決定の不確かさ
• 背景推定の不確かさ
• 検出条件による系統誤差

を区別する必要があります。

平均プロファイルの再現性を調べるには、ピーク単位または観測区間単位のbootstrapが有効です。

---

Q10. ショット解析とパワースペクトル解析は、どちらを使うべき？

どちらか一方を選ぶ必要はありません。

パワースペクトル解析は、変動パワーを周波数ごとに分解し、特徴的な時間スケールを系統的に調べることに向いています。

ショット解析は、特定の条件を満たすピークの前後で、

• 強度
• 硬度
• スペクトル
• 他波長の信号

がどう変化するかを直感的に調べることに向いています。

同じデータに両方を適用し、互いに矛盾しないか確認することが重要です。

---

参考文献・リンク

SciPy

• SciPy find_peaks documentation
  https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.find_peaks.html

• SciPy argrelmax documentation
  https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.argrelmax.html

• SciPy peak_prominences documentation
  https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.peak_prominences.html

• SciPy peak_widths documentation
  https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.peak_widths.html

• SciPy peak-finding source code
  https://github.com/scipy/scipy/tree/main/scipy/signal

ショット解析

• Negoro, H., Miyamoto, S., & Kitamoto, S. 1994,
  “Structure of X-Ray Shots of Cygnus X-1 in Its Low State,”
  ApJL, 423, L127
  NASA ADS

• Yamada, S. et al. 2013,
  “Rapid Spectral Changes of Cygnus X-1 in the Low/Hard State with Suzaku,”
  ApJL, 767, L34
  NASA ADS
  arXiv

• Focke, W. B., Swank, J. H., & Wai, L. 2005,
  “Time Domain Studies of X-Ray Shot Noise in Cygnus X-1,”
  ApJ, 633, 1085
  arXiv

• Bhargava, Y. et al. 2022,
  “Probing the shot behaviour in Cygnus X-1 using simultaneous AstroSat–NICER observation,”
  MNRAS, 512, 6067
  arXiv

---

最後に

本記事では、Pythonを用いて一次元データから複数のピークを検出する方法を、単なる関数の使い方ではなく、

• 離散データにおけるピークの定義
• 差分による検出
• plateauの問題
• argrelmax と find_peaks の違い
• height、distance、prominence、width の意味
• 検出したピークの重ね合わせ

という流れで説明しました。

ピーク検出では、アルゴリズムが「真のピーク」を自動的に発見してくれるわけではありません。

アルゴリズムが行っているのは、解析者が与えた条件に基づいて、局所最大値を選別することです。

したがって、最も重要なのは、

自分が何をピークと定義し、その定義によって、どのようなデータを選び、どのようなデータを捨てているのか

を理解することです。

ショット解析についても同様です。

ピークを重ね合わせて得られた平均形状は、天体の時間変動について有用な情報を与えます。しかし、それは常に、ピーク検出条件のもとで得られた条件付き平均です。

検出条件への依存性、模擬データ、ランダム時刻、別エネルギーバンド、別観測区間との比較を行って初めて、平均プロファイルに含まれるどの特徴が天体固有なのかを議論できるようになります。

find_peaks は便利な関数ですが、重要なのは関数を呼び出すことではありません。

その関数が何をピークとして選び、その選択が最終的な科学結果にどのような影響を与えるかを考えることが、ピーク解析の本質です。