統計数理（2023）
第 71 巻 第 1 号 25–45©2023 統計数理研究所

特集「諸科学における統計数理モデリングの拡がり」［原著論文］

射影演算子法による統計的時系列データ解析とその応用
前山 伸也1 ・三分一 史和2（受付 2022 年 6 月 13 日；改訂 9 月 13 日；採択 9 月 20 日）

要旨

射影演算子法は非平衡統計物理学の分野で開発された数理手法であり，興味ある目的変数の時間発展を，説明変数に対する相関項と無相関項に分離し，一般化 Langevin 形式で記述する．この手法は，データセットに対する実用的な統計的時系列データ解析として利用できる．
本研究では，その定式化を再訪しつつ，連続時間および離散時間システムでの任意の目的・説
明変数に対する拡張を行った．続いて，類似した時系列データ解析手法として構造ベクトル自モデルとの比較を行った．射影演算子法は任意の説明変数に対する定式化であ己回帰（SVAR）り SVAR モデルよりも広い適用範囲を持つこと，説明変数として自己回帰形を用いた場合はSVAR モデルと同等の相関を抽出可能であることを示した．さらに応用例として，プラズマ乱流による帯状流生成・維持過程を一般化 Langevin 方程式として捉える解釈を提案し，射影演算子法による時系列データ解析に基づく根拠を与えた．以上により，物理的解釈を与えるためのデータ解析，および，確率的時系列データのモデリングに対する射影演算子法の有用性を示した．開発した手法は Python からのシンプルな関数呼び出しで利用できるオープンソースコー
ドとして公開し，興味を持った読者は手軽に試すことができる．
キーワード：射影演算子法，連続/離散時間システム，乱流，時系列解析．
1.

はじめに

1.1 本研究の背景
森 - Zwanzig の射影演算子法は，非平衡統計物理の分野で開発された数理手法である．原著にあたる Zwanzig（1961）や Mori（1965）の論文では，主として多粒子系における遅い集団運動とミクロな揺動力を分離する目的で用いられた．ある種の粗視化手法であるとともに，決定論的な運動方程式から確率的な一般化 Langevin 方程式を得る方法論としても知られる．量子系の密度行列の時間発展を記述するマスター方程式，ガラス転移におけるモード結合理論（川崎,，分子動力学における散逸粒子動力学モデルの導出（Español, 1995）といった，多様な分2009）
野で射影演算子法の応用がなされてきた．射影演算子法から得られる記憶関数を解析的に導出することは一般に困難であり，実用的な表式を得るために，しばしば，摂動的に高次微小量を無視する，高次相関を低次相関の積として近似する，粒子間力を二体間力の和で表わしつつ時間スケールの分離を仮定するなど，それぞれの現象に応じたアドホックな仮定に依存した理論
1

名古屋大学大学院 理学研究科：〒 464–8602 愛知県名古屋市千種区不老町2

統計数理研究所：〒 190–8562 東京都立川市緑町 10–3

26

統計数理

第 71 巻 第 1 号 2023

構築が行われてきた．射影演算子法について，非平衡統計力学に焦点を当てて解説している文，藤坂（1998），宗像（1996）を挙げておく．
献として 鈴木（2000）
これらとはやや異なる視点から，Unresolved 変数の効果を取り入れた Resolved 変数の時間発展の最適予測という文脈で射影演算子法の定式化を再訪したのが Chorin et al.（2000, 2002）
である．以降，この論文を起点として，記憶関数のモデル化や Navier-Stokes 乱流への適用など，多数の論文が報告されている（Stinis, 2004; Chorin and Stinis, 2006; Hald and Stinis,．さらに最近では，
2007; Stinis, 2013, 2015; Parish and Duraisamy, 2017; Gouasmi et al., 2017）
分子動力学や乱流現象を対象とした数値シミュレーションデータ解析への応用（Meyer et al.,2020; Maeyama and Watanabe, 2020）
や，射影演算子法と外生変数を含む非線形自己回帰移動
平均モデル（NARMAX, nonlinear auto-regression moving average with exogenous inputs）の比，動的モード分解などの Koopman 学習との関係性
較（Chorin and Lu, 2015; Lin and Lu, 2021）が議論されるなど（Lin et al., 2021），統計的時系列データ解析としての射影演算子法の理論的理解が進展してきている．
連続時間システムか離散時間システムか，（ii）射影射影演算子法はその定式化によって，（i）
目的変数と説明変数の関係，（iv）単変量解析か多変量解析か，といった演算子の選び方，（iii）
分類ができるが，これらの具体的な区別は第 2 節の定式化の説明に譲る．
実用的な統計データ解析として射影演算子法を用いるために，記憶関数を近似なしにデータから直接評価している例は少ない．そこで，上述の文献を数値的評価手法によって再分類してで線形射影のアルゴリズムが提案され，Zwanzig（1980）
ではみよう．原理的には，Mori（1965）
高次項を加えることで非線形（有限ランク）
射影が構築できることが提案された．Chorin et al.
（2002）がおそらく時系列データ解析としての最初の実装例と思われる．時系列データから相互相関関数を計算し，記憶関数に関する式を台形公式で離散的に積分することで記憶関数を数値的に評価した．ただし，評価された記憶関数に関する図は明示されなかった．日本でも，森本人とその同僚により記憶関数の数値的評価がされたが（Okamura, 2006; Mori and Okamura,「記憶関数に関する式を反復法的に解いた」とのことで，それ以上2007），実装方法についてはの詳細な記述はなかった．単変量解析の場合はよく研究されている 1 つの未知変数についての第 2 種 Volterra 積分方程式に帰着するため，説明不要と判断したのかもしれない．近年では，
Meyer et al.（2020）
は反復法による記憶関数の評価手法を提案しており，同時期に，Maeyamaは任意関数を目的変数とする場合の定式化や，Laguerre 展開により記憶and Watanabe（2020）
項を Markov 埋め込みすることで時間積分する手法を提案している．続く Lin et al.（2021）では多変量解析や離散時間射影演算子法に対する実装がなされ，さらに最近の Lin et al.（2022）の研究では，森の線形射影のアルゴリズムではなく，回帰を用いた非線形射影への拡張がなされている．
1.2 本研究の目的と論文構成
本論文では，個別の物理現象に依らない，汎用的な統計的時系列データ解析手法としての射影演算子法の利用に焦点を当てる．この観点からは，1 変数や線形射影については多くの計算・応用例が報告されているが，多変量・非線形・任意関数に注目した解析は少なく，これらを複合的に扱った解析は十分に議論がなされていないというのが現状である．こうした現状に対して，本研究では以下の 3 つの問題

（1）任意目的変数に対する射影演算子法の定式化が十分に整理されていないこと（2）時系列データ解析に用いられる自己回帰モデルとの類似点が見られるが，どのような相違点があるか明らかでないこと

射影演算子法による統計的時系列データ解析とその応用

27（3）射影演算子法に基づく時系列データ解析をどのように物理解析に適用すれば意味ある結果を得られるかは問題次第だが，応用事例がまだ多くはないことを解決することを目的とした．
（1）に関連して，Lin et al.（2021）らの以降の論文構成は次の通りである．第 2 節では，問題定式化を拡張することで，任意関数に対する連続時間・離散時間射影演算子法を定式化し，それぞれに成り立つ一般化第 2 種揺動散逸定理を導く．また，これらが実用的な統計的時系列
データ解析アルゴリズムとして利用できることを説明する．第 3 節では，射影演算子法の利用（2）に関連して，構造ベクトル自己回帰モデルと射影演算子法に例を示す．第 3.1 節では，問題（3）に関連して，実物理現象への応用事例として長谷よる解析を比較する．第 3.2 節では，問題
川 - 若谷プラズマ乱流シミュレーション時系列データの解析・モデリングへの適用を行う．最後に，第 4 節で本研究を総括する．
2.

射影演算子法の定式化

