1. はじめに：なぜライブラリを選ぶ必要があるのか？
X線天文学における時系列解析では、観測ギャップ（不等間隔データ）やポアソンノイズが避けられません(図1)。そのため、FFT（高速フーリエ変換）の代わりにLomb-Scargle Periodogram が標準的に使われます。
しかし、Pythonにはこの実装が乱立しており、それぞれ計算速度や結果の「値」が異なります。 本記事では、主要な6つのライブラリを全て実装・検証し、その違いを定量的に明らかにしました。
2. 比較した6つのライブラリ
| No. | ライブラリ | 特徴 | 実装アルゴリズム |
|---|---|---|---|
| 1 | Scipy | 標準ライブラリ。最も手軽。 | $O(N^2)$ (Townsend) |
| 2 | Astropy | 天文標準。論文作成に最適。 | $O(N \log N)$ (Press & Rybicki) |
| 3 | Stingray | X線特化。Leahy正規化など実務機能が満載。 | $O(N \log N)$ (Astropyベース) |
| 4 | Gatspy | Astropyの前身。Scikit-learnライクなAPI。 | $O(N \log N)$ (Press & Rybicki) |
| 5 | P4J | 情報理論ベース。外れ値に強い「エントロピー最小化」。 | $O(K N^2)$ (Correntropy) |
| 6 | cuSignal | GPU爆速。CuPyに統合されたGPU版Scipy。 | $O(N^2)$ (Parallel Bruteforce) |
3. 定量評価：ベンチマーク結果
Google Colab (T4 GPU) にて、$N=5,000$ 点の擬似X線データ（周期信号+ポアソンノイズ+外れ値）を用いて検証を行いました。
3-1. 計算速度と鋭さの比較表
各ライブラリで同じデータを処理した際の結果一覧です
| ライブラリ | 平均実行時間 (ms) | 速度評価 | FWHM (Hz) | ピークの鋭さ |
|---|---|---|---|---|
| cuSignal | 4.2 ms | S+ (爆速) | 0.010 | 標準 |
| Astropy | 45.0 ms | A (高速) | 0.010 | 標準 |
| Gatspy | 48.0 ms | A (高速) | 0.010 | 標準 |
| Scipy | 620.0 ms | B (普通) | 0.010 | 標準 |
| Stingray | 750.0 ms | B- (やや重い) | 0.010 | 標準 |
| P4J | 1,200.0 ms | C (重い) | 0.002 | S (極めて鋭い) |
速度とNの比較
3-2. 結果の考察
A. 速度の王者：cuSignal (CuPy)
ScipyやAstropyと比較して 100倍以上高速 です。アルゴリズム自体は工夫のない $O(N^2)$ ですが、GPUの数千コアによる並列計算（Parallel Bruteforce）でねじ伏せています。全天サーベイやシミュレーションには必須です。
B. 鋭さの王者：P4J
表の中で唯一、FWHM（半値全幅）が極端に小さい（鋭い）ことがわかります。他のライブラリが「最小二乗法」に基づいているのに対し、P4Jは「カレントロピー（情報ポテンシャル）」を最大化します。これにより、ノイズや外れ値を自動的に無視してフィッティングするため、埋もれた信号を鋭く切り出せます。
C. 安定のAstropy/ Gatspy
$O(N \log N)$ の高速アルゴリズム（NFFT）を実装しており、CPU処理としては最速クラスです。StingrayはAstropyをラップしている分、わずかにオーバーヘッドがありますが、機能性を考えれば誤差範囲です。
ソフトウェアごとの実行結果
4. なぜ結果が変わるのか？
「同じLomb-Scargleなのに、なぜ結果が違うの？」 ここからは、実務でハマりやすい 3つの「アルゴリズムの罠」 について、Colabでの検証結果を交えて解説します。
1. 浮動平均 (Floating Mean) 問題：0Hzの亡霊
ScipyやGatspy（デフォルト設定）を使用した際、「低周波（0Hz付近）に桁違いの巨大なパワーが検出される」現象です。
【数理的背景】
古典的な Lomb-Scargle（Scipyの scipy.signal.lombscargle 等）は、観測データ $y(t_i)$ に対して以下の正弦波モデルを当てはめ、そのカイ二乗誤差を最小化します。
$$
y(t) = A \cos(\omega t) + B \sin(\omega t)
$$
このモデルは、データの平均値が $\mu = 0$ であることを暗黙に仮定しています。
しかし、実際のX線データ（光子カウント）には、必ず正のオフセット（バックグラウンド $C_{bg}$）が含まれます。
もしオフセットを含んだままこのモデルでフィッティングを行うと、アルゴリズムは定数項 $C_{bg}$ を表現するために、周期無限大（$\omega \to 0$）の成分 を無理やり使おうとします。
結果として、0Hz極限においてスペクトル密度が発散に近い挙動を示します（PSD正規化の場合）。
$$
\lim_{\omega \to 0} P(\omega) \propto C_{bg}^2 \quad (\text{Fixed Mean Model})
$$
【解決策：一般化Lomb-Scargle】
AstropyやStingrayは、Zechmeister & Kürster (2009) で提唱された**「一般化 Lomb-Scargle (Generalized Lomb-Scargle)」**を採用しており、定数項 $C$ を含む以下のモデルを使用します（Floating Mean）。
$$
y(t) = A \cos(\omega t) + B \sin(\omega t) + \mathbf{C}
$$
最小二乗法において $A, B$ と同時に $C$ も最適化パラメータとして扱われるため、データのオフセットは定数項 $C$ に吸収されます。これにより、パワースペクトル上の 0Hz 成分（DC成分）は抑制され、真の信号成分のみが検出されます。
-
教訓: 前処理で Detrend（平均引き）を行うか、Astropy (
fit_mean=True) を使用してください。
② 正規化 (Normalization) の哲学：値の大きさの意味
Stingrayだけパワーの値が数千に達し、他は1以下になる現象です。これは「誰に向けて作られたツールか」という設計思想の違いです。
【Astropy / Scipy：数学者の視点 (Standard)】
確率論的に扱いやすい Standard正規化 を採用しています。データの分散（Reference $\chi^2$）に対し、周期モデルによってどれだけ残差が減ったかを $0 \sim 1$ で評価します。
$$
P_{std}(\omega) = \frac{\chi^2_{ref} - \chi^2(\omega)}{\chi^2_{ref}}
$$
これは相関係数の二乗（$R^2$）に近い指標であり、単位系に依存しない汎用性を持ちます。
【Stingray：X線天文学者の視点 (Leahy)】
X線解析のデファクトである Leahy正規化 (Leahy et al. 1983) を採用しています。総光子数を $N_{ph}$、離散フーリエ変換を $DFT(\omega)$ とすると、以下のように定義されます。
$$
P_{Leahy}(\omega) = \frac{2}{N_{ph}} \left| \sum_{j=1}^{N} x_j e^{-i \omega t_j} \right|^2
$$
この正規化の最大の利点は、統計的性質が明確なことです。
もしデータが純粋なポアソンノイズ（白色雑音）のみで構成されている場合、そのパワースペクトル分布は**自由度2のカイ二乗分布（$\chi^2_2$）**に従います。
$$
\text{Prob}(P_{Leahy} > z) = e^{-z/2}
$$
つまり、パワーの値 $P$ を見るだけで、「$P=30$ ならば、これがノイズである確率は $e^{-15} \approx 3 \times 10^{-7}$ である」といった有意性検定が即座に可能になります。
③ グリッド解像度 (Binning)：スカロッピング損失
「鋭いピークがあるはずなのに、パワーが低く出る」現象です。
【原因：サンプリング定理の死角】
周波数探索を行う際の「独立な周波数グリッド」の間隔 $\Delta f$ は、観測時間を $T_{obs}$ として通常以下のように定義されます。
$$
\Delta f = \frac{1}{T_{obs}}
$$
Stingrayなどのデフォルト設定は、計算コスト削減のためにこのグリッドを使用します。
しかし、真の信号周波数 $f_{true}$ がグリッドの中間点に位置してしまった場合、検出されるパワーは最大で以下のように減衰します（Scalloping Loss）。
$$
\text{Loss} \approx 1 - \text{sinc}^2(0.5) \approx 36%
$$
【対策：オーバーサンプリング】
これを防ぐには、理論上の解像度よりも細かく周波数をサンプリングする必要があります。Astropyの autopower メソッドの nyquist_factor を上げることで、グリッド間隔を $\Delta f / n$ に細分化し、ピークの頂点を正確に捉えることができます。
- 教訓: パルサーのようなコヒーレントな（Q値が高い）信号を探すときは、デフォルト設定を信じず、必ず**オーバーサンプリング（5〜10倍程度）**を行ってください。
結論： 天文データ解析では、必ず Floating Mean 対応の Astropy を使いましょう。