1/17変光星の周期解析方法のレビュー
植村誠 （広島大学）

2/17

動機
• 2013年の連星・変光星・低温度星研究会＠
九州大で「スパースモデリングを使った変光星の周期解析」という発表をした。
• すると「結局のところ一番良い周期解析方法は何なんですか？」という質問を受けることがあった。

• 答えは「最大公約数的な手法はある（＝通常のフーリエ変換）が、『一番良い手法』はデータの性質（データ点数・S/N比・対象天体の情報）に依存する」だと思う。
• その辺りのことをレビューしてみます。
• 前半はフーリエ変換の復習です。

スパースモデリングによる周期解析の例

フーリエ変換とパワースペクトル推定
光度曲線（横軸：時間）3/17

パワースペクトル（横軸：周波数）• 通常は複素表記（直感的にわかりにくいので今回はパス）
フーリエ変換

フーリエ逆変換

Wikipedia フーリエ変換より• sin, cos の畳み込みとして線形変換（連立一次方程式）で書ける（フーリエ級数：次頁）＝光度曲線を様々な周波数の三角関数の重ね合わせで表現する。

4/17
周波数ごとの係数フーリエ変換に相当する行列光度曲線

ある周波数 νi のパワーは𝑃 𝜈𝑖 = 𝑎𝑖2 + 𝑏𝑖2

フーリエ逆変換に相当する行列

光度曲線
周波数ごとの係数• N個の光度曲線データ(y)から、2M個の周波数成分の係数を推定する問題。

5/17

サンプリング定理• ナイキスト周波数：𝑓𝑁𝑦𝑞 = 1/(2∆𝑡) 、つまり、サンプリング周波数の 1/2• N点の光度曲線データから、𝑓𝑁𝑦𝑞 以下の周波数領域のM=N/2点ずつの下式の係数 a,b が決定でき、それら a,b から光度曲線が完全に再構成できる。
• 連立一次方程式の変数の数と方程式の数が同じなら解ける、ということ（フーリエ級数の正規直交性）。
• ただし、一定のサンプリング間隔（ ∆𝑡 が一定）を想定
実際のデータはデータ間隔は非均一。問題（p.7)  p.12以降の手法

6/17

ナイキスト周波数とパワースペクトルの折り返し• 例：• 3つの周期 1.25日、3日、10日で変光している変光星がある。周波数にすると 0.8, 0.333, 0.1 /day。
• この星を毎日ある決まった時刻に観測する（ ∆𝑡 = 1 day)。ナイキスト
周波数は 1/2 = 0.5 /day。
• この時、光度曲線から得られるパワースペクトルには1(= 1/∆𝑡 )−0.8=0.2 /day = 5日の偽信号（エイリアス）が現れる。
観測されるパワースペクトル真のパワースペクトル

P(ν)

P(ν)

0.10.50.333(ナイキスト)0.8

ν0.10.50.333(ナイキスト)0.8ν

• 実際は「5日」の信号が偽物かどうかは判断できない。
• ただし、天体の既知の情報から、ナイキスト周波数以下には信号がなく、ナイキスト周波数以上のどの辺りの周波数に信号がある、ということがわかっていれば、その周波数の推定は可能。

窓関数・エイリアス・ホワイトノイズ

7/17• 現実のサンプリング間隔が非均一な光度曲線
 sin, cos畳み込み(p.4)でパワースペクトル計
算：何をしているのか？

• 等間隔データの欠損部分を計算していない＝明るさ（等級）「０」（ゼロ）と仮定していることと数理的には同じ(0-padding) エイリアス時間欠損がある場合の光度曲線• 窓関数：データのある部分を”1”、ない部分を”0”とした関数。（次頁）• 実際に得られるパワースペクトルは真のパワースペクトルと窓関数のパワースペクトルが畳み込まれたもの。
• 非均一なサンプリング間隔  複雑な窓関数パワースペクトル• 有限のサンプリング期間  幅を持つ信号  推定周波数の不定性（p.14）• サンプリング定理  完全再構成＝データのノイズまで説明する（パワースペ
クトルのホワイトノイズ）• エイリアスやホワイトノイズを含んだパワースペクトルはモデルの予測性能を下げる（過学習）

8/17

9/17

普通のフーリエ変換でパワースペクトル計算する時の注意点まとめ• ナイキスト周波数を意識しよう• 計算してるパワースペクトルの周波数帯域は大丈夫？
• 実際のデータ間隔は非均一でナイキスト周波数が明確に決まらないこともある• サンプリング定理を意識しよう• 光度曲線の点数に対して、周波数のパワーをやたら多く計算してないか？

• エイリアスと窓関数を意識しよう• その信号、偽物かもしれません• この点については p.12以降の手法が有効10/17

Phase Dispersion Minimization (PDM)• Stellingwerf (1978)• 光度曲線をある周期で折り畳んで(folding)、位相ごとの分散を計算。真の周期のところで分散が小さくなる。
• 統計量• 全データ点数N, 位相ビンの数 M, j番目のビンの分散 sj, データ個数 nj••

