JP2021064323A - 時系列データ解析装置及び時系列データ解析用プログラム - Google Patents

Abstract

【課題】解析対象のデータが離散系であろうと連続系のデータであろうと、適切に特にカオス性解析を行うことが可能な時系列データ解析装置を提供する。【解決手段】ｎを正の整数として、データ長がｎ＋１の時系列データ｛ξ0，ξ1，ξ2，・・・，ξn｝について、ξk∈Ｘiとなる確率分布ｐ（ｉ）、ξk∈Ｘi、ξk+1∈Ｘjとなる同時確率分布ｐ（ｉ，ｊ）を求める確率分布算出手段２０１と、条件付確率ｐ（ｊ｜ｉ）を求める条件付確率算出手段２０２と、ｍ個の測度ｂjの各クラスが｛ｂ1，ｂ2，ｂ3，・・・，ｂm｝であり、各クラスの最大値がｂmax、であるとき、確率ｑ（ｉ，ｊ）を求める確率算出手段２０３と、条件付確率に関する修正エントロピーｓ´（ｉ）を、確率分布ｐ（ｉ）と、条件付確率ｐ（ｊ｜ｉ）と、確率ｑ（ｉ，ｊ）を用いて求める修正カオス尺度算出手段２０４とを備える。【選択図】図５

Description

この発明は、時系列データ解析装置及び時系列データ解析用プログラムに関するものである。

時系列データがどのような特徴を有しているのかを把握するために、様々な手法が用いられる。例えば、時系列データについてカオス度合を把握するための定量化法として、リアプノフ指数と呼ばれる指標がよく用いられている。リアプノフ指数を求める場合は、データ生成源の情報（離散系の差分方程式や連続系の微分方程式等）が既知である必要がある。データ生成源の情報が未知である場合には、大量のデータから推定する手段が与えられているが、所謂「埋め込み次元の推定」に関する処理は、次元毎に系の様子を見ながら調整を繰り返す必要があり煩雑であるにもかかわらず、必ずしも一意に得られるものではないという問題がある（非特許文献５、６、７）。

近年になって情報理論に基づく「カオス尺度（Chaos Degree）」と称される指標が提案された（非特許文献１）。このカオス尺度とリアプノフ指数は極めて類似の挙動を示すことが知られている（非特許文献２、３、４）。図１には、ロジスティック写像（ａ：３．５〜４．０）のカオス尺度（分割数ｍ＝１００）とリアプノフ指数の値の変化を示している。この図１から、カオス尺度とリアプノフ指数カオス尺度は、類似の変化を行うことが明らかである。カオス尺度は、データのみから一意に計算可能であるという特徴を有している。

近年になって、本願発明者らは、カオス尺度とリアプノフ指数の数理学的関係性を明らかにした（非特許文献８、９、１０）。

＜カオス尺度の定義＞
データ長がｎ＋１の時系列データを
｛ξ₀，ξ₁，ξ₂，・・・，ξ_n｝・・・（１）
と表す。
データは、差分方程式τ

により、

として生成されているものとする。

上記のξ_kが含まれている区間Ｉをｍ個の区間に等分割する。この分割区間をＸ_i（ｉ＝１，２，・・・，ｍ）とすると、Ｉについて次の式（２）が成り立つ。

ここで、ξ_k∈Ｘ_iとなる確率分布ｐ（ｉ）と、ξ_k∈Ｘ_i，ξ_k+1∈Ｘ_jとなる同時確率分布ｐ（ｉ，ｊ）を次の式（３）、（４）の通りに算出する。

上記において＃［（条件式）］は、（条件式）を満たす数を意味する。

上記のように確率分布ｐ（ｉ）と同時確率分布ｐ（ｉ，ｊ）が求まると、カオス尺度Ｈは以下の式（５）または式（６）として定義される。

ただし、０ｌｏｇ０＝０とする。

上記のカオス尺度Ｈの最小値は０であるが、最大値は分割数ｍに依存してｌｏｇｍである。

カオス尺度Ｈの値は大きいほどカオス性が大きい（複雑性が高い）ことを意味する。即ち、カオス尺度Ｈに関する解釈は、次のようにまとめることができる。

上記のカオス尺度Ｈを、（１）単調増加データ、（２）周期データ、（３）ノイズが加わった周期データ、（４）正規分布に従うランダムデータ、（５）一様ランダムデータ、のそれぞれについて計算した結果、図２に示すようになった。

上記カオス尺度Ｈの計算は次の通りに行った。即ち、各データのデータ長をｎ＋１＝１０１、カオス尺度の分割数をｍ＝５、データの値域［０，１］をｍ＝５で等分割した区間を分割区間として計算を行った。この場合のカオス尺度の最大値はｌｏｇｍ＝ｌｏｇ５＝１．６０９５であった。よって（５）の一様ランダムデータのカオス尺度が最大値に近く、カオス性（複雑度）が高いことを確認できる。

図２には、上記（１）〜（５）のデータについてカオス尺度Ｈと共に分散ｓ²を計算して示してある。この図２からは、直感的に、（１）と（２）のデータが複雑とは言えず、（３）のデータがやや複雑さを帯びており、（４）と（５）のデータがランダムデータ故に複雑であると結論できる。この図２では、（１）〜（５）へ移行するに従ってデータの複雑さが増すように並べており、カオス尺度Ｈの値が徐々に大きくなっており直感的に複雑度が増していることが分かる。

一方、分散ｓ²は、以下の式で表わされる。

上記分散ｓ²は、暫し、「データのばらつき度合を示す」と説明される。分散は、データの平均値からの距離の二乗平均であって、データの乱雑さとは異なるものであることが図２に示されている値からも理解できる。即ち、データのデタラメさ、ランダムさ、乱雑さ、複雑さ、カオス性などを定量化して議論する場合には、分散を用いるよりもカオス尺度を用いる方が適切であることが分かる。

＜カオス尺度とリアプノフ指数＞
リアプノフ指数λは、不変測度（確率密度関数）ρを用いて次のように表わされる。

