最小二乗スペクトル解析

Periodicity computation method

データポイントの集合を二次関数で近似した結果

最小二乗スペクトル解析LSSA )は、フーリエ解析に似た、データサンプルへの正弦波の最小二乗近似に基づいて周波数スペクトルを推定する方法です[1] [2]科学で最もよく使用されるスペクトル手法であるフーリエ解析は、一般に、長くギャップのある記録で長周期ノイズを増幅しますが、LSSAはそのような問題を軽減します。[3]フーリエ解析とは異なり、LSSAを使用するためにデータが等間隔である必要はありません。

LSSAは1969年[4]と1971年[5]に開発され、ペトル・ヴァニチェクにちなんでヴァニチェク法ガウス・ヴァニチェク法とも呼ばれ、[6] [7]また、ニコラス・R・ロンブ[9]とジェフリー・D・スカーグル[10]による簡略化に基づいてロンブ法[3 ]またはロンブ・スカーグル・ピリオドグラム[2] [8]としても知られています

歴史的背景

フーリエ解析ピリオドグラム、そして正弦波の最小二乗フィッティングの間には密接な関係があることは古くから知られていました。[11] しかし、ほとんどの開発は等間隔のサンプルの完全なデータセットに限定されていました。1963年、アムステルダムのMathematisch CentrumのFreek JM Barningは、同様の手法を用いて不等間隔のデータを扱いました。 [12]これには、今日ではロンブ法と呼ばれるピリオドグラム解析と、そのようなピリオドグラムから決定された正弦波の選択された周波数の最小二乗フィッティングの両方が含まれていました。そして、これらは今日ではポストバックフィッティングによるマッチング・パーシュート[13]または直交マッチング・パーシュート[14]として知られる手順によって結び付けられました。

ニューブランズウィック大学カナダ人地球物理学者測地学者、ペトル・ヴァニーチェクは、1969年に等間隔および不等間隔のデータに対するマッチング・パースート法を提唱しました。彼はこれを「逐次スペクトル分析」と呼び、その結果を「最小二乗ピリオドグラム」と名付けました。[4]彼はこの手法を一般化し、単純な平均値を超える体系的な要素、例えば「予測される線形(二次、指数関数など)の、大きさが不明な永年的傾向」などを考慮できるようにしました。そして1971年に、様々なサンプルに適用しました。[5]

ヴァニチェクの厳密な最小二乗法は、1976年にシドニー大学のニコラス・R・ロンブによって簡略化され、ピリオドグラム解析との密接な関係が指摘されました[9] その後、不等間隔データのピリオドグラムの定義はNASAエイムズ研究センターのジェフリー・D・スカーグルによって修正・解析され、[10]わずかな変更を加えることで、個々の正弦波周波数をフィッティングするためのロンブの最小二乗式と同一になることが示されました。

スカーグルは、この論文は「新しい検出技術を導入するものではなく、観測時間が不均等な場合の、最も一般的に使用される技術であるピリオドグラムを用いた検出の信頼性と効率性を研究するものである」と述べ、さらに、ピリオドグラム解析と比較した正弦波の最小二乗近似に関して、この論文は「(提案された修正により)これら2つの方法が全く同等であることを、明らかに初めて証明した」と指摘している。[10]

プレス[3]はこの発展を次のように要約している。

これらの困難を軽減し、他の非常に望ましい特性も備えた、不均一にサンプリングされたデータに対するスペクトル分析のまったく異なる方法が、ロンブによって開発されました。この方法は、バーニングとヴァニセックによる以前の研究に部分的に基づいており、さらにスカーグルによって詳細化されました。

1989年、オンタリオ州キングストンのクイーンズ大学のマイケル・J・コーレンバーグは、スペクトルやその他の問題のほぼ最適な分解をより迅速に見つける「高速直交探索」法を開発しました。[15]これは後に直交マッチング追求として知られるようになった手法に似ています。

LSSAとその変種の開発

ヴァニーチェク法

線形回帰では、観測値()は、従属変数(y)と独立変数(x )の間の根底にある関係(青)からのランダム偏差()の結果であると仮定されます。そして、最小二乗法などのノルムフィッティングでは、データポイント(赤)はノルム的に最もよく適合する直線()によって表され、そこから常に「残差」()が残ります。