ここで考える統計的時系列データ解析の問題設定や変数の記号を定義しよう．統計的時系列
データ解析では，図 1 のような多数の変数の時系列データセットから，変数間の情報を抽出す
t
n）
と個々のサ
ることが目的である．目的変数のデータは f (l) (t) は連続時間 （または離散時間
ンプル l の次元を持つ．説明変数も同様に時間 t とサンプル l の次元を持ち，多変量のベクト
ル g (l) (t) = (g0(l) (t), g1(l) (t), . . .) である．私たちは，統計平均操作を用いて，f (l) (t) がどのように
g (l) (t) と相関するかを評価したい．同様の解析を複数の目的変数に繰り返し行う場合は，目的変数もベクトル f (l) (t) となる．統計平均操作 · · ·  は，データセットのサンプル l に関する平均（サンプル平均）により算出される．例えば，異なる時刻 t + s と t における相互時間相関は以下で与えられる．
1  (l)∗(l)fi (t + s)gj (t),
S
S−1

(2.1)

fi (t + s)gj (t) =l=0

ここで，S はサンプル数である．別の統計平均操作の定義としては，エルゴード性を仮定し長時間平均を利用する場合もあり得るが本論文では用いない．以下では，サンプル番号に関する上付き添え字 (l) は，曖昧さが生じない場合は表記の簡単のため省略する．
加えて，本論文中で用いる「相関」という語について先に定義を明記しておく．射影演算子法「相関」とは，単に相互時間相により抽出される，目的変数 f (t) が説明変数 g(t) に対して持つ）の意味ではなく，一般化 Langevin 形式（式（2.19），（2.24））に現れる Markov 係数 Ω
関（式（2.1）

図 1．目的変数と説明変数の時系列データセット．

28

統計数理

第 71 巻 第 1 号 2023と記憶関数 Γ を指す意味で用いる．目的変数から相関項（Markov 項と記憶項）を除いた残差を「無相関」
であるとは，任意の時刻 t に
無相関項 r(t) と呼ぶが，この項が説明変数 g(t) に対して
おいて無相関項 r(t) と説明変数の初期値 g(0) との相互時間相関（式（2.1））が 0 となることを意味する．これらの数式上の表現は第 2.1 節，第 2.2 節で顕わに導く．
2.1 連続時間射影演算子法
2.1.1 森 - Zwanzig 形式
現象を記述する物理量 u(t) = (u0 (t), u1 (t), . . .) は以下のような連立非線形常微分方程式系に従うとする．
(2.2)du(t)= N (u(t)),
dt

ここで，t は連続時間，N は任意の非線形作用素である．粒子多体系のような離散系のみならず，流体現象のような連続系においても，Galerkin 法のように基底関数を用いて展開する場合
や差分法のように離散格子点で表す場合を考えれば，上記は一般的な問題を記述できる．
目的変数や説明変数は，物理量の何らかの関数 f (t) = F (u(t)), g(t) = G(u(t)) としてあらわされる．実験観測的な例でいえば，計測できる気温 (=: f ) と風速 (=: g) の関連を調べる際に，これらは気体運動論の立場からは粒子分布関数 (=: u) のモーメント量として記述されるが，必ずしも分布関数そのものは計測できない（する必要がない）という状況に対応する．あるいは，高解像度の大規模数値シミュレーションの全情報 (=: u) を調べずとも，大きなスケールの変動(=: f , g) 間の関係のみ抽出したいなどという状況も当てはまるだろう．
以下の解析では，系全体が統計的定常状態にあると仮定しよう．つまり，任意関数 f (t), g(t)に対して，統計平均 · · ·  を用いて定義される時間相関関数は相対時間のみに依存し，(2.3)

f (t + s)g(t) = f (s)g(0),という時間並進対称性が成り立つ．このため，初期時刻 t = 0 は特別な意味を持たず，初期値
f (l) (0), g (l) (0) は単に定常確率分布から拾い上げられたあるサンプルの値に過ぎない．相対時
間のみが意味を持つという意味で，f (s)g ≡ f (t + s)g(t) や f g ≡ f (t)g(t) の表記を適宜用いる．
に倣って，以下の 3 つの演算子を定義す
射影演算子法の議論を進めるに当たり，岡村（2014）
る．一つ目は射影演算子 P である．射影演算子は冪等率 PP = P を満たす演算子であり，この段階ではまだこれ以上は特定しない．二つ目は補射影演算子 Q = 1 − P であり，P の冪等性から QQ = Q, PQ = QP = 0 が確認できる．三つめは Liouville 演算子 Λ であり（山本・中村,
，
1999）

du ∂
∂
(2.4)
Λ=
Ni (u)
,
·
=
dt ∂u
∂u
i
i
で定義される．これは物理量 u の任意の関数に作用し，例えば fi (t) = Fi (u(t)) について，
(2.5)
(2.6)

dfi (t)= Λfi (t),dtfi (t) = etΛ fi (0),という性質を満たす．ここで，etΛ は冪級数展開で定義される演算子指数関数であり，fi (0) から fi (t) まで時間を t だけ進める時間進展演算子とみなすことができる．
森 - Zwanzig 形式では，任意関数の時間発展式（2.6）を以下のように射影と補射影に分離する．

射影演算子法による統計的時系列データ解析とその応用

(2.7)

29

fi (t) = etΛ fi (0),= etΛ Pfi (0) + etΛ Qfi (0), te(t−s)Λ PΛesQΛ Qfi (0)ds + etQΛ Qfi (0).
= etΛ Pfi (0) +0

ここで，式（2.7）
第 1 段では初期値を射影と補射影に分離し，続く第 2 段では Dyson 公式
ttΛ(t−s)Λ
e = 0 dse
PΛesQΛ + etQΛ を適用している．これにより，初期値の補射影部分 Qfi (0) がと補射影に入り続ける部分（右辺第 3 項）を分離
時間発展とともに射影に入る部分（右辺第 2 項）している．特に後者を(2.8)

ri (t) ≡ etQΛ Qfi (0),と定義する．このとき，ri (0) = Qfi (0), dri (t)/dt = QΛri (t) であり，任意の時刻 t について
(2.9)

Pri (t) = 0,

として常に射影から除かれることが確認できる．このため，ri (t) は無相関項または直交項(orthogonal dynamics) と呼ばれる．あるいは，乱数的な挙動を期待してノイズ項や揺動項などと呼ばれることもあるが，あくまで ri (t) の重要な性質は，定義より常に補射影に入り続ける
こと（式（2.9））であり，必ずしも乱数的に見なせるとは限らない点に注意が必要である（岡村,．ri (t) の定義式（2.8）
を用いて，射影演算子による分解を行った式（2.7）を書き換えると
2014）
 t
(2.10)fi (t) = etΛ Pfi (0) +e(t−s)Λ PΛri (s)ds + ri (t),0を得る．式（2.10）を一般化 Langevin 形式と呼び，右辺第 1 項は Markov 項，第 2 項は過去の履歴に依存する非 Markov 型の項で記憶項と呼ばれる．特別な場合として，後述する森の線形射
影を射影 P として採用し，さらに説明変数の時間微分項を目的変数として用いる場合，よく知られた一般化 Langevin 方程式の形に帰着する．
2.1.2 射影演算子 P の選択
ここまでは射影と補射影の分解と記憶の獲得という一般的な形式論であるが，実用的な計算を進めるには，射影演算子 P を具体的に定義する必要がある．冪等率 PP = P を満たす任意の演算子が適用可能であるため幅広い選択肢が考えられるが，ここでは，目的変数 fi をある
説明変数 g に射影するクラスの演算子を考えることにする．Chorin et al.（2002）は，Zwanzig（1980）の非線形射影を再訪し，説明変数 g による目的変数 fi の最適推定を与える射影演算子
として，与えられた g の下での fi の条件付き期待値関数 E を射影演算子として採用した．
(2.11)
PE fi = E[fi |g] = fi Φfi |g (fi , g)dfi ,
ここで，Φfi ,g (fi , g) は fi , g の同時確率密度，Φg (g) = Φfi ,g (fi , g)dfi は g の周辺確率密度，
Φfi |g (fi , g) = Φfi ,g (fi , g)/Φg (g) は g が与えられた下での fi の条件付確率密度である．現象を記述する全自由度の物理量 u の確率密度 Φ(u) を用いて定義される場合もあるが，上記のように観測量 fi , g に関する確率密度のみに縮約できるため，fi , g の観測データが十分にあれば条件付き期待値関数を数値的に評価することが可能である．この関係を図 2 に模式的に示した．
射影として条件付期待値を用いる場合，PE fi は g を引数として fi の値を返す非線形関数として与えられる．
同様の観点から，Mori（1965）の射影演算子法を再訪すると，森の線形射影は上記の条件付
き期待値の線形近似として与えられる．