本願の発明者らは、カオス尺度の分割をΔｘ≡｜Ｘｉ｜（ｊ＝１，２，・・・，ｍ）と表したとき、カオス尺度はΔｘ→０の極限において、

の関係にあることを示した（非特許文献９、１１）。ここで、Ｄ（ｘ）は非負の関数である。よって、ＨΔ_x→₀≧λの関係にある。このＤ関数の詳細は非特許文献８、９に示す通りである。｜τ´（ｘ）｜におけるＤ関数の上限（Ｄ_sup）と下限（Ｄ_inf）を図３に示す。

このように、本願発明者らは、カオス尺度とリアプノフ指数の数理学的関係性を明らかにした（非特許文献８、９、１０）。

本発明は、上記のような事情に鑑み、カオス尺度とリアプノフ指数の差分を極めて小さくするように修正を行った新たなカオスに関する指標（「修正カオス尺度」と称する。）を用いて、与えられたデータのみからリアプノフ指数に極めて近似した値を一意に得るようにする。

また、データ生成源である写像が既知であり、かつ写像が保測写像である特殊な場合には、カオス尺度とＫＳエントロピーが同値を示すことが分かっている。そこで、本発明は、ＫＳエントロピーの推定にも適用可能となるものである。

大矢雅則，原利英，"数理物理と数理情報の基礎，"近代科学社，２０１６年 K.Inoue, M.0hya, K.Sato, "Application of chaos degree to some dynamical systems," Chaos, Solitons& Fractals, Volume 11, Issue9, July2000, Pages 1377-1385 K.Inoue, "Basic properties of entropic chaos degree in classical systems," Internationa1 journal on information 16(12(B)), December2013,8589-8596 井上哲，"カオス尺度における準周期軌道の取り扱い，"日本応用数理学会論文誌，ＶＯＬ．２５，Ｎｏ２，２００５年，ｐｐ．１０５−１１５ 岡田大樹，梅野健，"新たな非線形時系列解析の手法―移動最大リアプノフ指数線によるカオス解析―，"レーザー学会，Ｖｏｌ．４３，Ｎｏ６（２０１５）１９９３，ｐｐ．１１７−１３４． A. Wolf, J. B. Swift, H. L. Swinney and J. A. Vastano, "Determining Lyapunov Exponents from a Time Series," Physica D, Vol. 16, No3,1985, pp285-317. M. T. Rosenstein, J. J. Collins and C. J. De Luca, "A Practical Method for Calculating Largest Lyapunov Exponents for Small Data Sets," Physica D, Vol.65,1993, pp.117-134. 奧富秀俊，真尾朋行，"カオス尺度とリアプノフ指数の数理的関係性について，"電子情報通信学会，２０１７年度第３回ＣＣＳ研究会（２）， November,2017 奥富秀俊，真尾朋行，"カオス尺度の極限値の算出について，"日本応用数理学会２０１８年度年会，応用カオス（１）−３， September2018 真尾朋行，奧富秀俊，梅野健，"カオス尺度の計算方法について，"日本応用数理学会、２０１８年度年会，応用カオス（１）−４， September2018 奧富秀俊，"カオス尺度とＫＳエントロピーの関係についての考察，"電子情報通信学会，２０１８年度第１回ＣＣＳ研究会（２）， June,2018 真尾朋行，奧富秀俊，"心拍間隔データのカオス尺度と自律神経活動の関連について，"日本応用数理学会２０１６年度年会，応用カオス（２）−１， September2016 真尾朋行，奥富秀俊，"心拍変動のカオス性と副交感神経活動との関係について，"電子情報通信学会，第３回ＣＣＳ研究会（１）， November,2017

本発明は、解析対象のデータが離散系であろうと連続系のデータであろうと、適切に特にカオス性解析を行うことが可能な時系列データ解析装置を提供する。

実施形態に係る時系列データ解析装置は、ｎを正の整数として、データ長がｎ＋１の時系列データ｛ξ₀，ξ₁，ξ₂，・・・，ξ_n｝が、写像τによりξ_k＝τ（ξ_k-1）＝τ^k（ξ₀），ｋ＝１，２，・・・，ｎとして生成されるという仮定の基に、ξ_kが含まれる区間Ｉ≡［ａ，ｂ］をｍ個の区間に等分割した分割区間Ｘ_i（ｉ＝１，２，・・・，ｍ）とについて、

であり、ξ_k∈Ｘ_iとなる確率分布ｐ（ｉ）を求め、ξ_k∈Ｘ_i、ξ_k+1∈Ｘ_jとなる同時確率分布ｐ（ｉ，ｊ）を求める確率分布算出手段と、条件付確率ｐ（ｊ｜ｉ）を求める条件付確率算出手段と、ｍ個の測度ｂ_jの各クラスが｛ｂ₁，ｂ₂，ｂ₃，・・・，ｂ_m｝であり、各クラスの最大値がｂ_max、であるとき、測度ｂ_iの最大値ｂ_maxに対する割合をｑ（ｉ）としたとき、同時測度ｂ_i,jについては、最大値ｂ_maxに対する割合である確率ｑ（ｉ，ｊ）を求める確率算出手段と、条件付確率に関する修正エントロピーｓ´（ｉ）を、第１の情報量［−ｌｏｇｐ（ｊ｜ｉ）］から第２の情報量［−ｌｏｇｑ（ｉ，ｊ）］を引いた情報量の平均と定義するとき、前記条件付確率に関する修正エントロピーｓ´（ｉ）の平均として定義された修正カオス尺度を、前記確率分布算出手段により算出された確率分布ｐ（ｉ）と、前記条件付確率算出手段により求められた条件付確率ｐ（ｊ｜ｉ）と、前記確率算出手段により求められた確率ｑ（ｉ，ｊ）を用いて求める修正カオス尺度算出手段とを備えることを特徴とする。