ヴァニチェク法では、離散データセットは、標準的な線形回帰または最小二乗近似を用いて、段階的に決定された周波数の正弦波の加重和によって近似されます。[16] 周波数はバーニング法に類似した方法で選択されますが、最小二乗近似後の残差を最小化する周波数を選択することで、各連続する新しい周波数の選択をさらに最適化します(現在、事前バックフィッティングを伴うマッチング追求法として知られているフィッティング手法に相当します[13])。正弦波の数は、データサンプルの数以下でなければなりません(同じ周波数の正弦波と余弦波は別々の正弦波としてカウントします)。

データベクトルΦは、重みベクトルxを使用してサンプル時間で各関数を評価することによって行列Aに表された正弦波基底関数の加重和として表されます

ϕ A x , {\displaystyle \phi \approx {\textbf {A}}x,}

ここで重みベクトルxは、 Φを近似する際の二乗誤差の合計を最小にするように選択される。xの解は標準的な線形回帰を用いた閉形式である[17]

x = ( A T A ) 1 A T ϕ . {\displaystyle x=({\textbf {A}}^{\mathrm {T} }{\textbf {A}})^{-1}{\textbf {A}}^{\mathrm {T} }\phi .}

ここで、行列Aは、サンプル時間において評価される際に互いに独立な(必ずしも直交する必要はない)任意の関数の集合に基づくことができる。スペクトル解析に用いられる関数は、典型的には、対象となる周波数範囲にわたって均等に分布する正弦関数と余弦関数である。あまりにも狭い周波数範囲であまりにも多くの周波数を選択すると、関数の独立性が不十分となり、行列の条件が悪くなり、結果として得られるスペクトルは意味をなさなくなる。[17]

Aの基底関数が直交(つまり、相関がなく、列同士のドット積がゼロ)の場合、行列A T Aは対角行列です。すべての列のべき乗(要素の平方和)が同じ場合、その行列は単位行列に定数を掛けたものになるため、逆変換は簡単です。後者は、サンプル時間が等間隔で、正弦波として、サンプルあたり 0 から半周期までの周波数間隔で等間隔​​にペアになった正弦と余弦を選択した場合(サンプルあたり 1/N 周期間隔で、0 と最大周波数で正弦波位相が共にゼロになるものは省略)に当てはまります。この場合は離散フーリエ変換と呼ばれ、測定値と係数の観点から少し書き直されています。[17]

x = A T ϕ {\displaystyle x={\textbf {A}}^{\mathrm {T} }\phi } —スカラー係数内で、N 個の等間隔のサンプルと周波数の DFT ケース。

ロンブ法

ピリオドグラムで計算された2つの正弦波基底関数のパワースペクトル(振幅の2乗)

1976年にVaníček法の計算負荷を軽減しようとして[9](現在は問題ではない)、Lombは、同じ周波数の正弦基数と余弦基数の間の相関はペアワイズで小さいことが多いため、上記の簡略化を一般に用いることを提案した。これは、少なくとも間隔が狭くない場合は、正弦波のペア間の相関が小さいためである。この定式化は、基本的には従来のピリオドグラムの定式化であるが、間隔が不均等なサンプルで使用できるように改良されている。ベクトルxは基礎スペクトルのかなり良い推定値であるが、相関を無視するため、A xはもはや信号の良い近似値ではなく、この方法はもはや最小二乗法ではないが、文献では依然としてそのように呼ばれている。

スカーグルは、正弦波と余弦波の波形を直接データにドット積するだけでなく、標準的なピリオドグラムの公式を修正して、まず時間遅延を求めるようにした。これにより、この2つの正弦波はサンプル時間で互いに直交し、さらにこれら2つの基底関数の潜在的に不等なべき乗も調整され、周波​​数における電力のより正確な推定値が得られる。[3] [10]この手順により、スカーグルの修正ピリオドグラム法はロンブの方法と完全に同等になった。時間遅延は定義により次の式に等しい 。 τ {\displaystyle \tau } t j {\displaystyle t_{j}} τ {\displaystyle \tau }