統計数理

30

第 71 巻 第 1 号 2023

図 2．統計的データセットと条件付き期待値関数，射影演算子 P の関係．多数のサンプル（a）散布図，および，
データ l について説明変数 g (l) と目的変数 f (l) を評価することで（b）f , g の同時確率密度分布 Φf,g を得る．これらから算出した条件付き期待値 E[f |g]を図（a）中に実線で表す．非線形射影 PE は条件付き期待値関数で表される非線形曲線，
線形射影 PL はその線形近似，有限ランク射影 PF や回帰ベース射影 PR は非線形基底
または回帰による曲線近似に対応する．

PL fi = fi g ∗  · gg ∗ −1 · g =(2.12)

j

ここで，fi gj∗  =∗ −1fi gj∗ (gg ∗ −1 )jk gk ,kfi gj∗ Φfi ,g (fi , g)dfi dg を 2 変数相関と呼ぶ．上付き添え字 ∗ は複素共役を表し，gg  は逆行列である．式（2.12）は目的変数 fi の説明変数 g に対する多変量線形回帰での範囲で最適推定を与える．
あり，線形近似（図 2 の条件付き期待値関数の直線近似）
上述の条件付き期待値関数による射影 PE は理想的な最適推定であり，非線形関係も表せるが，実現象において限られた数の観測データから確率密度分布を構築することは難しい．一方，森の線形射影 PL は，線形関係であるため数値実装も簡便であり，得られた結果も解釈しやすいため広く用いられてきた．これら両者の中間に位置する射影演算子をさらに 2 つ紹介する．
．
一つは，有限ランク射影（finite-rank projection）と呼ばれる方法である（Chorin et al., 2002）
説明変数 g は u の任意の関数であったので，説明変数として g 2 , g 3 などの高次項を加えて森の線形射影を適用することで，有限の多項式基底による条件付き期待値の非線形近似を構築することができる．
(2.13)

PF fi = fi h∗  · hh∗ −1 · h,

ここで，h(g) は任意の g の高次項を含むという意味で，条件付き期待値の非線形近似（図 2 の曲線近似）となっているが，アルゴリズムとしては森の線形射影（2.12）と全く同等に算出できる．
もう一つは，回帰ベース射影 PR と呼ばれる方法である（Lin et al., 2022）．多項式，スプラ
イン，あるいはニューラルネットワークなどの種々の回帰を射影演算子として用いることで，有限のデータセットから fi と g の非線形関係（図 2 に対する回帰）を表現する．射影演算子法
への適用としては比較的新しい提案手法であるが，文献 Lin et al.（2022）では，回帰ベース射影と森の線形射影や有限ランク射影の比較を行った結果，回帰ベース射影を用いた方がいくつかのケースで時系列データの再現性が良いこと，回帰ベース射影のアルゴリズムにおいて線形回帰を用いた場合は森の線形射影と一致することなどが報告されている．
2.1.3 森の線形射影 PL に対する一般化 Langevin 形式と第 2 種揺動散逸定理
本論文では射影演算子として説明変数の初期値 g(0) に対する森の線形射影 P = PL（式（2.12））を採用する．繰り返しになるが，統計的定常性を仮定しているため初期時刻 t = 0 というのは単に便宜的な呼び方で，その起点からの相対時間のみが意味をもつ．なお，説明変数 g射影演算子法による統計的時系列データ解析とその応用

31に高次項が含まれるとして再定義すれば，有限ランク射影 PF も内包される．以下では表記の簡略化のため下付き添え字 (L) は省略する．このとき，無相関項 ri の性質（式（2.9））は，(2.14)

ri (t)g(0) = 0,となり，ri (t) は任意の時刻 t で説明変数の初期値 g(0) と相関を持たないことが分かる1) ．
森の線形射影（式（2.12））を用いて，一般化 Langevin 形式（2.10）の各項を書き下してみよう．
まず右辺第 1 項の Markov 項は(2.15)etΛ Pfi (0) =Ωij gj (t),j

(2.16)∗

Ω = f g  · gg ∗ −1 ,ここで，Ω は Markov 係数行列と呼ばれ，統計的定常性（2.3）により時間依存性を持たず同時刻
相関 f g ∗  = f (0)g ∗ (0) = f (t)g ∗ (t) のみで決まる定数行列であるため，時間進展演算子 etΛの影響を受けない．次に，式（2.10）
右辺第 2 項の記憶項は
 t t(t−s)Λ
(2.17)
ePΛri (s)ds = −Γij (s)gj (t − s)ds,00
dg ∗ (0)(2.18)
Γ(s) = − Λr(s)g ∗ (0) · gg ∗ −1 = r(s)· gg ∗ −1 ,dt
ここで，Γ(s) は記憶関数（行列）である．途中，Λr(s)g ∗ (0) + r(s)Λg ∗ (0) = 0 を用いた．式（2.18）は記憶関数と無相関項 ri (t) の間に成り立つ関係を示しており，一般化第 2 種揺動散逸定理と呼ばれる．特に，力学的時間発展を記述する文脈ではしばしば目的変数として説明変数の時間微分項 f = dg/dt が採用され，この場合 Γ(s) = r(s)r(0) · gg ∗ −1 というよく知られた形に帰着する．ただし，この定理が物理学における Langevin 方程式と同様にマクロな摩擦とミ
クロな揺らぎの関係を表わす第 2 種揺動散逸定理としての物理的意味を持つためには，r(t) がランダムな揺らぎに関する運動であることが保証される必要がある（岡村, 2014）．最終的に，森の線形射影を適用した場合の一般化 Langevin 形式（2.10）は
 t(2.19)
f (t) = Ω · g(t) −Γ(s) · g(t − s)ds + r(t),0となる．任意の説明変数 g(t) と目的変数 f (t) に対する定式化であり，Markov 係数 Ω や記憶
関数 Γ(s) は 2.3 節で後述するように，g(t) と f (t) の時系列データセットから直接評価できる．
Markov 項は同時刻の説明変数のみから目的変数を f (t) ∼ Ω · g(t) として推定する条件付き期
待値の線形近似であり，記憶項は g(t − s) の過去の履歴をもとに f (t) の推定値を改良する項と（2.14）の意味で，説明変数とは無相関に変動する部分を表わす．
言える．無相関項 r(t) は，式
2.2 離散時間射影演算子法
前節では連続時間システムに対する射影演算子法の定式化を行った．一方，社会現象などの応用例では離散時間システムも存在し，また実験計測や数値シミュレーションで得られるデー
タも離散時間点であることから，離散時間システムに対する射影演算子法も開発されてきた（Munakata, 1994; Darve et al., 2009）．導出の要は連続時間システム同様，時間発展とともに補射影に入り続ける部分を抽出することである．離散時間システムでは，現象を記述する物理量は以下の連立非線形差分方程式系に従うとする．
(2.20)u(n + 1) = N (u(n)),

32

統計数理

第 71 巻 第 1 号 2023

ここで，離散時間 n は整数であり，N は任意の非線形作用素である．離散時間システムでは，Liouville 演算子のような連続的な時間進展演算子は定義できない．代わりに，物理量 u の任意の関数 F (u) に作用し，新たな合成関数を生成する Koopman 演算子 M を以下で定義する（Arnold and Avez, 1968; アーノルド・アベズ, 1972）．
MF ≡ F ◦ N .

(2.21)

定 義 よ り ，任 意 の 関 数 f (n) = F (u(n)) に 対 し て Koopman 演 算 子 を 作 用 さ せ る と ，Mf (n) = MF (u(n)) = F (N (u(n))) = F (u(n + 1)) = f (n + 1) として，1 ステップだけ時間進展させる演算子と見なすことができる．
（- 2.10））と同様に，目的変数 fi (n) の時間発展を射影と補射連続時間システムの場合（式（2.7）
影に分離することで，以下の一般化 Langevin 形式を得る．
(2.22)fi (n) = Mn fi (0),= Mn Pfi (0) + Mn Qfi (0),= Mn Pfi (0) +nMn−m PMri (m − 1) + ri (n),m=1
ここで，離散時間 Dyson 公式 Mn = nm=1 Mn−m PM(QM)m−1 + (QM)n を用いた．また，無相関項を ri (n) = (QM)n Qfi (0) として定義し，任意の時間ステップ n について，
Pri (n) = 0,(2.23)