ロジスティック写像のカオス尺度とリアプノフ指数の値の変化を示す図。 カオス尺度Ｈを、（１）単調増加データ、（２）周期データ、（３）ノイズが加わった周期データ、（４）正規分布に従うランダムデータ、（５）一様ランダムデータ、のそれぞれについて計算した結果を示す図。 写像をτとしたとき、｜τ´（ｘ）｜におけるＤ関数の上限（Ｄ_sup）と下限（Ｄ_inf）を示す図。 本発明の実施形態に係る時系列データ解析装置のブロック図。 本発明の実施形態に係る時系列データ解析装置の内部機能を示すブロック図。 本発明の実施形態に係る時系列データ解析装置による計算法１を用いて修正カオス尺度を算出する動作のフローチャート。 ｍ分割した写像τの１区間||Ｘ_i||について、写像後の値がどのように変化するかを表す計算を行う計算工程の説明図。 本発明の実施形態に係る時系列データ解析装置による計算法２を用いて修正カオス尺度を算出する動作のフローチャート（前半）。 本発明の実施形態に係る時系列データ解析装置による計算法２を用いて修正カオス尺度を算出する動作のフローチャート（後半）。 本発明の実施形態に係る時系列データ解析装置による計算法２によりロジスティック写像から得られた修正カオス尺度と、従来のカオス尺度、更にリアプノフ指数の波形を示す図。 本発明の実施形態に係る時系列データ解析装置による計算法２によりロジスティック写像から得られた修正カオス尺度と、従来のカオス尺度、更にリアプノフ指数の波形を示す図。 計算法３により修正カオス尺度を算出する過程の説明を行うためのタイミングチャート。

以下添付図面を参照して、本発明の実施形態に係る時系列データ解析装置及び時系列データ解析用プログラムを説明する。各図において同一の構成要素には、同一の符号を付して重複する説明を省略する。図４に実施形態に係る時系列データ解析装置１００のブロック図を示す。時系列データ解析装置１００は、クラウドコンピュータ、サーバコンピュータ、パーソナルコンピュータ、その他のコンピュータにより構成することができる。

時系列データ解析装置１００は、ＣＰＵ１０１が主メモリ１０２のプログラムやデータに基づき演算を行うものである。ＣＰＵ１０１には、バス１０３を介して外部記憶装置１０４が接続されており、外部記憶装置１０４には、時系列データ解析用プログラムが記憶されている。ＣＰＵ１０１が外部記憶装置１０４から時系列データ解析用プログラムを主メモリ１０２へ読み出してこのプログラムを実行することにより時系列データ解析装置として機能する。

バス１０３には、外部記憶装置１０４以外に時系列データ供給部１０５が接続されている。時系列データ供給部１０５は、外部のセンサなどからリアルタイムで時系列データを取り込み保持するものとすることができ、或いは、外部の何らかの装置などが収集した時系列データを取り込み保持したものとすることができる。更に、収集した時系列データを記憶した媒体がセットされることにより、時系列データを保持し供給可能となっている装置であっても良い。更に、上記の構成を全て備えたものであっても良い。いずれにしても、ＣＰＵ１０１が時系列データ解析用プログラムを実行して時系列データの解析を行う場合には、時系列データはこの時系列データ供給部１０５から供給される。

バス１０３には、結果出力部１０６が接続されている。結果出力部１０６は、表示装置やプリンタなど、時系列データ解析装置１００において処理した結果を出力する装置とすることができる。また、結果出力部１０６は、時系列データ解析装置１００において処理した結果を記憶する媒体でもよく、更に、回線などを介して処理の依頼者（クライアント）へ処理結果を送信などする装置であっても良い。

外部記憶装置１０４に記憶されている時系列データ解析用プログラムが実行されることにより、図５に示される各手段が実現される。即ち。時系列データ解析装置１００は、図５に示されるように、確率分布算出手段２０１、条件付確率算出手段２０２、確率算出手段２０３、修正カオス尺度算出手段２０４を具備している。

本実施形態の時系列データ解析装置１００では、以下に示す＜修正カオス尺度＞を得るものであり、これにより式（７）に示したＤ関数に相当する部分を小さく抑えるものであり、即ち、カオス尺度とリアプノフ指数の差を小さく抑え、リアプノフ指数の推定値として機能させるものである。これによって、ＫＳエントロピーの推定値としても機能させるようにする。

確率分布算出手段２０１は、ｎを正の整数として、データ長がｎ＋１の時系列データ｛ξ₀，ξ₁，ξ₂，・・・，ξ_n｝が、写像τによりξ_k＝τ（ξ_k-1）＝τ^k（ξ₀），ｋ＝１，２，・・・，ｎとして生成されるという仮定の基に、ξ_kが含まれる、区間Ｉ≡［ａ，ｂ］をｍ個の区間に等分割した分割区間Ｘ_i（ｉ＝１，２，・・・，ｍ）とについて、

であり、ξ_k∈Ｘ_iとなる確率分布ｐ（ｉ）を求め、ξ_k∈Ｘ_i、ξ_k+1∈Ｘ_jとなる同時確率分布ｐ（ｉ，ｊ）を求めるものである。

既に、＜カオス尺度の定義＞において述べたように、確率分布ｐ（ｉ）は式（３）により求められる。ここでは、データ長がｎ＋１の上記時系列データ｛ξ₀，ξ₁，ξ₂，・・・，ξ_n｝を軌道と称する。この軌道ξ_kについて、ξ_k∈Ｘ_i（ｉ＝１，２，・・・，ｍ）を満たす測度（式（３）における右辺の分子）をｂ_iで表す。

上記で、全測度は、定義より次の通りである。

よって、各クラスの生起確率ｐ（ｉ）は、次のように記載することができる。

上記式（１０）は、既述の式（３）に等しい。確率分布算出手段２０１は、上記式（１０）に示した各クラスの生起確率ｐ（ｉ）を算出するものである。