tan 2 ω τ = j sin 2 ω t j j cos 2 ω t j . {\displaystyle \tan {2\omega \tau }={\frac {\sum _{j}\sin 2\omega t_{j}}{\sum _{j}\cos 2\omega t_{j}}}.}

すると、周波数におけるピリオドグラムは次のように推定されます。 ω {\displaystyle \omega }

P x ( ω ) = 1 2 [ [ j X j cos ω ( t j τ ) ] 2 j cos 2 ω ( t j τ ) + [ j X j sin ω ( t j τ ) ] 2 j sin 2 ω ( t j τ ) ] , {\displaystyle P_{x}(\omega )={\frac {1}{2}}\left[{\frac {\left[\sum _{j}X_{j}\cos \omega (t_{j}-\tau )\right]^{2}}{\sum _{j}\cos ^{2}\omega (t_{j}-\tau )}}+{\frac {\left[\sum _{j}X_{j}\sin \omega (t_{j}-\tau )\right]^{2}}{\sum _{j}\sin ^{2}\omega (t_{j}-\tau )}}\right],}

スカーグルの報告によれば、これは均等にサンプリングされたケースのピリオドグラムと同じ統計分布を示す。[10]

任意の周波数において、この方法はその周波数の正弦波に最小二乗近似を適用した場合と同じパワーを与える。[18] ω {\displaystyle \omega }

ϕ ( t ) = A sin ω t + B cos ω t . {\displaystyle \phi (t)=A\sin \omega t+B\cos \omega t.}

実際には、与えられたロンブピークが有意であるかどうかを判断することは常に困難であり、特にノイズの性質が不明な場合は困難です。例えば、ノイズの多い周期信号のロンブピリオドグラム解析で誤ったスペクトルピークが検出された場合、乱流データのノイズが原因である可能性があります。[19]フーリエ法では、パッチを当てたデータや編集されたデータを解析する際に、誤ったスペクトルピークが報告されることもあります。[7]

一般化ロンブ・スカーグル・ピリオドグラム

標準的なロンブ・スカーグル・ピリオドグラムは、平均値がゼロのモデルに対してのみ有効です。一般的には、ピリオドグラムを計算する前にデータの平均を差し引くことで近似値が求められます。しかし、モデル(フィッティングされた正弦波)の平均がゼロでない場合、この仮定は不正確です。一般化ロンブ・スカーグル・ピリオドグラムでは、この仮定がなくなり、平均値を明示的に求めます。この場合、フィッティングされる関数は[20]です。

ϕ ( t ) = A sin ω t + B cos ω t + C . {\displaystyle \phi (t)=A\sin \omega t+B\cos \omega t+C.}

一般化ロム・スカーグル・ピリオドグラムは、文献では浮動平均ピリオドグラムとも呼ばれています。[21]

コーレンベルグの「高速直交探索」法

オンタリオ州キングストンクイーンズ大学のマイケル・コーレンバーグは、スペクトル解析の正弦波成分などの過剰完備集合から疎な成分集合を選択する手法、高速直交探索(FOS)を開発した。数学的には、FOSは平均二乗誤差低減(MSER)プロセスにおいて、疎行列反転として実装されたわずかに修正されたコレスキー分解を使用する。[15] [22] 他のLSSA手法と同様に、FOSは離散フーリエ解析の大きな欠点を回避しているため、埋め込まれた周期性を正確に識別し、不等間隔データで優れた性能を発揮する。高速直交探索法は、非線形システム同定などの他の問題にも適用された。

パーマーのカイ2乗法

パーマーは、任意の数の高調波に対して最適な関数を見つける手法を開発し、非正弦高調波関数をより自由に見つけられるようになりました。[23]彼は、非一様標準誤差を持つ任意間隔のデータに対して、重み付き最小二乗法による 高速(FFTベース)解析手法を開発しました。この手法を実装したソースコードは入手可能です。[24] データは均一間隔の離散時間でサンプリングされないことが多いため、この手法では、時系列配列をサンプル時間で疎に埋め込むことでデータを「グリッド化」します。すべてのグリッドポイントには統計的重みがゼロに設定され、これはサンプル間の時間に無限のエラーバーを持つことと同等です。