として常に射影から除かれる部分を抽出した．
を採用する場合，一般化 Langevin 形式（2.22）は以下の射影演算子として森の線形射影（2.12）形式に書き下される．
(2.24)

f (n) = Ω · g(n) −nΓ(m) · g(n − m) + r(n),m=1

ここで，Markov 係数および記憶関数は(2.25)

Ω = f g ∗  · gg ∗ −1 ,(2.26)

Γ(m) = − Mr(m − 1)g(0) · gg ∗ −1 = r(m)g(1) · gg ∗ (1)−1 ,である2) ．無相関項の性質を表わす式（2.23）
は，式（2.14）と同様に，任意の時間ステップ n における無相関項 ri (n) が説明変数の初期値 g(0) に対して相関を持たないことを意味する．式（2.26）右辺は，離散時間射影演算子法において，任意の目的変数に対して拡張された一般化
第 2 種揺動散逸定理であり，筆者の知る限り本論文で新規に導かれた結果である．先行研究
Munakata（1994）; Lin et al.（2022）では，式（2.26）中辺を指して揺動散逸定理と呼ばれた．この場合，第 2.3 節で説明する通り記憶関数の構築自体は可能であるが，背後の全物理量に依存する Koopman 演算子 M を直接評価できないために，得られた記憶関数が揺動散逸定理を満g(t) のみで表されるため，たすかどうかという数値的検証はできない．一方，本導出では r(t)，射影演算子法により記憶関数や無相関項をデータセットから構築したのちに，記憶関数が式（2.26）右辺を満たすかという数値的検証に用いることができる．特別な場合として，目的変数として説明変数の次ステップの値 f (n) = g(n + 1) を用いる場合，Lin et al.（2021）で示された
離散時間システムにおける通常の一般化第 2 種揺動散逸定理 Γ(m) = r(m)r(0) · gg ∗ (1)−1に帰着する．

射影演算子法による統計的時系列データ解析とその応用

33

2.3 実用的時系列データ解析アルゴリズムとしての射影演算子法
当初述べた通り，私たちの目的は，図 1 のような複数サンプルから成る統計的時系列デー
タセットを元に，射影演算子法を適用して変数間の相関を抽出することである．具体的には，Markov 係数行列 Ω と記憶関数行列 Γ(t) を算出することで，目的変数 f (t) を説明変数 g(t) に対して相関を持つ Markov 項と記憶項，相関を持たない無相関項に分離して評価する．ここで，一般化 Langevin 形式を導く際に統計的定常状態を仮定しているため，適用に際しては，解析に用いるデータについて定常性の検定を行うことが推奨される．
Markov 係数 Ω については，データセットに対しサンプル平均をとって同時刻相関行列 f g
と gg を作成し，その逆行列 gg−1 を算出すれば，直ちに定義式（2.16）または（2.25）の通りに評価できる．
や（2.26）には未評価の無相関項 r(t) を含むためそのまま記憶関数については，定義式（2.18）
では評価できない．連続時間システムの場合，一般化 Langevin 形式（2.19）に dg(0)/dt · gg−1をかけて統計平均をとることで，記憶関数 Γ(t) に関する第 2 種 Volterra 積分方程式を得る．
 t(2.27)
Γ(t) = A(t) +Γ(s) · B(t − s)ds,0dg(2.28)
A(t) = [f (t) − Ω · g(t)]· gg−1 ,
dt


dg
B(t − s) = g(t − s)
(2.29)· gg−1 ,dt

ここで，行列 A(t) と B(t) はデータセットに対するサンプル平均から直ちに評価できる量では，時間積分をデータセットの時間刻み幅 Δt について 2 次精度台形公式を用ある．式（2.27）
 Δt
いて 0 Γ(s) · B(t − s)ds = [Γ(0) · B(Δt) + Γ(Δt)B(0)]Δt/2 + O(Δt2 ) と近似することによって，初期値解 Γ(0) = A(0) から逐次的に解くことができる（Maeyama and Watanabe, 2020）．
同様に，離散時間システムに対しては，n − 1 時間ステップ目に対する一般化 Langevin 形式（2.24）に Koopman 演算子 M を作用させ，g(0) · gg−1 をかけて統計平均をとることで，以下（Volterra summation equation）を得る．
の Volterra 和分方程式(2.30)

Γ(n) = C(n) +n−1Γ(m) · D(n − m),m=1(2.31)

C(n) = − [f (n) − Ω · g(n)]g · gg−1 ,(2.32)

D(n − m) = − g(n − m)g · gg−1 ,こちらも，データセットから評価される行列 C, D をもとに，n ≥ 1 に対して Γ(1) = C(1) から逐次的に計算できる．なお，もともと一般化 Langevin 形式は fi (t) の各成分 i について独立で（和分）
方程式についても必ずしも記憶関数 Γ(t) の行列全体をあったので，上記の Volterra 積分
一度に解く必要はなく，各行 i に対する列ベクトル Γij (t)(j = 0, 1, . . .) について独立に解くことができ，計算メモリ容量に対する制約を軽減できる．最後に，無相関項 r(t) については，すでに Markov 係数や記憶関数が評価できたので，各サンプル (l) に対して目的変数 f (l) (t) から
Markov 項と記憶項を除いた残差として求まる．
定式化の部分で触れたように，Markov 係数 Ω や記憶関数 Γ(t) は統計平均量であり，目的変数 f (t) が説明変数 g(t) に対してどのように相関を持つかをデータセットから抽出したものである．例えば，第 3.2 節では，乱流の非線形相互作用項のデータセットに対して射影演算子法
を適用し，評価された記憶関数を元に，記憶項が流れの駆動または減衰に働くか，記憶時間はどのように決定されるかなどの物理的性質を議論した．このことは，射影演算子法により

34

統計数理

第 71 巻 第 1 号 2023

Markov 係数や記憶関数といった時系列データセットの持つ相関を抽出することで，物理的解釈を得るためのデータ解析手法として利用できることを示す例である．
Maeyama and Watanabe（2020）の論文では単変量射影演算子法（f (t) も g(t) も一変数）（f (t) も g(t)の 実 装 を 行 っ た が ，第 3.1 節 の 解 析 を 行 う に 当 た っ て 多 変 量 射 影 演 算 子 法も多変数）の実装を行い，その解析プログラムをオープンソースコードとして公開した（https://github.com/smaeyama/mzprojection/）．これらは Python の関数として提供され，をインポートすることで利用できる．連続時間射影演算子法の場
モジュール（mzprojection.py）
合，時間刻み幅 delta_t，目的変数 f ，説明変数と初期時刻における時間微分 g, dgdt0 の時系列
データセットを用意した後に，
Omega, Gamma = mzprojection_multivariate(delta_t, g, dgdt0, f)の関数呼び出しで，Markov 係数 Omega，記憶関数 Gamma を取得する．離散時間射影演算子法の場合も同様に，解析したい目的変数と説明変数の時系列データセットに対し，
Omega, Gamma = mzprojection_multivariate_discrete_time(g, f)とすればよい．シンプルな関数呼び出しで射影演算子法の評価が行えるので，興味を持った読者はぜひ自身の持つデータセットに試してみてほしい．
2.4 時系列データモデリングへの応用
Markov 係数と記憶関数が目的・説明変数間の相関を抽出したものである一方で，無相関項
r(t) は，目的変数 f (t) の内，変数の初期値 g(0) に対して相関を持たずに時間変動する部分を抽出したものである．一般には背後のダイナミクス u(t) から決定される量であり，単純な乱数過程で近似できるとは限らない．しかし，もしこれを何らかのモデルで与えられるとするなら
ば，目的変数 f (t) を説明変数 g(t) に依存する Markov 項・記憶項と無相関項のモデル r(t) のみで閉じた形で表すことができる．例えば，第 3.1 節や第 3.2 節では，無相関項を有限の振幅・時定数を持つ乱数過程でモデル化した場合に，元のデータセットをどの程度再現できるかについて議論した．このように，射影演算子法は統計的時系列データの相関の解析だけではなく，時系列データの統計的性質の再現を目的とした一般化 Langevin モデリングへも応用可能である．
（2.19）
または（2.24）に基づいて再現性シミュレーションを行う場合，離一般化 Langevin 形式
散時間系では記憶関数の次数 n が大きくない場合には，説明変数の過去のデータをストレージ
し，記憶項を直接評価することができる．一方，連続時間系では，シミュレーションの時間刻
み幅が細かい場合，ストレージしておく過去の説明変数のデータ量が膨大になり，計算機のメ
モリ制約から計算困難に陥る．論文 Maeyama and Watanabe（2020）では，連続時間射影演算子法における記憶項の計算アルゴリズムが提案された．記憶関数が，