＜カオス尺度の定義＞においては、式（４）に、ξ_k∈Ｘ_i、ξ_k+1∈Ｘ_jとなる同時確率分布ｐ（ｉ，ｊ）を示した。ここでは、ξ_k∈Ｘ_i（ｉ＝１，２，・・・，ｍ）かつ、ξ_k+1∈Ｘ_j（ｉ＝１，２，・・・，ｍ）を同時に満たす測度（式（４）の右辺の分子）を、同時測度と称して、ｂ_i,jとすると、以下の式（１１）により表すことができる。

上記の同時測度ｂ_i,jは、測度ｂ_iのうち、次の時刻にｊ番目のクラスに移る測度であるから、次の式（１２）に示す関係にある。

上記により、同時確率ｐ（ｉ，ｊ）及び条件付確率ｐ（ｊ｜ｉ）は、次の式（１３）、（１４）により与えることができる。

上記の式（１３）により示される同時確率ｐ（ｉ，ｊ）は、確率分布算出手段２０１が求めるものであり、式（４）と等しいものである。このように上記条件付確率算出手段２０２は、上記確率分布算出手段２０１が求めた確率分布ｐ（ｉ）と同時確率分布ｐ（ｉ，ｊ）から条件付確率ｐ（ｊ｜ｉ）を求めるものである。

確率算出手段２０３は、ｍ個の測度ｂ_jの各クラスが｛ｂ₁，ｂ₂，ｂ₃，・・・，ｂ_m｝であり、各クラスの最大値がｂ_maxであるとき、測度ｂ_iの最大値ｂ_maxに対する割合をｑ（ｉ）としたとき、同時測度ｂ_i,jについては、最大値ｂ_maxに対する割合である確率ｑ（ｉ，ｊ）を求めるものである。この確率ｑ（ｉ，ｊ）については、＜エントロピー＞と＜修正エントロピー＞と＜条件付確率ｐ（ｊ｜ｉ）に関するエントロピー＞と＜条件付確率ｐ（ｊ｜ｉ）に関する修正エントロピー＞を説明する必要がある。

＜エントロピー＞
エントロピー（シャノン・エントロピー）ｓは、確率分布ｐ（ｉ）に関する情報量式（１５）の平均値として式（１６）、（１７）により求めることができる。

＜修正エントロピー＞
本実施形態では、修正エントロピーｓ´を定義する。式（８）のｂ_iについて最大測度をｂ_maxとする。即ち、式（１８）が成り立つ場合を条件とする。

測度ｂ_iの最大測度ｂ_maxに対する割合を、次の式（１９）に示す。

ｑ（ｉ），ｉ＝１，２，・・・，ｎに関する情報量は次の式（２０）である。

エントロピーは情報量［−ｌｏｇｐ（ｉ）］の平均（式（１６））であるのに対して、修正エントロピーｓ´を、情報量［−ｌｏｇｐ（ｉ）］から情報量［−ｌｏｇｑ（ｉ）］を差し引いた情報量の平均と定義する。即ち、修正エントロピーｓ´は、次の式（２１）、（２２）により示される。

＜条件付確率ｐ（ｊ｜ｉ）に関するエントロピー＞
エントロピーｓを式（１６）に示すように、確率ｐ（ｉ）に関するエントロピーと考えたとき、条件付確率ｐ（ｊ｜ｉ）に関するエントロピーｓ（ｉ）を以下の式（２３）、（２４）のように定義することができる。

＜条件付確率ｐ（ｊ｜ｉ）に関する修正エントロピー＞
上記の式（１８）に記述した最大測度ｂ_maxに関して、同時測度ｂ_i,jのｂ_maxに対する割合をｑ（ｉ，ｊ）とすると、次の式（２５）により示すことができる。

ｑ（ｉ，ｊ），ｉ＝１，２，・・・，ｍ，ｊ＝１，２，・・・，ｍに関する情報量は次の式（２６）に示す通りである。

条件付確率ｐ（ｊ｜ｉ）に関する修正エントロピーｓ´（ｉ）を情報量［−ｌｏｇｐ（ｊ｜ｉ）］から情報量［−ｌｏｇｑ（ｉ，ｊ）］を差し引いた情報量の平均と定義する。即ち、修正エントロピーｓ´（ｉ）は、次の式（２７）により示される。

確率算出手段２０３は、上記式（２５）により確率ｑ（ｉ，ｊ）を求めることができる。

即ち、式（１３）、（１４）より、次の式（２８）が成り立つ。

＜再びカオス尺度＞
カオス尺度Ｈについては、＜カオス尺度の定義＞において、式（６）に示した通りであるが、このカオス尺度Ｈは、条件付確率ｐ（ｊ｜ｉ）に関するエントロピーｓ（ｉ）として式（２３）、（２４）に示したものの平均と解釈でき、次の式（２９）〜（３１）に示す通りとなる。

＜修正カオス尺度＞
本実施形態では、修正カオス尺度Ｈ´を、条件付確率ｐ（ｊ｜ｉ）に関する修正エントロピーｓ´（ｉ）として式（２７）、（２８）に示したものの平均と定義する。即ち、修正カオス尺度Ｈ´を、次の式（３２）、（３３）、（３４）から求める。

修正カオス尺度算出手段２０４は、条件付確率に関する修正エントロピーｓ´（ｉ）を、第１の情報量［−ｌｏｇｐ（ｊ｜ｉ）］から第２の情報量［−ｌｏｇｑ（ｉ，ｊ）］を引いた情報量の平均と定義するとき、前記条件付確率に関する修正エントロピーｓ´（ｉ）の平均として定義された修正カオス尺度を、前記確率分布算出手段２０１により算出された確率分布ｐ（ｉ）と、前記条件付確率算出手段により求められた条件付確率ｐ（ｊ｜ｉ）と、前記確率算出手段により求められた確率ｑ（ｉ，ｊ）を用いて求めるものである。即ち、修正カオス尺度算出手段２０４は、上記式（３２）、（３３）、（３４）から修正カオス尺度Ｈ´を求める。