正しい周期でfoldingした時(𝑛𝑗 −1)𝑠𝑗2(𝑥𝑖 −𝑥)2 2, 𝑠 = 𝑛 −𝑀
𝑁−1𝑗22
Θ = 𝑠 /𝜎 の極小を探す

𝜎2 =誤った周期でfoldingした時• 特長：三角関数の和で表しにくい、非対称性の強い波形に対して有効（矮新星のスーパーハンプ、等）。
Stellingwerf (1978)

11/17食連星やパルス波形の場合• 三角関数の重ね合わせでは難しそうなので、PDM？
• 大きく変光していない位相のデータは周期決定にあまり役に立たない• 大きく変光している領域の極大・極小時刻を計算し（極付近を２次関数で近似する、等）、線形回帰などで周期を推定するのがベター
パルス or 食的な光度曲線

t1t2t3t4t5t6

時間

t6t5t4t3t2t1
123456 このデータと直線との差を縦軸に取ると所謂、O-C diagram12/17Cleanアルゴリズム
窓関数のパワースペクトル• 普通のフーリエ変換でパワースペクトルを作ってから、強い信号に窓関数のパワースペクトルを当てはめて差し引いていく（周期解析への応用はRoberts, Lehar, & Dreher 1987)

Dirty spectrum• 特長：エイリアスに強い• 信号をいくつ検出するかは別の方法で選択する必要がある（Cleanだけでは決まらない）cleaned spectrum• P.8 のような場合は誤った信号を検出してしまう（非常に意地悪な設定）
Hartman+07

Lomb-Scargle periodgram

13/17• ある周波数の三角関数を光度曲線に当てはめ、𝜒 2 を計算。横軸周波数に対して縦軸𝜒 2をプロットしたものを period gram とする• データが等間隔の場合は通常のパワースペクトルに一致する。
• データが不等間隔の場合は窓関数の効果によって、周波数によって出やすい信号・出にくい信号が発生する。それを補正したもの。
• 特長：窓関数によるパワーの「感度」補正ができる  真の信号を見つけやすい• 計算時間がやや掛かる

VanderPlas 2017

14/17LASSOによるパワースペクトル推定（スパースモデリング）• 光度曲線に欠損がある（＝非等間隔データ）とフーリエ逆変換の問題が解けない（解が無限に存在）• パワースペクトルのうち、ゼロ出ない信号強度の周波数はわずかしか存在しないと仮定できる場合• 光度曲線は複数の周期の重ね合わせである• パワースペクトルがスパース• 最少二乗項に加えて制約項として一次ノルムを加える  無数にある解のうち最もスパースな解が選択される• スパース度（＝信号の個数）は交差検証法などで決める• 特長：エイリアス・ホワイトノイズに強い。
• 参考
•
•

「天文屋のための How to スパースモデリング」（植村）Kato & Uemura (2012)

不定性の評価• ブートストラップ法が安全• ある程度のデータ量（できれば >100点）が必要• 乱数を使ってデータを再サンプリングする• 光度曲線データ100点あると、そこから乱数を用いて重複有りで100点サンプリングし、周期を求める。これをたくさん繰り返して、周期の平均・分散を計算する。

• 様々な手法で可能（特にLASSO推定では必須）• フーリエ変換・PDM• データの観測期間で決まる信号の「幅」が不定性• 信号の極大・極小付近を２次関数で近似、HWHMなどで不定性を評
価する• Lomb-Scargle• χ2 と光度一定モデルの χ2 が計算されるので、統計的検定の枠組みで不定性の評価が可能

15/17

（おまけ）
MCMCによる周期推定• マルコフ連鎖モンテカルロ法(MCMC)：モデルパラメータの（事後）確率分布を推定する• R, pythonの汎用ライブラリも充実(ex.
Stan)• データがnoisyで周期性があるか・ないかを知りたい時に有効かも（？）• 周期1つの三角関数を当てはめる。
𝑓 𝑡 = 𝐴 sin 2𝜋𝑡𝜈 + 𝜙 + 𝑓0
• 周期性「ある」「なし」の判断
はモデル選択理論を使う• 交差検証法や情報量規準• 特長１：パラメータの最適値だ
けでなく分布全体や不定性がわかる。
• 特長２：調べる周波数を手で設定する必要がない。
• 強いエイリアスがたくさん出る
場合は不向き

16/17

まとめ17/17• どういう時に何を使えば良いのか？
データが等間隔Yes

普通のフーリエ変換

No

データがnoisyで、受かるとしてもせいぜい信号1つ

Yes

普通のフーリエ変換＋MCMC？ (p.16)
・PyAstronomy/pdm.py
・京都のグループからNo

しっかり見えている信号の周期を知りたいNo

エイリアスが邪魔で、信号がどこに何個あるのか知りたい
Yes

LASSO (p.14)• Kato & Uemura (2012)

Yes

波形が非対称Yes

PDM(p.10)、極大・極小の線形回帰(p.11)

No

普通のフーリエ変換• データが良質で知りたい周期が1つだけなら LASSOやMCMCのような最近の手法を使うまでもない。
• Clean (p.12) はLASSOの「手動版」に近い。
• Lomb-Scargle (p.13) は窓関数の効果を多少補正するが、エイリアスは残る。