〈総説〉
線形・非線形時系列解析とその応用 (2)
千 はじめに 姉妹論文では，まず、幾つかの人工的時系列信号に対 する伝統的時系列解析， とりわけスペクトル解析の結 果や，近年の非線形時系列解析， とりわけカオス時系 列解析における軌道不安定性 (orbital instability) を中 心に紹介した.また，時系列信号を間ヲ|いた場合の力 学系の位相的特徴の変化について紹介した.それに対 して， この論文では最初にリカレンスプロットの概要 を紹介し，人工的時系列信号にそれを施した 2 ， 3 の 例を紹介する.つぎに幾つかの人口的時系列信号に対 するフラクタル次元の概要の紹介 とりわけ樋口の方 法 (Higuchi's method) を人口的時系列信号に適用し た結果を紹介する.後半のフラクタル次元は，カオス 時系列信号の自己相似性 (self-similarity) という，カ オスの軌道不安定性と並び重要な特徴の l つである2
人工的時系列データのリ力レンスプロットの 概要と分析結果 既に姉妹論文の考察のところで指摘したように，線 形時系列解析の l っとしてのパワースペクトルの算出 が許容されるためには，時系列信号が定常性を持って いなければならない.また，非線形時系列解析， とり わけカオス時系列解析の l っとしてのリアプノフ指数 の算出が許容されるためには，時系列信号がエルゴー ド的でなければならない これに対して， これらの前提を必要とせず時系列信 号の多様な性質を描き出すことが可能な方法として， 近年リカレンスプロット (recurrenceplots，略して RP) が注目されている RP は， Eckmann et al.(1987) が提 唱したもので， 1990 年代以降多くの研究がなされて いる(例えば， Casdagli, 1997; 大ら， 2002; Fletcher, 2000; 平田， 2011; Hirata & Aihara, 2010;賓来ら， 2002; 里子 直 日 1 A 図 1 :ホワイトノイズ時系列データのリカレンスプロットlkeguchi etal., 1997; Ngamga et al., 2008; Romano et al.,
2004; Zbilut et al., 1998; 山田 ・ 合原， 1999; Webber & Zbilut, 1994). リカレンスプロ ッ トの描画法はきわめて単純なもの で，時系列信号旬(1) ， v(2) , v(3) ,…, v(N) の 2 時点聞の 距離 dij
=
1
v(i)-v(j)1,
) 咽 a E A ( に基づいて， dijがしきい値より小さい時に 2 次元の 座標 (i， j) に点をフ。ロッ トしたり， dijの値をグラデ ーションカラーで表現してすべての距離情報をフ。ロ ッ トしたりする. この論文で利用しているカオスソフト SCT では，後者を用いている. 例えば，姉妹論文でのホワイトノイズ， 一様乱数， 及びロジスティック写像のリカレンスプロットのカラ ー画像を示したものが図 l から 3 である.なお， '-'ｭ で用いたロジスティック写像は，次式による-x(η+1
)
=
b*
x(η)(1-
x(η)) ， b=4. (2) *愛知学院大学心身科学部心理学科 (連絡先) 〒470-0195 愛知県日進市岩崎町阿良池 12 E-mail:[email protected]千野直仁
0
.
0
1
0
.
1
. H H JD
1
1
0
0
1
0
0
9
0
0
8
0
0
7
0
0
6
0
0
5
0
0
4
0
0
3
0
0
2
0
0
1
0
0
1
.
6
1
.
4
1
.
2
0
.
8
0
.
6
0
.
4
0
.
2
0
100 200 300400500 6007∞ 80090010001100 図 2 一様乱数時系列データのリカレンスフ。ロット ちなみに，カラー表示は平面上の度数の違いをカラ ーグラデュエーションで表すものである.なお， これ らの図では，左上三角部と右下三角部で色が異なるが， 前者は後者のデータの対数を取ったものである 上記 3 つのリカレンスプロットのうち，最後のロジ スティック写像の場合， 全 2 者と比べて，よく見ると (右下 3 角部の点の集合の中に)部分的に斜め右上向 きの集合(あるいは，短い上向きの対角線分)(
s
h
o
r
t
u
p
w
a
r
d
d
i
a
g
o
n
a
l
segments) が数多く見られる. これは， カオスの l つの特徴である(例えば，平田，2011
, p.1
5
1
-
1
5
2
;
W
e
b
b
e
r
&Zbilut, 1994,
p
.
968
,F
i
g
.
4
)
.
つぎに， リカレンスフロットでロジスティック写イ象 以上に短い上向きの対角線分が数多く見られるエノン 写像 (Hénon map) のそれを見てみよう.エノン写像は， 2 次元の差分方程式であり (Hぬon， 1976) ，次式で表 4 2 1 8 6 4 2 0 2 4 6 8 1 2 4 1 1 0 0 8 0 0 B O D E -ｭ 刈-される. x(η+ 1) y(η+ 1) 1-αx(η)2+
y(η),
bx(η)1
E
0 0
.
0
0
0
1
0
.
0
0
1
0
.
0
1
0
.
1
1
1
0
0
1
0
0
9
0
0
8
0
0
7
0
0
6
0
0
5
0
0
4
0
0
3
0
0
2
0
0
1
0
0
1 00 200 300 400 500 600700 800 9001 0001100 (3)D
i
j
1
.
4
1
.
2
0
.
8
0
.
6
0
.
4
0
.
2
0
図 3 :ロジスティ ック写像のリカレンスプロット ここで， α= 1. 4， b=0.3 である.まず，原信号であ るエノン写像の x-軸座標の時系列を図4 に示す. エノン写像のリカレンスフロ ットの「短い上向きの 対角成分」をより見やすくするために，カラーではな くモノクロ画像で示したのが，図5 である. これに対 して， ランダム性もカオス性もない周期時系列信号， 例えば姉妹論文 (1 )のサインカーブ時系列データの 場合はリカレンスプロ ットはどうなるであろうか. こ のデータのリカレンスプロットは図6 のようになる. 図 6 の右下三角部を見ると，サインカーブの特徴を反 映して，規則的な模様がきれいに描かれていることが わかる 100 20 日 300 400 500 600 700 800 900 100日 tin噌 t 図 4 :エノン写像時系列データ -102一一5
5
0
0
5
0
0
0
4
5
0
0
4
0
0
0
3
5
0
0
3
0
0
0
2
5
0
0
2
0
0
0
1
5
0
0
1
0
0
0
5
0
0
100200 300 400500 600 700 800 900 図5 エノン写像のリカレンスプロッ ト0
.
0
0
1
0
.
0
1
0
.
1
Dﾜ
500 10001500 2000 2500 3000 3500 4000 4500 5000 55003
2
.
5
2
1
.
5
0
.
5
。千野直仁 0.0001 0.001 0.01 0.1 280 260 240 220 200 180 160 140 120 100 80 60 40 20 Dij
3
2
.
5
2
1
.5
0
.
5
。 20 40 60 80 100 120 140 160180200 220 240 260 280 図 7 :間引き後のサインカーブ時系列データのリカレンスプロット 一方，同信号を 100 時点づっ間 51いた場合のリカレ ンスプロット結果を示したのが 図 7 である. 姉妹論 文で述べたように この信号は間引き前の信号では遅 延座標が閉軌道となっていたものがトーラス状に変化 しており，図 7 のリカレンスプロットはトーラス状軌 道のそれであることに注意したい. 一方，サインカーブが時間とともに減衰し平衡点 ( equilibrium) に収束していくようなカーブではリカ レンスプロットはどうなるであろうか.まず図 8 は， 減衰線形振動子 (damped linearoscillator) の特別なケ ースを描いたものである. ここでの減衰線形振動子は 次式で表される 2 階微分方程式である d2x dーす +2α三 +(α2
+
b2)x 二 0，
α=
0.1,
b=
1. (4) αγ dt 図 9 は，減衰線形振動子の時系列データのリカレン スプロット結果を示す. この図を見ると右下 3 角部分 のほとんどは赤一色となっているが，わずかに同部分 の左端に近いところにサインカーブ時系列のような穴 目模様が見られることがわかる. さらに， この穴目模 様は右に行くほど、小さくなっていることもわかる こ れは減衰線形振動子の振動の特性を反映しているとみ れよう. 最後に，よく知られた間欠性カオス (intermittentchaos) (例えば， Aizawa & Kohyama, 1984; Kohyama, 1984; Manneville, 1980; Manneville & Pomeau, 1980) を
図 10 に示す. この図は， Aizawa
&
Kohyama (1984, p.847) の修正ベルヌイ写像 (themodified Bemoulli map)
を表す次式の 2 次元差分方程式のパラメ ータのうち， B=3，C= 0 としたケースである. x(n+ 1) = x(η) + 2(B -1)(1-2C)x(n)B + C, x(η+1) = x(η) -2(B -1)(1-2
*
C)(l -x(η))B -C. (5) 修正ベルヌイ写像時系列信号の特徴は，図にあるよ うに激しく振動するバースト層 (the bursting phase) と振動のないラミナ一層 (the laminarphase) から成 る点である.なお，彼らによれば， (5) 式のパラメー タのうち C がゼロでない小さな値をとる場合には， この時系列はある限られた領域で定常的な時系列とな り ， rvスペクトル (ther
v spectrwn) が観測できると いう.図 11 は，この信号のリカレンスプロットを示す.3
人工的時系列データのフラクタル次元の概要 と分析結果 一般に，カオスの幾つかの特徴(例えば，合原編， 2000) のうち，軌道不安定性と並び，重要な特徴とし てアトラクタのフラクタ Jv性 (fractal property) ，自己 -104一一γ
v
J I 、守 、 B P J d 凋抽 y nJ ι4 . ， A Mun hua 凋守 合 Jι A HUA7ιa 凋守 n hU AMU 4Etn Jι eB''B ， AHυoHUAHunHunHunHunHU 内Hu--EEm
-- 4
100 200 300 400 500 600 図 8 :減衰線形振動子時系列データ 1 E 0 0.000010.0001 0.001 0.01 0.1 19
0
0
8
0
0
7
0
0
6
0
0
5
0
0
4
0
0
3
0
0
2
0
0
1
0
0
100 200 300 400 500 600 700 800 900 700 Dij 図 9 :減衰線形振動子時系列データのリカレンスプロット 800 900 time t 43
.
5
3
2
.
5
2
1
.
5
0
.
5
0
千野直仁 X1 nud n MU 『 ，， ph urhda 凋守 内ヨリ A1ι ・ Ee h HU AHUh Hu h H υn H UA HU A H UA HU A HU 500 1000 1500 ~OOO 2500 3000 3500 4 即日 time t 図 10 : 修正ベルヌイ写像時系列データ 1 E 0.000010.0001 0.01 0.1
Dﾜ
1
.
6
3
5
0
0
1
.
4
3
0
0
0
1
.
2
2
5
0
0
2
0
0
0
0
.
8
1
5
0
0
0
.
6
1
0
0
0
0
.
4
5
0
0
0
.
2
。 500 1000 15002000250030003500 図 11: 修正ベルヌイ写像データのリカレンスプロ ット 相似性があげられる.軌道不安定性に対して，後者は 図形や軌道の幾何学的特徴を表す. ある図形が自己相 似的である場合，その図形のどんな部分を取り出して そこを拡大しても，その図形の全体的な特徴がみられ る.同様なことが時系列信号の軌道についても言える 場合，その軌道は自己相似性を持つという. 自己相似性は，各種のフラクタル次元(丘actal dimensions) により表現される.同次元については， これまでにボックスカウント次元 (box counting dimension)，相関次元 (correlation dimension) ，樋口の 方法 (Higuchi， 1987) などが提案されている.ちなみに， ボックスカウント次元は，容量次元 (capacity dimenｭ sion)，グリッド次元 (grid dimension) などとも呼ば れることがある(仔日えば， Russell et al., 1980). まず，最も単純なフラクタル次元について紹介する とつぎのようである.すなわち，通常の図形の次元は， 例えば通常の線は l 次元であり同平面は 2 次元などの ように正数値となる. この直感的な次元数をより正確 に定義するには，次のような関係に注意するとよい すなわち，ある図形を 7 に縮小した図形を用いて，も との図形を覆い尽くすために最低限必要な縮小図形の 数を N(けとすれば，両者の間には，つぎの関係が成 -106り立つ.
N
(
r
)
=
r
D 20 ￡一 〉〉 15 れJ 由 お 10 4コ k u v A ) J O G --5 0 4 5 6 log k with base 2 図 12 :樋口の方法によるロジスティック写像の フラクタ Jレ性 例えば，線分を i に縮小した線分で，もとの線分を 再現するには縮小線分は 2 つ二 21必要である.すな わち， この場合 D=l であり， したがって，フラクタ ル次元は 1 となる. しかし，図形の中には線分から成 り立っているにも関わらず 次元が正数値とならない ものが存在する. よく知られたものに，コッホ曲線 (Koch curve) やシェルピンスキーのギャスケット (SierpinskiGasket) などがある.ちなみに，前者のフ ラクタル次元は1. 26 シェルピンスキーのギャスケッ トの同次元は1. 585 であるつぎに相関積分 (colTelation integral) について述べ
る.相関積分は， Grassberger & Procassia(1983a, b) に よるもので，相関次元を求めるためには相関積分を求 めなければならない.相関積分は m 次元空間上で の軌道の各点 v(i) ， i=l ，・・・ ， N に対して N
C
rn(γ)=J!弘元22ン [r-
1
v
(
i)-V
(
j
)
1
.
]
(7) 1子正3 ここで ， I(t) はヘビサイド階段関数 (Heaviside-step functio引で， 1inU 〆 'EI 〈、， l 、 一一 、 、lJ ノ 4 F U 〆， a ，、 町、 r l(
t 三 0)(
t
:
:
;
0
)
つまり，相関積分は m 次元空間上での軌道の各 点で距離 T 内に入る軌道の点の数の平均(の極限)で ある. Crn( γ) 氏 γ v(rn) ヲ (8) (6) の時， (8) 式のスケーリング指数 ν(m) は相関指数 (correlationexponent) と呼ばれる. これより， (8) 式 の両辺の対数を取ることにより得られる関係式を考慮、 して， logr を横軸に， logCη吋7 とき，あるスケ←ルの範囲町三 γ 三 r2 で ν(m) なる勾 配の直線になっていれば，その範囲内で自己相似性が 成り立つものと見れる 15 10 kd N 由的問 2ztE( 五 )JOG-10 5 -10 0 内t e e d a hU 6 m W L K G U O 10 12 図 13 :樋口の方法によるサインカーブ時系列の フラクタ Jv↑生 (7) 式の相関積分を一般化して，距離 T の超立方体を η 個考え，C
;
;
(
r
)
=
LP;'
(9) i=l とするとき， (9) 式が C;;(γ)α rVq(rn) う 、、 ，，， r A U ' E EA ( と書けたとする.上式の vq(m) は ， q 次の相関指数と呼ばれる (Grassberger
&
Procassia, 1983a). ここで，q=2 の場合，相関次元と呼ばれる. 相関次元も，既述のフラクタル次元同様，多くの場 合非整数の値をとるが 整数値をとる場合もある.例 えば，記述のロジスティック写像などの非可逆的写像 (noninvertiblemaps) では，相関次元は整数値となるが， エノン写像などの散逸的写像 (dissipative maps) では 同次元は非整数となる(例えば， Sprott, 2003). ちな みに，ロジスティック写像では相関次元は 1 ，エノン 写像では 1.220 土 .036 となることが知られている
千野直仁 ., Xn として，新たな時系列を次のように生成す る: XmぅXm+k ， Xm+2k ぅ ，X
m
+ 円!I':.]k' m 二 1 ぅぅ k つぎに，樋口の方法ではの長さを以下のよう に定義する: L ん 1 111 1 〆 1111山一山
十 Z + 中 山一
T
hL
凶 f it -、 III t L ここで， (11)式の[・l 記号は，ガウスの整数 (Gaussian integer) である 最後に，われわれは時間間隔 k ごとに，上記 Lm(k) の k 組に亘る平均長 L(k) を計算する.もし， これが， L(k) α k-D， (12) ならば，信号は次元 D のフラクタルであると呼ぶ. 最後に，相関次元の場合と同様， (1 2) 式の両辺の対 数を取り， logk を横軸に， logL( めを縦軸にフ。ロット した時， これが-D なる勾配の直線になっていれば， 信号は自己相似性が成り立っと見れる. この時，当該 信号は次元 D のフラクタルと言える なお，フラクタル次元 D は，時系列信号のパワー スペクトル P(f)cx FCt における α との聞につぎのよう な関係があることが知られている (Berry， 1970) : D=
(5-α)/2ぅ 1<α<3. (13) Berry (1 979) は，上式での D の違いにより， 1 D=l (α=3) の時，周辺的フラクタル (marginal 仕actal) で，ほとんどなめらかなカーブ， 2 D=1.5 (α=2) の時 ブラウン的フラクタノレ (Brownian fractal) で 1 次元ブラウン運動 (Brownian motion) の時， 3 D=2 (α=1) の時，極担なフラクタル (extreme 自'ac叫で， j ゆらぎ (1 noise) の時， の 3 者を区別している.また，樋口 (1989 )は D土 1 (α 三 3) ， Mandelbrot (1 977) は D=2 (0 三 α< 1) として いる 以降，上記の樋口のフラクタル次元を ， D を相関次 元 D2 と区別するために DH と表記することにする. なお，上記の関係から，フラクタル次元で、は白色ノイ ズと十を区別できない また Cりt et a.l(2013) は， 時系列信号が周期的な場合， log2L(劫 vs.log2 k のグラ フには k の大きな区間で振動が見られることを報告し ている.さらに，彼らは， (1)周期的要素が入り込むこ とを避けるには ， k (max) をできるかぎり小さく取る， (2)低周波数の周期的成分を同定したければ， k (mαx) をできるかぎり大きく取る ことを勧めている. 図 12 は，第 2 節で述べたロジスティック写像につ いて，樋口の方法による log2 L( 劫 vs. log2k のグラフ を示したものである.図から明らかなように，グラフ は，ほぼ直線となっており， この時系列信号がフラク タル特性を有していることがわかる.ちなみに，樋口 のフラクタル次元 DH の値は 2.00035 とほぼ 2 となっ ている.なお，姉妹論文で示したこの時系列信号のパ ワースペクトルがホワイトノイズなどと同桔ミフラット になっているので ロジスティック写像では上記 α=0 と見れる. したがって， ここでの DH は Mandelbrot の 指摘するフラクタル次元の値と一致する 同じく姉妹論文で取り上げたホワイトノイズ，及び 一様乱数信号でも，同様な結果が得られる.すなわち， それらの log2 L(k)vs. log2k のグラフは共にほぼ直線 となり DH の{直も )1頃に， 1.99806, 1.99959 となり， 両者ともほぼ 2 である. つぎに，同じく第 2 節で述べたサインカープ時系列 について log2 L(k)vs. log2k のグラフを示したのが， 図 13 である.図から明らかなように， この場合図は k の大きな領域で振動しており， この結果は Coyt et al. (2013) の指摘と一致している. サインカ←ブ時系列は， フラクタ }vt生を持たないことがわかる. References 合原一幸編 (2000). カオス時系列解析の基礎と応用 産業 図書Aihara, K., Ikeguchi, T., & Matsumoto, G. (1998). Deterministic nonlinear dynamics of a forced osciIIation experimentaIIy observed witha squid giant axon. International Jourηαlof Chaos Theory and itsApplicαtions ， 3, 5-20.
Aizawa, Y., & Kohyama, T.(1984). Asymtrotically nonｭ stationary chaos. Progressof Theoretical Physics, 71, 847 850. 大聖一郎・和田充雄・山口明宏・広奥暢 (2002). バイオメ カニズム 16 ， 285-298. Casdagli, M. C.(1997).RecuITence plots revisited. Physica D, 108,12-44.
Coyt, G. G., Diosdado, A. M., Balderas, 1.A., deI Rio Correa, J.
L., & Brown, F.A. (2013). Revista Mexicana de Fisica S,
59, 1-6.
Eckmann, J. -P., 0 lIiffson Kamphol叫， S., & RueIle, D. (1987) Recurrence plots of dynamical systems.Europhysics Letters, 4, 973-977.
the Acoustical Society0/ Amerika, 108, 82 卜826.
Grassberger, P., & Procaccia, .1(1983a). Characterization of strange attractors. PhysicalReνiew Letters, 50, 346-349.
Grassberger, P., & Procaccia, T.(l983b). Measuring the strangeness of strance attractors. Physica D, 9, 189-208.
Hénon, M. (1976). A two-dimensional mapping with a strange attractor. Communications inMα thematical Physics, 50,
69-77.
Ikeguchi, T., & Aihara, K. (1997). Lyapunov spectralanalysison random data. Jnternational Journal 0/ Bifurcationand Chaos, 7,1267-1282. Higuchi, T. (1988). Approach to an irregular til11e series on the basis ofthefractaltheory. Physica D, 31, 277-283. 樋口知之(1 989). 時系列のフラクタル解析 統計数理 37う 209-233. 平田祥人 (2011). リカレンスフ。ロット:時系列の視覚化を 超えて 数理解析研究所講究録1768 ， 150-162
Hirata, Y.,
&
Aih創刊K. (2010). Devaney's chaos on recurrence plots. PhysicalReνiew E, 82, 036209. 賓来俊介，山田泰司，合原一幸 (2002). 同方向リカレンス プロットによる決定論的解析電学論 C, 122,141-147. Kohyama, T. (1984). Non-stationarity of chaoticl110tions in an area preservingl11apping. Progress 0/ Theoretical Physics, 71, 1104-11 07.Manneville, P. (1980). Interl11ittency, self-similarity andlIf
spectrum in dissipative dynamical systems. Journalde
Physique, 41,1235-1243.
Manneville, P., & Pomeau, Y.(1980). Different ways to turbulence in dissipative dynal11ical systems. Physica D, 1, 219-226.
Ngamga, E.J., Nandi, A., Ramaswamy, R., Romano, M. C., Thiel, M., & Kurths, J.(2008).Recurrences of strange
attractors.Pramαna， 70,1039-1045.
Romano, M. C., Thiel, M., Kur・ths， J., & von Bloh, W. (2004). Multivariate recurrence plots, Physics Letters A, 330, 214 223
Russell, D. A., Hanson, J.0.,
&
Ott, E. (1980). Dimension of strangeattractors.Physicα1 Reνiew Letters, 45, 1175-1178. Sprott, J.C. (2003).Chαos and Time-SeriesAnalysis. New York:Oxford University Press.
Webber, C. L. Jr., & Zbilut, J. P. (1994). Dynamical assessment of physiological systems and states using recurrence plot strategies. Journal0/ AppliedPhysiology, 76, 965-973.
山田泰司，合原一幸 (1999) リカレンスフ。ロットと 2点間
距離分布における非定常時系列解析 電子情報通信学
会論文誌， J82A, 1016-1028.
Zbilut, J. P., Giuliani, A., Webber C.L., Jr., & Colosimo, A. (1998).Recurrence quantificationanalysis instructureｭ function relationships of proteins: an overview of a general methodology applied to the case of TEM-l ﾟ-lactamase Protein Engineering, 11,87-93.
Bulletin ofThe FacultyofPsychological& Physical Science, No. 11,101-110,2015
L
i
n
e
a
r
and
N
o
n
l
i
n
e
a
r
Time
S
e
r
i
e
s
A
n
a
l
y
s
e
s
and t
h
e
i
r
A
p
p
l
i
c
a
t
i
o
n
s
(
2
)
Naohito CHINO Abstract
Recurrenceplot and fractal dimensionanalysisare introduced, both ofwhich can beapplicabletovarious time series signals. A1though autocorrelationsandpowerspectra are valid only if they are stationary, they are well known indicators fo1'detectingproperties of these signals in the linear time series analysis.Similar匂，Lyapunov
exponents aremeanmg乱11 if the time series under study is ergodic, which arefrequentlyusedtoexamine
whether it is chaotic in the nonlinear time series analysis. 1n contrast
, recun・ence
plot andfractal dimensions have no 01'few assumptions in app1ication to time series signals. Resu1ts of the application of recurrence plot tosome artificialtime series signals show that, for example, it can certainly be an indicator of chaos. Fractal dimension analysis which is used to detect geometrical property of chaos is alsoa relatively robust method and it is valid when it isappliedto time seriessignalswithself-similarity. Applications of Higuchi's method, which is one of the methodsfor detectingfractal dimension, to two artificialtimeseriessignals show thathis method enables us to detect fractal property of thesesignalswhen properlyapplied.
Keywords and phrases: chaos, recu汀enceplot, fractal dimensions, Higuchi's method