＜Ｈ´の計算法１＞
＜Ｈ´の計算法１＞は、写像τと不変測度（確率密度関数）ρが既知である場合に適用される。この計算法１による時系列データ解析装置或いは時系列データ解析用プログラムによって行われる処理を、図６のフローチャートを参照して説明する。初期設定として、写像をτ、分割数をｍ、写像の定義域をＩ＝［ａ，ｂ］とする（Ｓ１１）。

図７には、ｍ分割した写像τの１区間||Ｘ_i||について、写像後の値がどのように変化するかを表す計算を行う計算工程説明図が示されている。ステップＳ１１に続いて、修正カオス尺度Ｈ´を求める演算を行うためのレジスタｈ、レジスタｉ、レジスタｊの値をリセットする（ｈは０，ｉ，ｊは１）（Ｓ１２）。

ステップＳ１３からステップＳ１８までにおいてはｉをｍまで１づつ歩進させて処理を行うステップである。ステップＳ１３においてｉを１歩進させると、ステップＳ１４では次の式（３５）によりｐ（ｉ）を求める。

次のステップＳ１５からステップＳ１７までにおいてはｊをｍまで１づつ歩進させて処理を行うステップである。ステップＳ１５においてｊを１歩進させると、ステップＳ１６では次の式（３６）によりｐ（ｊ｜ｉ）を求め、更に式（３７）によりｑ（ｉ，ｊ）を求める。

更に、ステップＳ１６においては、修正カオス尺度Ｈ´を得るための途中計算として次の式（３８）によるｕを求める計算と、式（３９）によるｈを求める計算を行う。このｕを求める計算と、ｈを求める計算は、修正カオス尺度算出手段２０４が行う。なお、式（３８）、（３９）においては、矢印「←」の右辺の値を左辺にセットすることを示す。

ステップＳ１６が終了すると、ステップＳ１７において、ｊ＝ｍとなったかを検出し、ｊ＝ｍとなっていなければステップＳ１５へ戻ってｊを１歩進させると、ステップＳ１６へと進む処理を繰り返す。このステップＳ１５からステップＳ１７の処理は、図７に示すｐ（ｊ｜ｉ）とｑ（ｉ，ｊ）との計算を示す図において、同じＸ_i内においてｊを歩進させて||Ｘ_j||のｊの値を歩進させて演算を行う処理であり、Ｘ_iは１つの区間に留まるものである。

このようにして処理を繰り返すうちにステップＳ１７を通過し、ステップＳ１８へ処理が進む。ステップＳ１８においては、ｉ＝ｍとなったかを検出し、ｉ＝ｍとなっていなければステップＳ１３へ戻ってｉを１歩進させると、ステップＳ１４へと進む処理を繰り返す。このステップＳ１３からステップＳ１８の処理は、図７に示すｐ（ｊ｜ｉ）とｑ（ｉ，ｊ）との計算を示す図において、||Ｘ_i||のｉの値を歩進させて演算を行う処理であって、Ｘ_iは現在の１区間から次の１区間へ移動して処理が行われる。１つＸ_iが区間を移動すると、移動した１区間においてステップＳ１５からステップＳ１７の処理が行われる。

上記のようにしてステップＳ１３からステップＳ１８の処理を繰り返すうちにステップＳ１８を通過すると、ステップＳ１９においては、既に求まっているｈを修正カオス尺度Ｈ´として結果出力部１０６から出力する。

＜Ｈ´の計算法２＞
＜Ｈ´の計算法２＞は、写像τと不変測度（確率密度関数）ρが未知であり、与えられた時系列データから修正カオス尺度Ｈ´を計算する場合に適用される。＜Ｈ´の計算法２＞では、写像τと不変測度（確率密度関数）ρが未知であることが前提であるから＜Ｈ´の計算法１＞において使用した式（３５）〜式（３７）ではなく、ｐ（ｉ）、ｐ（ｊ｜ｉ）、ｑ（ｉ，ｊ）を、与えられた時系列データから求めることになる。即ち、確率分布算出手段２０１は、データ長がｎ＋１の時系列データ｛ξ₀，ξ₁，ξ₂，・・・，ξ_n｝から、確率分布ｐ（ｉ）、条件付確率ｐ（ｊ｜ｉ）、確率ｑ（ｉ，ｊ）を求める。

この計算法２による時系列データ解析装置或いは時系列データ解析用プログラムによって行われる処理を、図８、図９のフローチャートを参照して説明する。初期設定として、分割数をｍ、細分割数をＭ、時系列データの定義域をＩ＝［ａ，ｂ］とし、Δｘ＝（ｂ−ａ）／ｍとする（Ｓ２１）。更に、レジスタｈとレジスタｉ，ｊ，ｋの初期化を行い（Ｓ２２）、カウンタｃ₁［］，ｃ₂［］［］，ｃ₃［］［］の初期化が行われる（Ｓ２３）。次に、データ長がｎ＋１の時系列データ｛ξ₀，ξ₁，ξ₂，・・・，ξ_n｝の取り込みが行われる（Ｓ２４）。

ステップＳ２５からステップＳ２７までにおいてはｋをｎ−１まで１づつ歩進させて処理を行うステップである。ステップＳ２５においてｋを１歩進させ、ステップＳ２６ではｐ（ｉ）、ｐ（ｊ｜ｉ）、ｑ（ｉ，ｊ）を求めるために、カウンタｃ₁［］，ｃ₂［］［］，ｃ₃［］［］の値を確定する。

具体的には、次の式（４０）、（４１）、（４２）からν₁、ν₂、ν₃を求める。

上記の通りの計算からν₁、ν₂、ν₃が求められると、次の式（４３）、（４４）、（４５）からカウンタｃ₁［］，ｃ₂［］［］，ｃ₃［］［］の値が求められる。