∞
p

Cp t
(2.33)
Γ(t) =
e−t/τ ,p! τ
p=0と し て ，係 数 行 列 (Cp )ij お よ び 時 定 数 τ を 用 い て 展 開 で き る 場 合 ，記 憶 項 M (t) ≡tt− 0 Γ(s) · g(t − s)ds = − 0 Γ(t − s) · g(s)ds は，(2.34)

M (t) =∞


Cp · M p (t),p=0(2.35)1
M p (t) = −
p!

 t
0

t−s
τ

p

e−(t−s)/τ g(s)ds,

射影演算子法による統計的時系列データ解析とその応用(2.36)⎧M 0 (t)⎪⎨−g(t) −,dM p (t)τ=⎪dt⎩ M p−1 (t) − M p (t) ,τ35(for p = 0)(for p ≥ 1)として表わされる．つまり，記憶項を計算するために説明変数 g(t) の過去の履歴をすべてストレージしておく必要はなく，近似に必要な次数 p だけの変数 M p に対して Markov 埋め込みさ（2.36）に基づいて時間発展させることで，式（2.34）から記れる．初期値 M p (t = 0) = 0 から式憶項を計算できる．
3.

解析・モデリングへの応用例

この節では，射影演算子法の解析・モデリングへの適用例を示す．第 3.1 節では，第 1.2 節に関して，射影演算子法と構造ベクトル自己回帰モデル（SVAR, Structural
で挙げた問題（2）との比較を行い，その類似点・相違点を明らかにする．続く第Vector Auto-Regressive model）
3.2 節では，問題（3）に関して，物理現象解析の例として長谷川 - 若谷プラズマ乱流への適用によって得られた成果を述べる．
3.1 離散時間射影演算子法と構造ベクトル自己回帰モデルの比較
3.1.1 構造ベクトル自己回帰モデルの数値シミュレーション
射影演算子法により得られる一般化 Langevin 形式（2.19）または式（2.24）は，目的変数が説明
変数の履歴と無相関項により与えられるという観点から，外生変数型自己回帰モデル（ARX，Auto-Regressive model with Exogenous input）との類似性がある．自己回帰型のモデルにおい（Akaike, 1973）によるモデル選択やては赤池情報量規準（AIC，Akaike Information Criterion）
Granger 因果性検定（Granger, 1969）などがこれまでによく研究されており，多変量自己回帰を導入するとインパルス応答分析や周波数ドメイ
モデル（VAR，Vecrot Auto-Regressive model）
ンでのパワー寄与率などの解析を行うことが可能である．時系列データの計測において，計測
対象のシステムの変動周期が時間的に緩慢であったり，計測条件の問題でシステムの変動周期に対して十分に高いサンプリング周波数を設定できない場合には変量間の瞬時的な作用も考慮
モデルが用いられ（Kilian and Lütkepohl, 2017; 沖本, 2010; 村尾, 2019），する構造 VAR（SVAR）
経済学や脳科学など数多くの分野で応用されている．
本節では，類似性を持つ自己回帰型モデルとの比較を通じて，射影演算子法の理解を深めることを目的とした．そのために，良く研究されている VAR モデルの中でも，より適用範囲の広い SVAR モデルとの比較を行った．SVAR モデルは以下の式で与えられる．
(3.1)

Y (n) =p


B(m)Y (n − m) + μ(n),μ(n) ∼ W.N.(Σ),

m=0

ここで，離散時間 n は整数であり，μ(n) は予測誤差で分散共分散行列 Σ を持つ白色正規乱数（I は恒等行列）と置き，式（3.1）の両辺に V −1 を掛けると，である．V = I − B(0)(3.2)

Y (n) =pV −1 B(m)Y (n − m) + V −1 μ(n),

m=1
p

=



A(m)Y (n − m) + ξ(n),m=1となり，係数に各変量の同時点間の相関を含まない誘導形が得られる．ξ(n) の分散共分散行
列は Cov = V −1 Σ(V −1 )T (上付き T は転置行列) であり，一般に非対角行列となるため誤差項

36

統計数理

第 71 巻 第 1 号 2023

図 3．設定した変量間の因果性（実線：過去からの影響，破線:同時点における影響）．
表 1．SVAR シミュレーションパラメータ．

ξ(n) は相関を持つことになる．B(0) はデータから推定した誘導形の Cov にコレスキー分解を施すことにより求められるが，識別性の問題により上あるいは下三角成分（対角成分を除く）
しか推定することができないことが知られている．つまり変量間の瞬時的作用についてはその方向性が定まらないことになるので，解析結果の解釈においては先験的な知見を導入する必要がある．シミュレーションでは 5 つの振動子 Y (n) = (Y0 (n), . . . , Y4 (n)) を仮定した．各々の振動子は異なる減衰率で 0.1Hz で振動し図 3 のような因果性での結合を設定した．シミュレーションデータの生成においてはサンプリング周波数を 1.0Hz とし，表 1 で示すような係数
行列 B(0), B(1), B(2) を用いた．μ(n) = (μ0 (n), . . . , μ4 (n)) は分散 μ2E = 0.1 の独立正規乱数とした．シミュレーションから得られた Y (n) の定常確率分布に従う時系列データから連続 10
時間ステップのデータを 199,980 サンプル切り出すことにより，時系列データセット Y (l) (n)(n = 0, 1, . . . , 9; l = 0, 1, . . . , 199979) を作成した．
3.1.2 射影演算子法による時系列データ解析と一般化 Langevin モデリング
こうして生成されたシミュレーションデータセットに対し，説明変数として g(n) := Y (n)，目的変数として次ステップの値 f (n) := Y (n + 1) を選び，離散時間射影演算子法を適用して，

射影演算子法による統計的時系列データ解析とその応用

37

図 4．SVAR モデルと射影演算子法による解析結果の比較．（a）（b）（c）SVAR モデルシミュ
レーションにおいて表 1 で与えた係数行列 A(n) および誤差項 ξ(n) 分散共分散行列，（e）（f）データセットに離散時間射影演算子法を適用することで得られた Markov 係（d）
数 Ω と記憶関数 Γ(n) および無相関項 r(n) の分散共分散行列．

Markov 係数 Ω や記憶関数 Γ(n) を評価した．SVAR モデルの誘導形（3.2）と離散時間射影演算子法の式（2.24）を見比べると，A(1) ↔ Ω，A(2) ↔ −Γ(1) の類似性が見て取れる．これらの係数の比較を図 4 に示す．射影演算子法により評価された Markov 係数 Ω と記憶関数 Γ(n) がシミュ
レーション上与えた係数行列をおおよそ捉えていることが確認できる．同様に，ξ(n) ↔ r(n)として，誘導形のノイズ項と射影演算子法により評価される無相関項の間にも類似性が見て取れる．図 4 から，射影演算子法により抽出された無相関項が，シミュレーション上与えたノイズ項の相互相関を捉えていることが確認できる．この例は，射影演算子法を統計的時系列デー
タセットに適用することで，SVAR モデルと同等の解析が可能であることを示している．
さらに，データセットから評価された Markov 係数 Ω と記憶関数 Γ(n)，無相関項の相互相
関 rr の値を用いて，目的変数を f (n) = y(n + 1)，説明変数を g(n) = y(n) とした一般化
Langevin モデル（式（2.19））に基づいて時間発展を計算し，元データに対する再現性テストを行った．この際，無相関項は相互相関を満たす正規乱数として与えた．図 5 は元データと再現性テストにより得られた変数を比較したもので，変数の確率密度分布および相互時間相関といった統計的性質をよく再現している．この例は，射影演算子法が統計モデリングにも応用できることを示している．
3.1.3 射影演算子法と自己回帰分析の相違点
以上の解析では両者の類似性が際立つ結果となったが，定式化やデータ解析の実装上は以下のような相違点があることを明記しておく．
• 定式化の違い．自己回帰モデルに基づく時系列データ分析において，自己回帰モデル式のパラメータは予測誤差平方和を最小にする最小二乗法や予測誤差分散の期待値を最小にす
る Yule-Walker 法などを用いて推定される．一方，射影演算子法は，説明変数に対する射