アプリケーション

パラメータの異なる値に対するベータ分布

LSSA の最も便利な機能は、データを操作したり、存在しないデータを作り出したりする ことなく、不完全な記録をスペクトル分析できることです。

LSSAスペクトルの振幅は、時系列分散に対する周波数または周期の寄与を表します[4]一般的に、このように定義されたスペクトル振幅は、出力の有意水準レジームを容易にします。[25]あるいは、ヴァニーチェクスペクトルのスペクトル振幅はdB で表すこともできます[26]ヴァニーチェクスペクトルのスペクトル振幅はβ分布に従うことに注意してください[27]

ヴァニチェクのLSSAの逆変換は可能であり、順方向変換を行列として記述することで最も簡単に確認できます。逆行列(行列が特異でない場合)または擬似逆行列は逆変換になります。選択された正弦波がサンプル点で相互に独立しており、その数がデータ点の数と等しい場合、逆変換は元のデータと完全に一致します。[17] このような逆手順は、ピリオドグラム法では知られていません。

実装

LSSAはMATLABコード1ページ未満で実装できます[28]要するに:[16]

「最小二乗スペクトルを計算するには、m個のスペクトル値を計算する必要があります。これは、異なる周波数の[スペクトルパワー]を取得するために、毎回最小二乗近似をm回実行することを含みます。」

すなわち、所望の周波数セット内の各周波数について、データサンプルに対応する時刻において正弦関数余弦関数を評価し、データベクトルと正弦波ベクトルのドット積をとって適切に正規化する。ロンブ/スカーグル・ピリオドグラムとして知られる手法に従い、ドット積の前に正弦成分と余弦成分を直交化するために、各周波数の時間シフトを計算する。[17]最後に、これら2つの振幅成分からべき乗を計算する。データが時間的に均一に配置され、選択された周波数が有限データレコード上の整数サイクル数に対応する場合、 この同じプロセスは離散フーリエ変換を実現する。

この手法は、各正弦波成分がデータ点と直交していない場合でも、それらを独立して、あるいは文脈を無視して扱う。これはヴァニーチェクの独自の手法である。さらに、行列方程式を解き、指定された正弦波周波数間のデータ分散全体を分割することで、完全な同時最小二乗法、あるいは文脈を考慮した最小二乗法を実行することも可能である。[17]このような行列最小二乗法は、MATLABにおいてバックスラッシュ演算子としてネイティブに利用可能である[29]

さらに、同時法またはインコンテキスト法は、独立法またはアウトオブコンテキスト法(ロンブによるピリオドグラム法も同様)とは対照的に、データサンプルの数よりも多くの成分(正弦と余弦)を当てはめることができないため、次のようになります。[17]

「...選択された周波数によってフーリエ成分(三角関数)の一部が互いにほぼ線形依存するようになり、悪条件またはほぼ特異なNが生成された場合、深刻な影響が生じる可能性があります。このような悪条件を回避するには、推定する周波数の異なるセット(たとえば、等間隔の周波数)を選択するか、N(つまり、非対角ブロック)の相関を無視して、個々の周波数について逆最小二乗変換を個別に推定する必要があります...」

一方、ロンブのピリオドグラム法では、標準的なピリオドグラムと同様に、任意の数、あるいは密度の周波数成分を使用することができる。つまり、周波数領域を任意の係数でオーバーサンプリングすることができる。[3]しかし、上述のように、ロンブの単純化と最小二乗基準からの逸脱は、彼の手法に重大な誤差をもたらし、誤ったスペクトルピークをもたらすことにも留意する必要がある。[19]

フーリエ変換離散フーリエ変換などのフーリエ解析では、データにフィッティングされる正弦波はすべて相互に直交しているため、単純なコンテキスト外ドット積ベースの基底関数への投影とコンテキスト内の同時最小二乗フィッティングの間に区別はありません。つまり、異なる周波数の直交正弦波間の分散を最小二乗分割するために逆行列は必要ありません。[30]過去には、等間隔のサンプルを含む完全なデータレコードが利用できる場合、フーリエの高速フーリエ変換実装の処理効率が高いため、多くの人が選択した方法であり、ギャップのあるレコードの解析にもフーリエファミリーの手法が使用されていましたが、そのためには、存在しないデータを操作したり、作り上げたりする必要がありました。