上記ステップＳ２６の処理が終了すると、ステップＳ２７においてｋがｎ−１となったかを検出し、ｎ−１となっていなければ、ステップＳ２５へ戻りｋを１歩進してステップＳ２６の処理を行い、ステップＳ２７へ進みｋがｎ−１となったかを検出する処理を繰り返す。このようにして、処理を繰り返すとステップＳ２７においてｋ＝ｎ−１の成立が検出されると、図９のフローチャートに示す処理が行われる。即ち、ステップＳ３１からステップＳ３６までにおいてはｉをｍまで１づつ歩進させて処理を行うステップである。

ステップＳ３１においてはｉをｍまで１づつ歩進させると、ｐ（ｉ）へｃ₁［ｉ］／ｎをセットし、ステップＳ３３へ進む。ステップＳ３３からステップＳ３５まではｊをｍまで１づつ歩進させて処理を行うステップである。ステップＳ３３においてｊを１歩進させ、ｐ（ｉ，ｊ）を求め、ｐ（ｊ｜ｉ）を求める。ステップＳ２６における処理によってｃ₁［ｉ］，ｃ₂［ｉ］［ｊ］は、それぞれ次の式（４６）、（４７）によって表されるように求まっている。

上記からｐ（ｉ）は既に述べた通り、ｃ₁［ｉ］／ｎがセットされた式（４８）により、ｐ（ｉ，ｊ）、ｐ（ｊ｜ｉ）は、それぞれ式（４９）、（５０）により表される式により計算される。

このステップＳ３４では、時系列データからｑ（ｉ，ｊ）についても求める。ｑ（ｉ，ｊ）を求める意味は、分割区間Ｘ_i内にデータが存在する範囲の割合を求めることである。そこで、分割区間Ｘ_iを更にＭ等分した細分割区間Ｙ_iを用いる。即ち、区間Ｉをｍ×Ｍ等分した細分割区間をＹ_i（ｉ＝１，２，・・・，ｍ×Ｍ）とする。Ｉは、次の式（５１）により表すことができる。

また、上記細分割区間をＹ_iと分割区間Ｘ_iの関係は以下の式（５２）の通りである。

具体的な手法を以下に述べる。既にステップＳ２３において説明した通り、カウンタｃ₃［ｉ］［ｊ］を用意し、初期化すると次の式（５３）のようになる。

次に、ξ_k∈Ｘ_iの基で、ξ_k+1∈Ｙ_jとなる測度を次の式（５４）により示す。

上記の式（５４）によりカウンタを更新する（ｃ₃［ｉ］［ｊ］の最大値が１である点に注意する。）。従って、

と、上記式（５３）により定義されるｔ（ｉ，ｊ）は、ξ_k∈Ｘ_iの基で、ξ_k+1∈Ｘ_jを満たすＸ_jにおいて、Ｘ_jに内にデータが存在する範囲の割合を意味する。なお、ｔ（ｉ，ｊ）は、式（３４）における第２項のｂ_i,jに相応する。ｔ（ｉ，ｊ）の最大値がＭであることから、Ｍは式（３４）における第２項のｂ_i,jに相応する。以上により、ｑ（ｉ，ｊ）は、上記ｔ（ｉ，ｊ）を用いて次の式（５６）に示されるように記載することができる。

以上によって式（４９）、（５０）、（５６）を用いて得たｐ（ｉ，ｊ）、ｐ（ｊ｜ｉ）、ｑ（ｉ，ｊ）を式（３３）へ代入し修正カオス尺度Ｈ´を算出し、修正カオス尺度Ｈ´として結果出力部１０６から出力することができる。

以上のようにして＜Ｈ´の計算法２＞によりロジスティック写像から得られた修正カオス尺度と、従来のカオス尺度、更にリアプノフ指数の波形を図１０、図１１に示す。図１０では、分割数ｍが２０、細分割数Ｍが２０、データ長を定義するｎが１００００である場合を示す。また、図１１では、分割数ｍが１００、細分割数Ｍが１００、データ長を定義するｎが１００００００である場合を示す。いずれの場合も修正カオス尺度がリアプノフ指数に近似しており、時系列データ解析に極めて優れたものであることを示している。

＜Ｈ´の計算法３＞
本実施形態に係るデータ解析手法は、処理が所謂軽量にできているため、リアルタイムで得られるデータに対してリアプノフ指数による解析に極めて近似した解析を行うことが可能である。時系列データがリアルタイムに得られている場合に、修正カオス尺度Ｈ´を得る計算法を＜Ｈ´の計算法３＞と称し、以下に説明を行う。図１２のタイミングチャートに示されるように、時系列データは最上部の右横向き矢印により示されている通りリアルタイムで常時到来するものとする。

初めに時刻ｔ₀から時刻ｔ₁までの時間間隔Ｔ（＝ｔ₁−ｔ₀）における時系列データＤＡ１を抽出して第１番目の抽出データとする。この第１番目の抽出データは時刻ｔ₁の直後に本実施形態に係る時系列データ解析用プログラムに渡され、処理が行われて第１番目の修正カオス尺度Ｈ₁´が得られ、その後に結果出力部１０６から出力される。

次に、時刻ｔ₀＋Δｔから時刻ｔ₁＋Δｔまでの時間間隔Ｔ（＝ｔ₂−ｔ₁）における時系列データＤＡ２を抽出して第２番目の抽出データとする。この第２番目の抽出データは時刻ｔ₁＋Δｔの直後に本実施形態に係る時系列データ解析用プログラムに渡され、処理が行われて第２番目の修正カオス尺度Ｈ₂´が得られ、その後に結果出力部１０６から出力される。

次に、時刻ｔ₀＋２Δｔから時刻ｔ₁＋２Δｔまでの時間間隔Ｔ（＝ｔ₃−ｔ₂）における時系列データＤＡ３を抽出して第３番目の抽出データとする。この第３番目の抽出データは時刻ｔ₁＋３Δｔの直後に本実施形態に係る時系列データ解析用プログラムに渡され、処理が行われて第３番目の修正カオス尺度Ｈ₃´が得られ、その後に結果出力部１０６から出力される。