38

統計数理

第 71 巻 第 1 号 2023

図 5．SVAR モデルと射影演算子法による再現性テストの比較．SVAR モデルシミュレー（PDF，probability
ションの元データから得られた（a）変数 Y0 (n) の確率密度分布
distribution function）と（b）変数間の相互時間相関 Y0 (n)Yj ．射影演算子法により（a）変数 y0 (n) の PDF と評価した係数を用いた一般化 Langevin モデルから得られた（b）変数間の相互時間相関 y0 (n)yj ．

影と補射影の分離に基づいた目的変数の関数展開であり，無相関項は必ずしも最小化されない．
• 適用範囲の違い．SVAR モデルは式（3.1）
で与えられ，例えば予測誤差自身が有限の時間相関を持つ場合などは未知入力項を考慮したモデルへの拡張が必要となり，また，説明変数と目的変数が異なる場合なども適切なモデル選択が必要である．一方，射影演算子法はこ
れらのような場合でも広く適用できる．
• 記憶項の評価長さの選び方の違い．自己回帰分析においては，モデル方程式の時間遅れの次数 p は解析者が設定する量であり，次数によって評価されるモデル係数も変わりう
る．適切な次数は，データの再現性とモデル自由度の最小化の観点から，AIC などの基準（2.30）にによって決定される．一方，射影演算子法における記憶関数は式（2.27）または式
従って構成されるため，その上限時刻の取り方に依らず一意に定まる．記憶項の持つ特徴
時間は，（選択した射影演算子 P で表現できる範囲で）データ間の相関が失われるまでの時間であり，背後の物理的ダイナミクスを反映したものであると考えられる．
• 第 2 種揺動散逸定理の有無．自己回帰分析においては予測誤差 μ(t) はモデル係数 B(m) と無関係であり，時間相関を持たない白色乱数になりえる．一方，射影演算子法では，第 2
種揺動散逸定理により無相関項は記憶関数と関連付けられるため，記憶時間が有限であれば無相関項も同等の時間相関を持つ．
こうした違いがありつつも，本節での比較結果から，射影演算子法の目的・説明変数を自己回帰形に選んで解析に適用した場合は，SVAR における推定パラメータと同等の Markov 係数・記憶関数といった相関が抽出可能であることが示された．
3.2 応用例 2. 二次元長谷川 - 若谷プラズマ乱流
3.2.1 本節の物理的背景と射影演算子法を適用する動機
より現実的な物理現象解析の事例として，プラズマ乱流現象への適用について述べる．非平衡状態にあるプラズマ中では，密度・温度勾配といった熱力学的力を駆動源とする微視的不安定性により，複雑な流れを持つ乱流が生じる．プラズマ乱流では，流れが高周波・短波長で乱雑に振舞うように見えても，同時に，より時間・空間的に秩序だった構造を自発的に形成する射影演算子法による統計的時系列データ解析とその応用39ことがある．例えば，磁化プラズマ中で乱流が作りだす帯状流（ゾーナルフロー，Zonal flows）
は，核融合プラズマの閉じ込め改善に寄与するとして広く研究されている（藤澤, 2005）．帯状
流は閉じ込め磁場に沿って一様な構造によって特徴づけられる対称性の高い流れであり，大気準地衡流における Rossby 波が作り出すジェット気流とも共通する自己組織化現象の一例として知られる．
乱流による帯状流の形成機構に関する既存理論として，微視的不安定性により成長する波の振幅が一定以上になると，Kelvin-Helmholtz 不安定性などの 2 次的不安定性により崩壊し，その際に帯状流を励起するといった機構が準線形解析により示されている（Rogers et al., 2000）．
しかし，このような解析はあくまで過渡的な振る舞いであり，準定常的な飽和した乱流状態を記述するものではない．別のよく用いられる方法はエネルギーバランス解析であり，帯状流は線形駆動項を持たないため，時間平均的には，乱流との非線形相互作用項が帯状流の粘性散逸項と釣り合うことになる．つまり，単純な時間平均での理解では，乱流は帯状流の駆動に働いていると言える．一方で，帯状流と乱流との非線形相互作用項の時々刻々の変化を観察してみると，時間的に正負に激しく変動し，瞬時的には駆動にも減衰にも働く．実験や数値シミュ
レーションで乱流と帯状流が共存した状態が観測されている通り，こうした一見乱雑にみえる
時間変動も単にランダムな運動ではなく，帯状流を生成・維持する働きを持つはずである．その実効的な相互作用を抽出するには，過渡的な準線形解析や時間平均値での解析だけでは十分でないことが分かる．
射影演算子法は，確率的時系列データの間の相関を抽出する．そこで，上述のような一見乱雑にみえる乱流と帯状流の非線形相互作用の解析に利用することで，既存の時間平均解析では見えてこなかった動的相関が評価できるとの着想に至った．
3.2.2 長谷川 - 若谷方程式
プラズマ物理では，粒子速度分布関数を用いる運動論方程式から 1 流体に簡約化された電磁流体方程式まで，解析する物理現象の時間・空間スケールに応じて様々な基礎方程式が用いられる（石澤, 2014）．ここでは，微視的不安定性による乱流の駆動と非線形相互作用を介した帯状流の形成を表わす簡潔な方程式として，長谷川 - 若谷方程式（Wakatani and Hasegawa, 1984）に基づく解析を行う．一様背景磁場 B = B0 ẑ ，2 次元 (x, y) 局所乱流近似の仮定の下で，無次元化された長谷川 - 若谷方程式は以下で与えられる．
(3.3)
(3.4)

∂ω
− ∇φ × ẑ · ∇ω = C(φ − ω) − ν∇4 ω,∂t∂n∂φ− ∇φ × ẑ · ∇n + κ= C(φ − n) − ν∇4 n,∂t∂y

ここで，t は連続時間であり，φ(x, y, t), ω(x, y, t) = ∇2 φ, n(x, y, t) は静電ポテンシャル，渦度，プラズマの密度の揺らぎを表す．κ は背景プラズマ密度勾配を表すパラメータ，C ∝ ∇2 は断熱パラメータと呼ばれる有限の電気抵抗による電子の Boltzmann 応答からのずれを表す量で（3.3）（，3.4）は，背景密度あり，本研究では ∇ ∝ ∂y を仮定して，C = −ca ∂y2 により与えた．式
勾配を駆動源として不安定化する抵抗性ドリフト波乱流と，乱流から形成される y 方向に一様な流れである帯状流の相互作用のダイナミクスを記述する．そこで，帯状流の生成・崩壊メカ（3.3）の帯状流成分 (k = kx x̂, ky = 0) にニズムを分析するために，Fourier 変換した渦度方程式
関する式に着目した解析を行った．
(3.5)
(3.6)

dωk
= Ik − νkx4 ωk ,dtIk = −
δk  +k  ,k k  × ẑ · k  φk  ωk  .
k

k 

40

統計数理

第 71 巻 第 1 号 2023