参照

参考文献

  1. ^ Cafer Ibanoglu (2000). 変光星を天体物理学の必須ツールとして. Springer. ISBN 0-7923-6084-2
  2. ^ ab D. スコット・バーニー、デイヴィッド・オスパー、ギレルモ・ゴンザレス (2006). 観測天文学. ケンブリッジ大学出版局. ISBN 0-521-85370-2
  3. ^ abcde Press (2007). Numerical Recipes (第3版). Cambridge University Press. ISBN 978-0-521-88068-8
  4. ^ abc P. Vaníček (1969年8月1日). 「最小二乗法による近似スペクトル解析」(PDF) .天体物理学と宇宙科学. 4 (4): 387– 391. Bibcode :1969Ap&SS...4..387V. doi :10.1007/BF00651344. OCLC  5654872875. S2CID  124921449.
  5. ^ ab P. Vaníček (1971年7月1日). 「最小二乗法によるスペクトル解析のさらなる発展と特性」(PDF) .天体物理学と宇宙科学. 12 (1): 10– 33. Bibcode :1971Ap&SS..12...10V. doi :10.1007/BF00656134. S2CID  109404359.
  6. ^ J. Taylor; S. Hamilton (1972年3月20日). 「Vaníček法によるスペクトル解析のいくつかの検証」.天体物理学と宇宙科学. 17 (2): 357– 367. Bibcode :1972Ap&SS..17..357T. doi :10.1007/BF00642907. S2CID  123569059.
  7. ^ ab M. Omerbashich (2006年6月26日). 「Sepkoski大要のGauss-Vanicekスペクトル解析:新たなライフサイクルはなし」. Computing in Science & Engineering . 8 (4): 26– 30. arXiv : math-ph/0608014 . Bibcode :2006CSE.....8d..26O. doi :10.1109/MCSE.2006.68.
  8. ^ Hans PA Van Dongen (1999). 「生物学的リズムの探索:不等間隔データのピリオドグラムにおけるピーク検出」. Journal of Biological Rhythms . 14 (6): 617– 620. doi :10.1177/074873099129000984. PMID  10643760. S2CID  14886901.
  9. ^ abc Lomb, NR (1976). 「不等間隔データの最小二乗頻度解析」.天体物理学と宇宙科学. 39 (2): 447– 462. Bibcode :1976Ap&SS..39..447L. doi :10.1007/BF00648343. S2CID  2671466.
  10. ^ abcde Scargle, JD (1982). 「天文学的時系列解析の研究 II - 不均一間隔データのスペクトル解析の統計的側面」.アストロフィジカル・ジャーナル. 263 : 835. Bibcode :1982ApJ...263..835S. doi :10.1086/160554.
  11. ^ デイヴィッド・ブラント(1931年)『観察の組み合わせ』(第2版)ケンブリッジ大学出版局。
  12. ^ バーニング、FJM (1963)。 「12 のクサザルの光曲線の数値解析」。オランダ天文研究所の紀要17 : 22。ビブコード:1963BAN....17...22B。
  13. ^ ab Pascal Vincent; Yoshua Bengio (2002). 「カーネルマッチング追跡」(PDF) .機械学習. 48 ( 1–3 ): 165–187 . doi : 10.1023/A:1013955821559 .
  14. ^ YC Pati、R. Rezaiifar、およびP.S. Krishnaprasad、「直交マッチング追求:ウェーブレット分解への応用を伴う再帰関数近似」、 Proc. 27th Asilomar Conference on Signals, Systems and Computers、 A. Singh編、Los Alamitos, CA, USA、IEEE Computer Society Press、1993年
  15. ^ ab Korenberg, MJ (1989). 「システム同定と時系列解析のための堅牢な直交アルゴリズム」.生物サイバネティクス. 60 (4): 267– 276. doi :10.1007/BF00204124. PMID  2706281. S2CID  11712196.
  16. ^ ab Wells, DE, P. Vaníček, S. Pagiatakis, 1985.「最小二乗スペクトル解析の再考」測量工学部技術報告書84、ニューブランズウィック大学、フレデリクトン、68ページ、[1]で入手可能。
  17. ^ abcdefg Craymer, MR、「最小二乗スペクトル、その逆変換および自己相関関数:理論と測地学への応用」、博士論文、トロント大学、カナダ (1998)。
  18. ^ ウィリアム・J・エメリー、リチャード・E・トムソン(2001年)『物理海洋学におけるデータ分析法』エルゼビア社、ISBN 0-444-50756-6
  19. ^ ab Zhou, W.-X.; Sornette, D. (2001年10月). 「重裾相関ノイズにおける周期性と対数周期性の統計的意義」International Journal of Modern Physics C . 13 (2): 137– 169. arXiv : cond-mat/0110445 . Bibcode :2002IJMPC..13..137Z. doi :10.1142/S0129183102003024. S2CID  8256563.
  20. ^ M. Zechmeister; M. Kürster (2009年3月). 「一般化ロンブ・スカーグル・ピリオドグラム.浮動平均法とケプラー式ピリオドグラムのための新しい形式論」. Astronomy & Astrophysics . 496 (2): 577– 584. arXiv : 0901.2573 . Bibcode :2009A&A...496..577Z. doi :10.1051/0004-6361:200811296. S2CID  10408194.
  21. ^ Andrew Cumming、Geoffrey W. Marcy、R. Paul Butler (1999年12月). 「リック惑星探査:検出可能性と質量閾値」. The Astrophysical Journal . 526 (2): 890– 915. arXiv : astro-ph/9906466 . Bibcode :1999ApJ...526..890C. doi :10.1086/308020. S2CID  12560512.
  22. ^ Korenberg, Michael J.; Brenan, Colin JH; Hunter, Ian W. (1997). 「高速直交探索によるラマンスペクトル推定」. The Analyst . 122 (9): 879– 882. Bibcode :1997Ana...122..879K. doi :10.1039/a700902j.
  23. ^ Palmer, David M. (2009). 「不規則サンプルデータの周期探索のための高速カイ2乗法」. The Astrophysical Journal . 695 (1): 496– 502. arXiv : 0901.1913 . Bibcode :2009ApJ...695..496P. doi :10.1088/0004-637X/69​​5/1/496. S2CID  5991300.
  24. ^ 「David Palmer: 高速カイ二乗周期探索」。
  25. ^ Beard, AG, Williams, PJS, Mitchell, NJ & Muller, HG 惑星波と潮汐変動の特殊気候学、J Atm. Solar-Ter. Phys. 63 (09), p.801–811 (2001)。
  26. ^ Pagiatakis, S. 最小二乗スペクトルにおけるピークの確率的意義、J of Geodesy 73、p.67-78 (1999)。
  27. ^ Steeves, RR「最小二乗スペクトルにおけるピークの有意性の統計的検定」測地測量論文集、エネルギー・鉱山・資源省、測量・地図作成、オタワ、カナダ、p.149-166 (1981)
  28. ^ リチャード・A・ミュラー、ゴードン・J・マクドナルド (2000). 『氷河期と天文学的原因:データ、スペクトル解析、メカニズム』(第1版). シュプリンガー・ベルリン・ハイデルベルク. Bibcode :2000iaac.book.....M. ISBN 978-3-540-43779-6OL  20645181M。ウィキデータ Q111312009。
  29. ^ ティモシー・A・デイビス;カーミット・シグモン (2005)。 MATLAB入門書。 CRCプレス。ISBN 1-58488-523-8
  30. ^ ダレル・ウィリアムソン (1999). 離散時間信号処理:代数的アプローチ. シュプリンガー. ISBN 1-85233-161-5
Retrieved from "https://en.wikipedia.org/w/index.php?title=Least-squares_spectral_analysis&oldid=1318304067#The_Lomb.E2.80.93Scargle_periodogram"