更に、同様に第４番目の抽出データから第４番目の修正カオス尺度Ｈ₄´が得られ、以下同様に処理が続けられる。即ち、本実施形態では、データ長がｎ＋１の時系列データ｛ξ₀，ξ₁，ξ₂，・・・，ξ_n｝をサンプリング期間をずらして複数回サンプリングし、上記修正カオス尺度算出手段２０４は、サンプリングした複数回分の時系列データ毎に修正カオス尺度算出を行う。斯くして、本実施形態では、次々に所要時間修正カオス尺度Ｈ₄´が得られ、それぞれの時間単位でデータ解析を行うことができる。

なお、上記実施の形態では、ｋ番目の抽出データとｋ＋１番目の抽出データとがオーバーラップする部分があるが、ｋ番目とｋ＋１番目の間の抽出時間間隔Δｔにおいて大量のデータが到来するような場合には、オーバーラップ部分を設けることなく計算を行っても良いことは勿論である。

サンプリング期間をずらして複数回サンプリングして、サンプリングした複数回分の時系列データ毎に修正カオス尺度算出を行う（第１の処理態様という）のではなく、修正カオス尺度算出手段２０４は、データ長がｎ＋１の時系列データ｛ξ₀，ξ₁，ξ₂，・・・，ξ_n｝を１度得て修正カオス尺度算出を行うよう（第２の処理態様という）にしても良い。また、第１の処理態様と第２の処理態様と双方を有し、データの長さや計算を行う時間帯に応じて、更には、オペレータの指示により、上記の２態様の一方を採用するようにしても良い。

本願の発明者らは、非特許文献１２において、心拍間隔データのカオス尺度と自律神経活動の関連について記述し、非特許文献１３において、心拍変動のカオス性と副交感神経活動との関係について記述した。このように心拍変動のカオス性と眠気や怠慢等の生理状態とは相関が見られることから、本実施形態の解析装置や解析用プログラムは、生体の生理状態について解析する場合に極めて有用な修正カオス尺度を時系列データから得ることができるものである。

電子機器等においては、電源電圧の不規則な揺らぎ（カオス的な変化）は、電圧供給先である内部の電子部品乃至素子等の抵抗値等に係る電気的特性が不規則に変動している可能性を示唆するものである。従って、電子機器等の電源電圧を時系列データとして本実施形態に係る解析装置や解析用プログラムを用いて解析を行うことによって電子部品乃至素子等の異常や故障の検出並びに予測へつながることが期待される。

更に、機械学習やディープラーニングは、入力データの特徴Ａに対し人間や自然界或いは機械・装置において生じる事象Ｂとの対応関係を過去のデータを用いて学習しておき、新たな入力データｎｅｗＡに対して最も生じる可能性が高い事象ＢＢを自動的に割り出すものと言える。そこで、本実施形態に係る解析装置や解析用プログラムを用いて入力データの特徴Ａを解析しておき、これに対する事象Ｂを対応付けて、機械学習やディープラーニングへ与えて学習させ、その後に到来する入力データについての、本実施形態に係る解析装置や解析用プログラムを用いて得た特徴から最も生じる可能性が高い事象の推定などに適用でき得るものと考えられる。

まとめるならば、機械学習を含む一般的なデータ解析に用いられる数学的手法は、データの規則性を見出すことを主目的として発展してきた背景がある。これに対し、本実施形態に係る解析装置や解析用プログラムのように規則性とは逆の事象である不規則性（カオス性）を把握して事象に対する特徴とする手法は一般的に行われているものではない。従って、本実施形態に係る解析装置や解析用プログラムを用いて得た不規則性の特徴を用いるものは、発生事象の予測や発生している事象の検出など幅広い応用が期待できる。

１００ 時系列データ解析装置
１０２ 主メモリ
１０３ バス
１０４ 外部記憶装置
１０５ 時系列データ供給部
１０６ 結果出力部
２０１ 確率分布算出手段
２０２ 条件付確率算出手段
２０３ 確率算出手段
２０４ 修正カオス尺度算出手段

Claims (12)

1. ｎを正の整数として、データ長がｎ＋１の時系列データ｛ξ₀，ξ₁，ξ₂，・・・，ξ_n｝が、写像τによりξ_k＝τ（ξ_k-1）＝τ^k（ξ₀），ｋ＝１，２，・・・，ｎとして生成されるという仮定の基に、ξ_kが含まれる区間Ｉ≡［ａ，ｂ］をｍ個の区間に等分割した分割区間Ｘ_i（ｉ＝１，２，・・・，ｍ）とについて、
   

   

   であり、ξ_k∈Ｘ_iとなる確率分布ｐ（ｉ）を求め、ξ_k∈Ｘ_i、ξ_k+1∈Ｘ_jとなる同時確率分布ｐ（ｉ，ｊ）を求める確率分布算出手段と、
   条件付確率ｐ（ｊ｜ｉ）を求める条件付確率算出手段と、
   ｍ個の測度ｂ_jの各クラスが｛ｂ₁，ｂ₂，ｂ₃，・・・，ｂ_m｝であり、各クラスの最大値がｂ_max、であるとき、測度ｂ_iの最大値ｂ_maxに対する割合をｑ（ｉ）としたとき、同時測度ｂ_i,jについては、最大値ｂ_maxに対する割合である確率ｑ（ｉ，ｊ）を求める確率算出手段と、
   条件付確率に関する修正エントロピーｓ´（ｉ）を、第１の情報量［−ｌｏｇｐ（ｊ｜ｉ）］から第２の情報量［−ｌｏｇｑ（ｉ，ｊ）］を引いた情報量の平均と定義するとき、前記条件付確率に関する修正エントロピーｓ´（ｉ）の平均として定義された修正カオス尺度を、前記確率分布算出手段により算出された確率分布ｐ（ｉ）と、前記条件付確率算出手段により求められた条件付確率ｐ（ｊ｜ｉ）と、前記確率算出手段により求められた確率ｑ（ｉ，ｊ）を用いて求める修正カオス尺度算出手段と
   を備えることを特徴とする時系列データ解析装置。