非線形相互作用項 Ik は，三波結合条件 k  + k  = k を満たす乱流揺動 k  , k  が帯状流 k の時間（3.5）
発展に影響を与えることを示している．また，準定常乱流 (dt |ωk |2  = 0) を考えると，式から導かれるエンストロフィーの保存則より，非線形項は粘性 ν による散逸と釣り合うように，Re[ωk∗ Ik ] = νkx4 |ωk |2 > 0 となることがわかる．つまり前述の通り，単純な時間平均の描像からは，乱流相互作用による非線形項は帯状流の駆動に働くという理解が得られる．
3.2.3 連続時間射影演算子法の適用による乱流・帯状流相互作用の抽出
一方で，実際に渦度の帯状流成分 ωk やそれに対する非線形項 Ik の時系列データをプロットしてみると，乱流の時間スケールで激しく振動するデータが得られる（図 1 は実際にこれらの実部を図示した例であった）．そこで，目的変数を非線形項 f (t) := Ik (t)，説明変数を帯状流の渦度 g(t) := ωk (t) として連続時間射影演算子法を適用することで，非線形項の動的性質を抽出（a）に示すように，非線形項を，帯状流と相関してゆっくりと変動するすることを試みた．図 6と，相関を持たずに振動する無相関相関部分（Markov 項はこの場合小さく，主として記憶項）項に分離できている．乱流成分 k  , k  のダイナミクスで決まる非線形項はほぼ帯状流と無相関に振動している f (t) ∼ r(t) が，有限の記憶項を持つことが分かる．評価された記憶関数を図 6（b）に示す．記憶関数はほぼ正の実部を持つ減衰関数であり，指数関数 Γ(t) = γ/τ exp(−t/τ )でよく近似される．この非線形項と帯状流の相関を失うまでの記憶の時定数 τ ∼ 1 は，帯状流の変動の時間スケール（おおよそ図 6 の記憶項 M (t) と同程度の時間スケール ∼ 10）に比べて短いため，記憶関数を Γ(s) ∼ 2γδ(s) としてデルタ関数で表す Markov 近似で解釈できる．このとき，非線形項は Ik (t) ∼ −γωk (t) + r(t) として表され，γ > 0 であることから，記憶項が負の係数項として寄与することが示された．物理的には，一見乱雑に変動する乱流との非線形相互作
用項が帯状流と相関を持ち，その相関部分（記憶項）は帯状流の減衰に働いていることを意味す（3.5）を解く
る．加えて，無相関項を乱数モデルで置き換えた Ik の一般化 Langevin モデルと式ことで再現性テストを行った結果，帯状流のエネルギースペクトルをおおよそ再現することが確認された．こうして，記憶項は帯状流を減衰させるように働く一方，無相関項は帯状流の確率的駆動として働くことが示された．
上記の解析結果は，乱流と帯状流間の非線形相互作用に対する以下の様な物理解釈を示唆すると考えられる．まず，帯状流に対する相関項（記憶項）は乱流 - 帯状流相互作用を表し，乱流存在下では帯状流が典型的に阻害される傾向にあることを表している．一方，帯状流に対する

図 6．長谷川 - 若谷プラズマ乱流シミュレーションの時系列データセットから，帯状流に対する乱流非線形相互作用項 f (t) = Ik (t) と帯状流渦度 g(t) = ωk (t) に対して射影演算子
法を適用した解析結果 (kx = 0.3)．（a）射影演算子法による非線形項 f (t) の Markov
t
項 Ωg(t)，記憶項 M (t) = − 0 Γ(t − s)g(s)ds，無相関項 r(t) への分解を，ある 1 サ
ンプルについて例示．（b）評価された記憶関数 Γ(t) の実部・虚部．一般化第 2 種揺動
散逸定理（式（2.18））が数値的に成立していることも併せて確認された．

射影演算子法による統計的時系列データ解析とその応用

41

無相関項は乱流 - 乱流相互作用を表し，乱流渦の位相が偶然そろった時に起こる帯状流の確率的な駆動または減衰を表している．単純な時間平均解析ではなく，射影演算子法により非線形相互作用の動的性質を抽出することで得られたこれらの結果は，乱流 - 帯状流共存系における相互作用を一般化 Langevin 描像として捉える新しい解釈を提供する．
4.

おわりに

本研究では，統計的時系列データ解析としての射影演算子法について焦点を当てた．その定式化，および，時系列データ解析・モデリングへの適用という観点から，本研究で取り組む 3
つの問題を第 1.2 節で定義した．
（1）について，射影演算子として森の線形射影（有限ランク
第 2 節では，定式化に関する問題
射影）を用いた場合に，目的変数を任意関数とする定式化の拡張を行い，連続時間および離散時間射影演算子法に対する拡張された一般化第 2 種揺動散逸定理を新たに導出した．
第 3.1 節では，自己回帰モデルとの類似点・相違点の問題（2）について，SVAR モデルに基づ
く数値シミュレーションデータに射影演算子法を適用することで，両者の比較を行った．射影
演算子法により得られた記憶関数が自己回帰モデル係数を上手く捉えていること，導出された
一般化 Langevin 形式に基づいて元データの統計的性質の再現が可能であるという結果が得られた．射影演算子法は任意の目的・説明変数に対する定式化であるため SVAR モデルよりも広い適用範囲を持つこと，射影演算子法を自己回帰的に用いた場合は SVAR モデルと同等の相関を抽出可能であることを持つことが示された．
（3）に関連して，2 次元長谷川 - 若谷プラズマ乱流シ
第 3.2 節では，物理現象への適用性の問題
ミュレーションから得られる時系列データセットに射影演算子法を適用した．帯状流の生成・維持過程の理解を目的として，乱流との非線形相互作用項と帯状流との相関を解析した．その結果，記憶項が帯状流の減衰に働き，無相関項が帯状流の確率的駆動に働く一般化 Langevin 描
像として解釈できることを示した．以上の結果は，射影演算子法の実物理現象の応用事例の拡大，一般化 Langevin 描像としての物理解釈を与えるデータ解析としての有効性，確率的時系
列データ再現のためのモデリングへの応用性を示している．
本研究で説明した射影演算子法は，特定の物理現象に依らない，統計的時系列データセットに対する汎用的な解析手法である．本文中でも述べた通り，射影演算子法はモデルの当てはめではなく，射影と補射影の分離に基づく関数展開であり，時系列データの持つノイズの性質に依らない柔軟な解析が可能である．射影演算子法の解析プログラムはオープンソースとして提供されているため，興味を持った読者は，目的変数と説明変数の時系列データセットを用意す
れば，直ちに Markov 係数や記憶関数の評価が行える．目的変数や説明変数をどのように選ぶ
か，得られた一般化 Langevin 形式からどのような物理を抽出するかは解析者の嗜好に依るところだが，第 3 節で示したいくつかの事例が着想の起点になれば幸いである．
謝辞本研究の一部は JSPS 科研費 JP20K03892，JP21H04874，JP19K12212 の研究助成の下で行われました．本共同研究は統計数理研究所統計思考院共同研究スタートアップ 2021-11「射影演算子法に基づく時系列データの統計」の枠組みの下で行われました．
ソースコードおよびデータの可用性
射 影 演 算 子 法 の 汎 用 解 析 プ ロ グ ラ ム は オ ー プ ン ソ ー ス と し て ，https://github.com/

42

統計数理

第 71 巻 第 1 号 2023

smaeyama/mzprojection/ で公開されています．第 3.1 節に関する SVAR モデルシミュレーショ
ンコードおよび解析プログラムは，https://github.com/smaeyama/maeyama_toukeisuri_2023/
で公開されます．

注．
1)

注意深く区別すると，統計的定常性あるいは時間並進対称性（式（2.3））が成り立つ u(t) の任意関数 f (t) や g(t) とは異なり，無相関項はある時刻の u(t) だけでは決定されず，初期
時刻 t = 0 に補射影に入っていた部分 r(0) = Qf (0) を起点として，補射影に入り続ける部分として定義されているため，r(t) については統計的時間並進対称性が成立しない．式（2.14）は，ri (t) と初期値 g(0) に対する統計平均操作として定義される．
2)
記 憶 関 数 の 表 式 か ら Koopman 演 算 子 M を 除 去 す る 変 形 に つ い て 補 足 す る ．離よ り ，r(n) =散時間システムにおける無相関項と記憶関数の定義（ 式（2.26）中 辺 ）QMr(n − 1) = Mr(n − 1) − PMr(n − 1) = Mr(n − 1) + Γ(n) · g(0) であり，g(1)をかけて統計平均をとることで，r(n)g(1) = Γ(n) · gg(1) という式（2.26）最右辺の関係式を得る．ここで，無相関項の性質（式（2.23））r(n − 1)g(0) = 0 に M を作用させて Mr(n − 1)g(0) = Mr(n − 1)g(1) = 0 となることを用いた．連続時間シ
ステムにおける Liouville 演算子 (=時間微分演算子)Λ は，関数の積に作用した場合
Λ(F G) = (ΛF )G + F (ΛG) という分配則が成り立つが，離散時間システムにおけるKoopman 演算子 M は，関数の積もやはり関数 F (u)G(u) = (F G)(u) であることから，
M(F G) = (F G) ◦ N (u(n)) = (F G)(u(n + 1)) = F (u(n + 1))G(u(n + 1)) = (MF )(MG)となることに注意．