2. 前記条件付確率算出手段は、前記確率分布算出手段が求めた確率分布ｐ（ｉ）と同時確率分布ｐ（ｉ，ｊ）から条件付確率ｐ（ｊ｜ｉ）を求めることを特徴とする請求項１に記載の時系列データ解析装置。
3. 写像τと、不変測度（確率密度関数）ρが既知である場合に、前記確率分布算出手段は、確率分布ｐ（ｉ）を
   

   

   により、
   前記条件付確率算出手段は、条件付確率ｐ（ｊ｜ｉ）を、
   

   

   により、
   前記確率分布算出手段は、確率ｑ（ｉ，ｊ）を、
   

   

   により、算出することを特徴とする請求項１に記載の時系列データ解析装置。
4. 前記確率分布算出手段は、
   データ長がｎ＋１の時系列データ｛ξ₀，ξ₁，ξ₂，・・・，ξ_n｝から、確率分布ｐ（ｉ）、確率ｑ（ｉ，ｊ）を求めることを特徴とする請求項１または２に記載の時系列データ解析装置。
5. データ長がｎ＋１の時系列データ｛ξ₀，ξ₁，ξ₂，・・・，ξ_n｝を１度得て、前記修正カオス尺度算出手段は、修正カオス尺度算出を行うことを特徴とする請求項１乃至４のいずれか１項に記載の時系列データ解析装置。
6. データ長がｎ＋１の時系列データ｛ξ₀，ξ₁，ξ₂，・・・，ξ_n｝をサンプリング期間をずらして複数回サンプリングし、前記修正カオス尺度算出手段は、サンプリングした複数回分の時系列データ毎に修正カオス尺度算出を行うことを特徴とする請求項１乃至４のいずれか１項に記載の時系列データ解析装置。
7. コンピュータを
   ｎを正の整数として、データ長がｎ＋１の時系列データ｛ξ₀，ξ₁，ξ₂，・・・，ξ_n｝が、写像τによりξ_k＝τ（ξ_k-1）＝τ^k（ξ₀），ｋ＝１，２，・・・，ｎとして生成されるという仮定の基に、ξ_kが含まれる区間Ｉ≡［ａ，ｂ］をｍ個の区間に等分割した分割区間Ｘ_i（ｉ＝１，２，・・・，ｍ）とについて、
   

   

   であり、ξ_k∈Ｘ_iとなる確率分布ｐ（ｉ）を求め、ξ_k∈Ｘ_i、ξ_k+1∈Ｘ_jとなる同時確率分布ｐ（ｉ，ｊ）を求める確率分布算出手段、
   条件付確率ｐ（ｊ｜ｉ）を求める条件付確率算出手段、
   ｍ個の測度ｂ_jの各クラスが｛ｂ₁，ｂ₂，ｂ₃，・・・，ｂ_m｝であり、各クラスの最大値がｂ_max、であるとき、測度ｂ_iの最大値ｂ_maxに対する割合をｑ（ｉ）としたとき、同時測度ｂ_i,jについては、最大値ｂ_maxに対する割合である確率ｑ（ｉ，ｊ）を求める確率算出手段、
   条件付確率に関する修正エントロピーｓ´（ｉ）を、第１の情報量［−ｌｏｇｐ（ｊ｜ｉ）］から第２の情報量［−ｌｏｇｑ（ｉ，ｊ）］を引いた情報量の平均と定義するとき、前記条件付確率に関する修正エントロピーｓ´（ｉ）の平均として定義された修正カオス尺度を、前記確率分布算出手段により算出された確率分布ｐ（ｉ）と、前記条件付確率算出手段により求められた条件付確率ｐ（ｊ｜ｉ）と、前記確率算出手段により求められた確率ｑ（ｉ，ｊ）を用いて求める修正カオス尺度算出手段
   として機能させることを特徴とする時系列データ解析用プログラム。
8. 前記コンピュータを前記条件付確率算出手段として、
   前記コンピュータが前記確率分布算出手段として求めた確率分布ｐ（ｉ）と同時確率分布ｐ（ｉ，ｊ）から条件付確率ｐ（ｊ｜ｉ）を求めるように機能させることを特徴とする請求項７に記載の時系列データ解析用プログラム。
9. 写像τと、不変測度（確率密度関数）ρが既知である場合に、前記コンピュータを前記確率分布算出手段として、確率分布ｐ（ｉ）を
   

   

   により、
   前記コンピュータを前記条件付確率算出手段として、条件付確率ｐ（ｊ｜ｉ）を、
   

   

   により、
   前記コンピュータを前記確率分布算出手段として、確率ｑ（ｉ，ｊ）を、
   

   

   により、算出するように機能させることを特徴とする請求項７に記載の時系列データ解析用プログラム。
10. 前記コンピュータを前記確率分布算出手段として、
    データ長がｎ＋１の時系列データ｛ξ₀，ξ₁，ξ₂，・・・，ξ_n｝から、確率分布ｐ（ｉ）、確率ｑ（ｉ，ｊ）を求めるように機能させることを特徴とする請求項７または８に記載の時系列データ解析用プログラム。
11. データ長がｎ＋１の時系列データ｛ξ₀，ξ₁，ξ₂，・・・，ξ_n｝を１度得て、前記コンピュータを前記修正カオス尺度算出手段として、修正カオス尺度算出を行うように機能させることを特徴とする請求項７乃至１０のいずれか１項に記載の時系列データ解析用プログラム。
12. データ長がｎ＋１の時系列データ｛ξ₀，ξ₁，ξ₂，・・・，ξ_n｝をサンプリング期間をずらして複数回サンプリングし、
    前記コンピュータを前記修正カオス尺度算出手段として、サンプリングした複数回分の時系列データ毎に修正カオス尺度算出を行うように機能させることを特徴とする請求項７乃至１０のいずれか１項に記載の時系列データ解析用プログラム。

Images