参

考

文

献

Akaike, H. (1973). Information theory and an extension of the maximum likelihood principle, Proceedings of the 2nd International Symposium on Information Theory (eds. B. N. Petrov and F.
Csaki), Akadimiai Kiado, Budapest.
Arnold, V. I. and Avez, A. (1968). Ergodic Problems of Classical Mechanics, Benjamin Inc., New Yorkand Amsterdam.
アーノルド, ウラミジール, アベズ, アンドレ (1972).『古典力学のエルゴード問題』（吉田 耕作 訳）, 吉岡書店, 京都.
Chorin, A. J. and Lu, F. (2015). Discrete approach to stochastic parametrization and dimension reduction in nonlinear dynamics, Proceedings of the National Academy of Sciences, 112(32), 9804–9809.
Chorin, A. J. and Stinis, P. (2006). Problem reduction, renormalization, and memory, Communicationsin Applied Mathematics and Computational Science, 1(1), 1–27.
Chorin, A. J., Hald, O. H. and Kupferman, R. (2000). Optimal prediction and the Mori-Zwanzig representation of irreversible processes, Proceedings of the National Academy of Sciences, 97(7),2968–2973.
Chorin, A. J., Hald, O. H. and Kupferman, R. (2002). Optimal prediction with memory, Physica D:
Nonlinear Phenomena, 166(3), 239–257.
Darve, E., Solomon, J. and Kia, A. (2009). Computing generalized Langevin equations and generalized
Fokker-Planck equations, Proceedings of the National Academy of Sciences, 106(27), 10884–
10889.
Español, P. (1995). Hydrodyanamics from dissipative particle dynamics, Physical Review E, 52(2),1734–1742.
藤坂博一 (1998).『非平衡系の統計力学』, 産業図書, 東京.

射影演算子法による統計的時系列データ解析とその応用

43

藤澤彰英 (2005). ゾーナルフロー研究の現状と展望, プラズマ・核融合学会誌, 81(12), 971–1006.
Gouasmi, A., Parish, E. J. and Duraisamy, K. (2017). A priori estimation of memory eﬀects in reducedorder models of nonlinear systems using the Mori-Zwanzig formalism, Proceedings of the Royal
Society A: Mathematical, Physical and Engineering Sciences, 473(2205), 20170385.
Granger, C. (1969). Investigating causal relations by econometric models and cross-spectral methods,
Econometrica, 37(3), 424–438.
Hald, O. H. and Stinis, P. (2007). Optimal prediction and the rate of decay for solutions of the Euler equations in two and three dimensions, Proceedings of the National Academy of Sciences,104(16), 6527–6532.
石澤明宏 (2014). 二流体およびジャイロ流体方程式系とその MHD 不安定性解析への応用, プラズマ・核融合学会誌, 90(4), 213–227.
川崎恭治 (2009).『非平衡と相転移 — メソスケールの統計物理学 — 』, 朝倉書店, 東京.
Kilian, L. and Lütkepohl, H. (2017). Structural Vector Autoregressive Analysis, Cambridge University
Press, Cambridge.
Lin, K. K. and Lu, F. (2021). Data-driven model reduction, Wiener projections, and the KoopmanMori-Zwanzig formalism, Journal of Computational Physics, 424, 109864.
Lin, Y. T., Tian, Y., Livescu, D. and Anghel, M. (2021). Data-driven learning for the Mori–Zwanzig
formalism: A generalization of the Koopman learning framework, SIAM Journal on Applied Dynamical Systems, 20(4), 2558–2601.
Lin, Y. T., Tian, Y., Perez, D. and Livescu, D. (2022). Regression-based projection for learning MoriZwanzig operators, http://dx.doi.org/10.48550/ARXIV.2205.05135.
Maeyama, S. and Watanabe, T.-H. (2020). Extracting and modeling the eﬀects of small-scale ﬂuctuations on large-scale ﬂuctuations by Mori–Zwanzig projection operator method, Journal of the
Physical Society of Japan, 89(2), 024401.
Meyer, H., Pelagejcev, P. and Schilling, T. (2020). Non-Markovian out-of-equilibrium dynamics: A
general numerical procedure to construct time-dependent memory kernels for coarse-grained observables, EPL (Europhysics Letters), 128(4), 40001.
Mori, H. (1965). Transport, collective motion, and Brownian motion, Progress of Theoretical Physics,
33(3), 423–455.
Mori, H. and Okamura, M. (2007). Dynamic structures of the time correlation functions of chaotic
nonequilibrium ﬂuctuations, Physical Review E, 76, 061104.
Munakata, T. (1994). Nonequilibrium statistical-mechanical approach to discrete-time dynamics, Physical Review E, 50, 2351–2354.
宗像豊哲 (1996).『物理統計学』, 朝倉書店, 東京.
村尾 博 (2019).『R で学ぶ VAR 実証分析時系列分析の基礎から予測まで』, オーム社, 東京.
Okamura, M. (2006). Validity of the essential assumption in a projection operator method, PhysicalReview E, 74, 046210.
岡村 誠 (2014). 射影演算子法の基礎とカオス・乱流への応用, 日本流体力学会誌 ながれ, 33(1), 11–22.
沖本竜義 (2010).『経済・ファイナンスデータの計量時系列分析』, 朝倉書店, 東京.
Parish, E. J. and Duraisamy, K. (2017). A dynamic subgrid scale model for large eddy simulations
based on the Mori–Zwanzig formalism, Journal of Computational Physics, 349, 154–175.
Rogers, B. N., Dorland, W. and M., K. (2000). Generation and stability of zonal ﬂows in iontemperature-gradient mode turbulence, Physical Review Letters, 85(25), 5336–5339.
Stinis, P. (2004). Stochastic optimal prediction for the Kuramoto–Sivashinsky equation, Multiscale
Modeling & Simulation, 2(4), 580–612.
Stinis, P. (2013). Renormalized reduced models for singular PDEs, Communications in Applied Mathematics and Computational Science, 8(1), 39–66.
Stinis, P. (2015). Renormalized Mori-Zwanzig-reduced models for systems without scale separation,
Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 471(2176),

44

統計数理

第 71 巻 第 1 号 2023

20140446.
鈴木増雄 (2000).『統計力学（現代物理学叢書）』, 岩波書店, 東京.
Wakatani, M. and Hasegawa, A. (1984). A collisional drift wave description of plasma edge turbulence,
The Physics of Fluids, 27(3), 611–618.
山本義隆, 中村孔一 (1999).『解析力学 I, II』, 朝倉書店, 東京.
Zwanzig, R. (1961). Memory eﬀects in irreversible thermodynamics, Physical Review, 124(4), 983–992.
Zwanzig, R. (1980). Problems in nonlinear transport theory, Systems far from Equilibrium (ed. L.
Garrido), 198–225, Springer, Berlin, Heidelberg.

Proceedings of the Institute of Statistical Mathematics Vol. 71, No. 1, 25–45 (2023)

45

Statistical Analysis of Time-series Data Using the Projection Operator
Method and Its Application
Shinya Maeyama1 and Fumikazu Miwakeichi2
1 Department of Physics, Nagoya University
2 The Institute of Statistical Mathematics

The projection operator method is a mathematical method developed in nonequilibrium statistical physics, where the time evolution of the response variable is split into
correlated and uncorrelated terms for the explanatory variables and are described as a
generalized Langevin equation. The projection operator method is practically applicable for the statistical analysis of time-series data. In the present study, we revisited its
formulation and extended it to arbitrary explanatory and response variables in continuous and discrete-time systems. We also compared the method with a structural vector
auto-regressive (SVAR) model as a time-series analysis having similar structures. The
projection operator method has a wider application range than the SVAR model andcan extract equivalent correlations with the SVAR model when auto-regressive explanatory and response variables are employed. As an application example, we applied the
projection operator method to plasma turbulence phenomena and explained zonal-flow
generation/maintenance processes in a generalized Langevin description. These examples
demonstrate the validity of the projection operator method in data analysis for physics
interpretation and modeling for reproducing statistical time-series data. We provide thedeveloped method in this study as an open-source Python code, allowing readers to use
it by calling a simple function.

Key words: Projection operator method, continuous/discrete-time system, turbulence, time series analysis.