(57)【特許請求の範囲】
【請求項１】
ｎを正の整数として、データ長がｎ＋１の時系列データ｛ξ0，ξ1，ξ2，・・・，ξn｝が、写像τによりξk＝τ（ξk-1）＝τk（ξ0），ｋ＝１，２，・・・，ｎとして生成されるという仮定の基に、ξkが含まれる区間Ｉ≡［ａ，ｂ］をｍ個の区間に等分割した分割区間Ｘi（ｉ＝１，２，・・・，ｍ）とについて、
【数１】
0006869514.tif 000046
であり、ξk∈Ｘiとなる確率分布ｐ（ｉ）を求め、ξk∈Ｘi、ξk+1∈Ｘjとなる同時確率分布ｐ（ｉ，ｊ）を求める確率分布算出手段と、
条件付確率ｐ（ｊ｜ｉ）を求める条件付確率算出手段と、
ｍ個の測度ｂjの各クラスが｛ｂ1，ｂ2，ｂ3，・・・，ｂm｝であり、各クラスの最大値がｂmax、であるとき、測度ｂiの最大値ｂmaxに対する割合をｑ（ｉ）としたとき、同時測度ｂi,jについては、最大値ｂmaxに対する割合である確率ｑ（ｉ，ｊ）を求める確率算出手段と、
条件付確率に関する修正エントロピーｓ´（ｉ）を、第１の情報量［−ｌｏｇｐ（ｊ｜ｉ）］から第２の情報量［−ｌｏｇｑ（ｉ，ｊ）］を引いた情報量の平均と定義するとき、前記条件付確率に関する修正エントロピーｓ´（ｉ）の平均として定義された修正カオス尺度を、前記確率分布算出手段により算出された確率分布ｐ（ｉ）と、前記条件付確率算出手段により求められた条件付確率ｐ（ｊ｜ｉ）と、前記確率算出手段により求められた確率ｑ（ｉ，ｊ）を用いて求める修正カオス尺度算出手段と
を備えることを特徴とする時系列データ解析装置。
【請求項２】
前記条件付確率算出手段は、前記確率分布算出手段が求めた確率分布ｐ（ｉ）と同時確率分布ｐ（ｉ，ｊ）から条件付確率ｐ（ｊ｜ｉ）を求めることを特徴とする請求項１に記載の時系列データ解析装置。
【請求項３】
写像τと、不変測度（確率密度関数）ρが既知である場合に、前記確率分布算出手段は、確率分布ｐ（ｉ）を
【数２】
0006869514.tif 000047
により、
前記条件付確率算出手段は、条件付確率ｐ（ｊ｜ｉ）を、
【数３】
0006869514.tif 000048
により、
前記確率分布算出手段は、確率ｑ（ｉ，ｊ）を、
【数４】
0006869514.tif 000049
により、算出することを特徴とする請求項１に記載の時系列データ解析装置。
【請求項４】
前記確率分布算出手段は、
データ長がｎ＋１の時系列データ｛ξ0，ξ1，ξ2，・・・，ξn｝から、確率分布ｐ（ｉ）、確率ｑ（ｉ，ｊ）を求めることを特徴とする請求項１または２に記載の時系列データ解析装置。
【請求項５】
データ長がｎ＋１の時系列データ｛ξ0，ξ1，ξ2，・・・，ξn｝を１度得て、前記修正カオス尺度算出手段は、修正カオス尺度算出を行うことを特徴とする請求項１乃至４のいずれか１項に記載の時系列データ解析装置。
【請求項６】
データ長がｎ＋１の時系列データ｛ξ0，ξ1，ξ2，・・・，ξn｝をサンプリング期間をずらして複数回サンプリングし、前記修正カオス尺度算出手段は、サンプリングした複数回分の時系列データ毎に修正カオス尺度算出を行うことを特徴とする請求項１乃至４のいずれか１項に記載の時系列データ解析装置。
【請求項７】
コンピュータを
ｎを正の整数として、データ長がｎ＋１の時系列データ｛ξ0，ξ1，ξ2，・・・，ξn｝が、写像τによりξk＝τ（ξk-1）＝τk（ξ0），ｋ＝１，２，・・・，ｎとして生成されるという仮定の基に、ξkが含まれる区間Ｉ≡［ａ，ｂ］をｍ個の区間に等分割した分割区間Ｘi（ｉ＝１，２，・・・，ｍ）とについて、
【数５】
0006869514.tif 000050
であり、ξk∈Ｘiとなる確率分布ｐ（ｉ）を求め、ξk∈Ｘi、ξk+1∈Ｘjとなる同時確率分布ｐ（ｉ，ｊ）を求める確率分布算出手段、
条件付確率ｐ（ｊ｜ｉ）を求める条件付確率算出手段、
ｍ個の測度ｂjの各クラスが｛ｂ1，ｂ2，ｂ3，・・・，ｂm｝であり、各クラスの最大値がｂmax、であるとき、測度ｂiの最大値ｂmaxに対する割合をｑ（ｉ）としたとき、同時測度ｂi,jについては、最大値ｂmaxに対する割合である確率ｑ（ｉ，ｊ）を求める確率算出手段、
条件付確率に関する修正エントロピーｓ´（ｉ）を、第１の情報量［−ｌｏｇｐ（ｊ｜ｉ）］から第２の情報量［−ｌｏｇｑ（ｉ，ｊ）］を引いた情報量の平均と定義するとき、前記条件付確率に関する修正エントロピーｓ´（ｉ）の平均として定義された修正カオス尺度を、前記確率分布算出手段により算出された確率分布ｐ（ｉ）と、前記条件付確率算出手段により求められた条件付確率ｐ（ｊ｜ｉ）と、前記確率算出手段により求められた確率ｑ（ｉ，ｊ）を用いて求める修正カオス尺度算出手段
として機能させることを特徴とする時系列データ解析用プログラム。
【請求項８】
前記コンピュータを前記条件付確率算出手段として、
前記コンピュータが前記確率分布算出手段として求めた確率分布ｐ（ｉ）と同時確率分布ｐ（ｉ，ｊ）から条件付確率ｐ（ｊ｜ｉ）を求めるように機能させることを特徴とする請求項７に記載の時系列データ解析用プログラム。
【請求項９】
写像τと、不変測度（確率密度関数）ρが既知である場合に、前記コンピュータを前記確率分布算出手段として、確率分布ｐ（ｉ）を
【数６】
0006869514.tif 000051
により、
前記コンピュータを前記条件付確率算出手段として、条件付確率ｐ（ｊ｜ｉ）を、
【数７】
0006869514.tif 000052
により、
前記コンピュータを前記確率分布算出手段として、確率ｑ（ｉ，ｊ）を、
【数８】
0006869514.tif 000053
により、算出するように機能させることを特徴とする請求項７に記載の時系列データ解析用プログラム。
【請求項１０】
前記コンピュータを前記確率分布算出手段として、
データ長がｎ＋１の時系列データ｛ξ0，ξ1，ξ2，・・・，ξn｝から、確率分布ｐ（ｉ）、確率ｑ（ｉ，ｊ）を求めるように機能させることを特徴とする請求項７または８に記載の時系列データ解析用プログラム。
【請求項１１】
データ長がｎ＋１の時系列データ｛ξ0，ξ1，ξ2，・・・，ξn｝を１度得て、前記コンピュータを前記修正カオス尺度算出手段として、修正カオス尺度算出を行うように機能させることを特徴とする請求項７乃至１０のいずれか１項に記載の時系列データ解析用プログラム。
【請求項１２】
データ長がｎ＋１の時系列データ｛ξ0，ξ1，ξ2，・・・，ξn｝をサンプリング期間をずらして複数回サンプリングし、
前記コンピュータを前記修正カオス尺度算出手段として、サンプリングした複数回分の時系列データ毎に修正カオス尺度算出を行うように機能させることを特徴とする請求項７乃至１０のいずれか１項に記載の時系列データ解析用プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、時系列データ解析装置及び時系列データ解析用プログラムに関するものである。
【背景技術】
【０００２】
時系列データがどのような特徴を有しているのかを把握するために、様々な手法が用いられる。例えば、時系列データについてカオス度合を把握するための定量化法として、リアプノフ指数と呼ばれる指標がよく用いられている。リアプノフ指数を求める場合は、データ生成源の情報（離散系の差分方程式や連続系の微分方程式等）が既知である必要がある。データ生成源の情報が未知である場合には、大量のデータから推定する手段が与えられているが、所謂「埋め込み次元の推定」に関する処理は、次元毎に系の様子を見ながら調整を繰り返す必要があり煩雑であるにもかかわらず、必ずしも一意に得られるものではないという問題がある（非特許文献５、６、７）。
【０００３】
近年になって情報理論に基づく「カオス尺度（Chaos Degree）」と称される指標が提案された（非特許文献１）。このカオス尺度とリアプノフ指数は極めて類似の挙動を示すことが知られている（非特許文献２、３、４）。図１には、ロジスティック写像（ａ：３．５〜４．０）のカオス尺度（分割数ｍ＝１００）とリアプノフ指数の値の変化を示している。この図１から、カオス尺度とリアプノフ指数カオス尺度は、類似の変化を行うことが明らかである。カオス尺度は、データのみから一意に計算可能であるという特徴を有している。
【０００４】
近年になって、本願発明者らは、カオス尺度とリアプノフ指数の数理学的関係性を明らかにした（非特許文献８、９、１０）。
【０００５】
＜カオス尺度の定義＞
データ長がｎ＋１の時系列データを
｛ξ0，ξ1，ξ2，・・・，ξn｝・・・（１）
と表す。
データは、差分方程式τ
【数９】
0006869514.tif 000002
により、
【数１０】
0006869514.tif 000003
として生成されているものとする。
【０００６】
上記のξkが含まれている区間Ｉをｍ個の区間に等分割する。この分割区間をＸi（ｉ＝１，２，・・・，ｍ）とすると、Ｉについて次の式（２）が成り立つ。
【数１１】
0006869514.tif 000004
【０００７】
ここで、ξk∈Ｘiとなる確率分布ｐ（ｉ）と、ξk∈Ｘi，ξk+1∈Ｘjとなる同時確率分布ｐ（ｉ，ｊ）を次の式（３）、（４）の通りに算出する。
【数１２】
0006869514.tif 000005
上記において＃［（条件式）］は、（条件式）を満たす数を意味する。
【０００８】
上記のように確率分布ｐ（ｉ）と同時確率分布ｐ（ｉ，ｊ）が求まると、カオス尺度Ｈは以下の式（５）または式（６）として定義される。
【数１３】
0006869514.tif 000006
ただし、０ｌｏｇ０＝０とする。
【０００９】
上記のカオス尺度Ｈの最小値は０であるが、最大値は分割数ｍに依存してｌｏｇｍである。
【数１４】
0006869514.tif 000007
カオス尺度Ｈの値は大きいほどカオス性が大きい（複雑性が高い）ことを意味する。即ち、カオス尺度Ｈに関する解釈は、次のようにまとめることができる。
【数１５】
0006869514.tif 000008
【００１０】
上記のカオス尺度Ｈを、（１）単調増加データ、（２）周期データ、（３）ノイズが加わった周期データ、（４）正規分布に従うランダムデータ、（５）一様ランダムデータ、のそれぞれについて計算した結果、図２に示すようになった。
【００１１】
上記カオス尺度Ｈの計算は次の通りに行った。即ち、各データのデータ長をｎ＋１＝１０１、カオス尺度の分割数をｍ＝５、データの値域［０，１］をｍ＝５で等分割した区間を分割区間として計算を行った。この場合のカオス尺度の最大値はｌｏｇｍ＝ｌｏｇ５＝１．６０９５であった。よって（５）の一様ランダムデータのカオス尺度が最大値に近く、カオス性（複雑度）が高いことを確認できる。
【００１２】
図２には、上記（１）〜（５）のデータについてカオス尺度Ｈと共に分散ｓ2を計算して示してある。この図２からは、直感的に、（１）と（２）のデータが複雑とは言えず、（３）のデータがやや複雑さを帯びており、（４）と（５）のデータがランダムデータ故に複雑であると結論できる。この図２では、（１）〜（５）へ移行するに従ってデータの複雑さが増すように並べており、カオス尺度Ｈの値が徐々に大きくなっており直感的に複雑度が増していることが分かる。
【００１３】
一方、分散ｓ2は、以下の式で表わされる。
【数１６】
0006869514.tif 000009
【００１４】
上記分散ｓ2は、暫し、「データのばらつき度合を示す」と説明される。分散は、データの平均値からの距離の二乗平均であって、データの乱雑さとは異なるものであることが図２に示されている値からも理解できる。即ち、データのデタラメさ、ランダムさ、乱雑さ、複雑さ、カオス性などを定量化して議論する場合には、分散を用いるよりもカオス尺度を用いる方が適切であることが分かる。
【００１５】
＜カオス尺度とリアプノフ指数＞
リアプノフ指数λは、不変測度（確率密度関数）ρを用いて次のように表わされる。
【数１７】
0006869514.tif 000010
【００１６】
本願の発明者らは、カオス尺度の分割をΔｘ≡｜Ｘｉ｜（ｊ＝１，２，・・・，ｍ）と表したとき、カオス尺度はΔｘ→０の極限において、
【数１８】
0006869514.tif 000011
の関係にあることを示した（非特許文献９、１１）。ここで、Ｄ（ｘ）は非負の関数である。よって、ＨΔx→0≧λの関係にある。このＤ関数の詳細は非特許文献８、９に示す通りである。｜τ´（ｘ）｜におけるＤ関数の上限（Ｄ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
【発明の概要】
【発明が解決しようとする課題】
【００２１】
本発明は、解析対象のデータが離散系であろうと連続系のデータであろうと、適切に特にカオス性解析を行うことが可能な時系列データ解析装置を提供する。
【課題を解決するための手段】
【００２２】
実施形態に係る時系列データ解析装置は、ｎを正の整数として、データ長がｎ＋１の時系列データ｛ξ0，ξ1，ξ2，・・・，ξn｝が、写像τによりξk＝τ（ξk-1）＝τk（ξ0），ｋ＝１，２，・・・，ｎとして生成されるという仮定の基に、ξkが含まれる区間Ｉ≡［ａ，ｂ］をｍ個の区間に等分割した分割区間Ｘi（ｉ＝１，２，・・・，ｍ）とについて、
【数１９】
0006869514.tif 000012
であり、ξk∈Ｘiとなる確率分布ｐ（ｉ）を求め、ξk∈Ｘi、ξk+1∈Ｘjとなる同時確率分布ｐ（ｉ，ｊ）を求める確率分布算出手段と、条件付確率ｐ（ｊ｜ｉ）を求める条件付確率算出手段と、ｍ個の測度ｂjの各クラスが｛ｂ1，ｂ2，ｂ3，・・・，ｂm｝であり、各クラスの最大値がｂmax、であるとき、測度ｂiの最大値ｂmaxに対する割合をｑ（ｉ）としたとき、同時測度ｂi,jについては、最大値ｂmaxに対する割合である確率ｑ（ｉ，ｊ）を求める確率算出手段と、条件付確率に関する修正エントロピーｓ´（ｉ）を、第１の情報量［−ｌｏｇｐ（ｊ｜ｉ）］から第２の情報量［−ｌｏｇｑ（ｉ，ｊ）］を引いた情報量の平均と定義するとき、前記条件付確率に関する修正エントロピーｓ´（ｉ）の平均として定義された修正カオス尺度を、前記確率分布算出手段により算出された確率分布ｐ（ｉ）と、前記条件付確率算出手段により求められた条件付確率ｐ（ｊ｜ｉ）と、前記確率算出手段により求められた確率ｑ（ｉ，ｊ）を用いて求める修正カオス尺度算出手段とを備えることを特徴とする。
【図面の簡単な説明】
【００２３】
【図１】ロジスティック写像のカオス尺度とリアプノフ指数の値の変化を示す図。
【図２】カオス尺度Ｈを、（１）単調増加データ、（２）周期データ、（３）ノイズが加わった周期データ、（４）正規分布に従うランダムデータ、（５）一様ランダムデータ、のそれぞれについて計算した結果を示す図。
【図３】写像をτとしたとき、｜τ´（ｘ）｜におけるＤ関数の上限（Ｄsup）と下限（Ｄinf）を示す図。
【図４】本発明の実施形態に係る時系列データ解析装置のブロック図。
【図５】本発明の実施形態に係る時系列データ解析装置の内部機能を示すブロック図。
【図６】本発明の実施形態に係る時系列データ解析装置による計算法１を用いて修正カオス尺度を算出する動作のフローチャート。
【図７】ｍ分割した写像τの１区間||Ｘi||について、写像後の値がどのように変化するかを表す計算を行う計算工程の説明図。
【図８】本発明の実施形態に係る時系列データ解析装置による計算法２を用いて修正カオス尺度を算出する動作のフローチャート（前半）。
【図９】本発明の実施形態に係る時系列データ解析装置による計算法２を用いて修正カオス尺度を算出する動作のフローチャート（後半）。
【図１０】本発明の実施形態に係る時系列データ解析装置による計算法２によりロジスティック写像から得られた修正カオス尺度と、従来のカオス尺度、更にリアプノフ指数の波形を示す図。
【図１１】本発明の実施形態に係る時系列データ解析装置による計算法２によりロジスティック写像から得られた修正カオス尺度と、従来のカオス尺度、更にリアプノフ指数の波形を示す図。
【図１２】計算法３により修正カオス尺度を算出する過程の説明を行うためのタイミングチャート。
【発明を実施するための形態】
【００２４】
以下添付図面を参照して、本発明の実施形態に係る時系列データ解析装置及び時系列データ解析用プログラムを説明する。各図において同一の構成要素には、同一の符号を付して重複する説明を省略する。図４に実施形態に係る時系列データ解析装置１００のブロック図を示す。時系列データ解析装置１００は、クラウドコンピュータ、サーバコンピュータ、パーソナルコンピュータ、その他のコンピュータにより構成することができる。
【００２５】
時系列データ解析装置１００は、ＣＰＵ１０１が主メモリ１０２のプログラムやデータに基づき演算を行うものである。ＣＰＵ１０１には、バス１０３を介して外部記憶装置１０４が接続されており、外部記憶装置１０４には、時系列データ解析用プログラムが記憶されている。ＣＰＵ１０１が外部記憶装置１０４から時系列データ解析用プログラムを主メモリ１０２へ読み出してこのプログラムを実行することにより時系列データ解析装置として機能する。
【００２６】
バス１０３には、外部記憶装置１０４以外に時系列データ供給部１０５が接続されている。時系列データ供給部１０５は、外部のセンサなどからリアルタイムで時系列データを取り込み保持するものとすることができ、或いは、外部の何らかの装置などが収集した時系列データを取り込み保持したものとすることができる。更に、収集した時系列データを記憶した媒体がセットされることにより、時系列データを保持し供給可能となっている装置であっても良い。更に、上記の構成を全て備えたものであっても良い。いずれにしても、ＣＰＵ１０１が時系列データ解析用プログラムを実行して時系列データの解析を行う場合には、時系列データはこの時系列データ供給部１０５から供給される。
【００２７】
バス１０３には、結果出力部１０６が接続されている。結果出力部１０６は、表示装置やプリンタなど、時系列データ解析装置１００において処理した結果を出力する装置とすることができる。また、結果出力部１０６は、時系列データ解析装置１００において処理した結果を記憶する媒体でもよく、更に、回線などを介して処理の依頼者（クライアント）へ処理結果を送信などする装置であっても良い。
【００２８】
外部記憶装置１０４に記憶されている時系列データ解析用プログラムが実行されることにより、図５に示される各手段が実現される。即ち。時系列データ解析装置１００は、図５に示されるように、確率分布算出手段２０１、条件付確率算出手段２０２、確率算出手段２０３、修正カオス尺度算出手段２０４を具備している。
【００２９】
本実施形態の時系列データ解析装置１００では、以下に示す＜修正カオス尺度＞を得るものであり、これにより式（７）に示したＤ関数に相当する部分を小さく抑えるものであり、即ち、カオス尺度とリアプノフ指数の差を小さく抑え、リアプノフ指数の推定値として機能させるものである。これによって、ＫＳエントロピーの推定値としても機能させるようにする。
【００３０】
確率分布算出手段２０１は、ｎを正の整数として、データ長がｎ＋１の時系列データ｛ξ0，ξ1，ξ2，・・・，ξn｝が、写像τによりξk＝τ（ξk-1）＝τk（ξ0），ｋ＝１，２，・・・，ｎとして生成されるという仮定の基に、ξkが含まれる、区間Ｉ≡［ａ，ｂ］をｍ個の区間に等分割した分割区間Ｘi（ｉ＝１，２，・・・，ｍ）とについて、
【数２０】
0006869514.tif 000013
であり、ξk∈Ｘiとなる確率分布ｐ（ｉ）を求め、ξk∈Ｘi、ξk+1∈Ｘjとなる同時確率分布ｐ（ｉ，ｊ）を求めるものである。
【００３１】
既に、＜カオス尺度の定義＞において述べたように、確率分布ｐ（ｉ）は式（３）により求められる。ここでは、データ長がｎ＋１の上記時系列データ｛ξ0，ξ1，ξ2，・・・，ξn｝を軌道と称する。この軌道ξkについて、ξk∈Ｘi（ｉ＝１，２，・・・，ｍ）を満たす測度（式（３）における右辺の分子）をｂiで表す。
【数２１】
0006869514.tif 000014
【００３２】
上記で、全測度は、定義より次の通りである。
【数２２】
0006869514.tif 000015
よって、各クラスの生起確率ｐ（ｉ）は、次のように記載することができる。
【数２３】
0006869514.tif 000016
上記式（１０）は、既述の式（３）に等しい。確率分布算出手段２０１は、上記式（１０）に示した各クラスの生起確率ｐ（ｉ）を算出するものである。
【００３３】
＜カオス尺度の定義＞においては、式（４）に、ξk∈Ｘi、ξk+1∈Ｘjとなる同時確率分布ｐ（ｉ，ｊ）を示した。ここでは、ξk∈Ｘi（ｉ＝１，２，・・・，ｍ）かつ、ξk+1∈Ｘj（ｉ＝１，２，・・・，ｍ）を同時に満たす測度（式（４）の右辺の分子）を、同時測度と称して、ｂi,jとすると、以下の式（１１）により表すことができる。
【数２４】
0006869514.tif 000017
【００３４】
上記の同時測度ｂi,jは、測度ｂiのうち、次の時刻にｊ番目のクラスに移る測度であるから、次の式（１２）に示す関係にある。
【数２５】
0006869514.tif 000018
上記により、同時確率ｐ（ｉ，ｊ）及び条件付確率ｐ（ｊ｜ｉ）は、次の式（１３）、（１４）により与えることができる。
【数２６】
0006869514.tif 000019
上記の式（１３）により示される同時確率ｐ（ｉ，ｊ）は、確率分布算出手段２０１が求めるものであり、式（４）と等しいものである。このように上記条件付確率算出手段２０２は、上記確率分布算出手段２０１が求めた確率分布ｐ（ｉ）と同時確率分布ｐ（ｉ，ｊ）から条件付確率ｐ（ｊ｜ｉ）を求めるものである。
【００３５】
確率算出手段２０３は、ｍ個の測度ｂjの各クラスが｛ｂ1，ｂ2，ｂ3，・・・，ｂm｝であり、各クラスの最大値がｂmaxであるとき、測度ｂiの最大値ｂmaxに対する割合をｑ（ｉ）としたとき、同時測度ｂi,jについては、最大値ｂmaxに対する割合である確率ｑ（ｉ，ｊ）を求めるものである。この確率ｑ（ｉ，ｊ）については、＜エントロピー＞と＜修正エントロピー＞と＜条件付確率ｐ（ｊ｜ｉ）に関するエントロピー＞と＜条件付確率ｐ（ｊ｜ｉ）に関する修正エントロピー＞を説明する必要がある。
【００３６】
＜エントロピー＞
エントロピー（シャノン・エントロピー）ｓは、確率分布ｐ（ｉ）に関する情報量式（１５）の平均値として式（１６）、（１７）により求めることができる。
【数２７】
0006869514.tif 000020
【００３７】
【数２８】
0006869514.tif 000021
【００３８】
＜修正エントロピー＞
本実施形態では、修正エントロピーｓ´を定義する。式（８）のｂiについて最大測度をｂmaxとする。即ち、式（１８）が成り立つ場合を条件とする。
【数２９】
0006869514.tif 000022
【００３９】
測度ｂiの最大測度ｂmaxに対する割合を、次の式（１９）に示す。
【数３０】
0006869514.tif 000023
ｑ（ｉ），ｉ＝１，２，・・・，ｎに関する情報量は次の式（２０）である。
【数３１】
0006869514.tif 000024
【００４０】
エントロピーは情報量［−ｌｏｇｐ（ｉ）］の平均（式（１６））であるのに対して、修正エントロピーｓ´を、情報量［−ｌｏｇｐ（ｉ）］から情報量［−ｌｏｇｑ（ｉ）］を差し引いた情報量の平均と定義する。即ち、修正エントロピーｓ´は、次の式（２１）、（２２）により示される。
【数３２】
0006869514.tif 000025
【００４１】
＜条件付確率ｐ（ｊ｜ｉ）に関するエントロピー＞
エントロピーｓを式（１６）に示すように、確率ｐ（ｉ）に関するエントロピーと考えたとき、条件付確率ｐ（ｊ｜ｉ）に関するエントロピーｓ（ｉ）を以下の式（２３）、（２４）のように定義することができる。
【数３３】
0006869514.tif 000026
【００４２】
＜条件付確率ｐ（ｊ｜ｉ）に関する修正エントロピー＞
上記の式（１８）に記述した最大測度ｂmaxに関して、同時測度ｂi,jのｂmaxに対する割合をｑ（ｉ，ｊ）とすると、次の式（２５）により示すことができる。
【数３４】
0006869514.tif 000027
ｑ（ｉ，ｊ），ｉ＝１，２，・・・，ｍ，ｊ＝１，２，・・・，ｍに関する情報量は次の式（２６）に示す通りである。
【数３５】
0006869514.tif 000028
条件付確率ｐ（ｊ｜ｉ）に関する修正エントロピーｓ´（ｉ）を情報量［−ｌｏｇｐ（ｊ｜ｉ）］から情報量［−ｌｏｇｑ（ｉ，ｊ）］を差し引いた情報量の平均と定義する。即ち、修正エントロピーｓ´（ｉ）は、次の式（２７）により示される。
【数３６】
0006869514.tif 000029
確率算出手段２０３は、上記式（２５）により確率ｑ（ｉ，ｊ）を求めることができる。
【００４３】
即ち、式（１３）、（１４）より、次の式（２８）が成り立つ。
【数３７】
0006869514.tif 000030
【００４４】
＜再びカオス尺度＞
カオス尺度Ｈについては、＜カオス尺度の定義＞において、式（６）に示した通りであるが、このカオス尺度Ｈは、条件付確率ｐ（ｊ｜ｉ）に関するエントロピーｓ（ｉ）として式（２３）、（２４）に示したものの平均と解釈でき、次の式（２９）〜（３１）に示す通りとなる。
【数３８】
0006869514.tif 000031
【００４５】
＜修正カオス尺度＞
本実施形態では、修正カオス尺度Ｈ´を、条件付確率ｐ（ｊ｜ｉ）に関する修正エントロピーｓ´（ｉ）として式（２７）、（２８）に示したものの平均と定義する。即ち、修正カオス尺度Ｈ´を、次の式（３２）、（３３）、（３４）から求める。
【数３９】
0006869514.tif 000032
【００４６】
修正カオス尺度算出手段２０４は、条件付確率に関する修正エントロピーｓ´（ｉ）を、第１の情報量［−ｌｏｇｐ（ｊ｜ｉ）］から第２の情報量［−ｌｏｇｑ（ｉ，ｊ）］を引いた情報量の平均と定義するとき、前記条件付確率に関する修正エントロピーｓ´（ｉ）の平均として定義された修正カオス尺度を、前記確率分布算出手段２０１により算出された確率分布ｐ（ｉ）と、前記条件付確率算出手段により求められた条件付確率ｐ（ｊ｜ｉ）と、前記確率算出手段により求められた確率ｑ（ｉ，ｊ）を用いて求めるものである。即ち、修正カオス尺度算出手段２０４は、上記式（３２）、（３３）、（３４）から修正カオス尺度Ｈ´を求める。
【００４７】
＜Ｈ´の計算法１＞
＜Ｈ´の計算法１＞は、写像τと不変測度（確率密度関数）ρが既知である場合に適用される。この計算法１による時系列データ解析装置或いは時系列データ解析用プログラムによって行われる処理を、図６のフローチャートを参照して説明する。初期設定として、写像をτ、分割数をｍ、写像の定義域をＩ＝［ａ，ｂ］とする（Ｓ１１）。
【００４８】
図７には、ｍ分割した写像τの１区間||Ｘi||について、写像後の値がどのように変化するかを表す計算を行う計算工程説明図が示されている。ステップＳ１１に続いて、修正カオス尺度Ｈ´を求める演算を行うためのレジスタｈ、レジスタｉ、レジスタｊの値をリセットする（ｈは０，ｉ，ｊは１）（Ｓ１２）。
【００４９】
ステップＳ１３からステップＳ１８までにおいてはｉをｍまで１づつ歩進させて処理を行うステップである。ステップＳ１３においてｉを１歩進させると、ステップＳ１４では次の式（３５）によりｐ（ｉ）を求める。
【数４０】
0006869514.tif 000033
【００５０】
次のステップＳ１５からステップＳ１７までにおいてはｊをｍまで１づつ歩進させて処理を行うステップである。ステップＳ１５においてｊを１歩進させると、ステップＳ１６では次の式（３６）によりｐ（ｊ｜ｉ）を求め、更に式（３７）によりｑ（ｉ，ｊ）を求める。
【数４１】
0006869514.tif 000034
【００５１】
更に、ステップＳ１６においては、修正カオス尺度Ｈ´を得るための途中計算として次の式（３８）によるｕを求める計算と、式（３９）によるｈを求める計算を行う。このｕを求める計算と、ｈを求める計算は、修正カオス尺度算出手段２０４が行う。なお、式（３８）、（３９）においては、矢印「←」の右辺の値を左辺にセットすることを示す。
【数４２】
0006869514.tif 000035
【００５２】
ステップＳ１６が終了すると、ステップＳ１７において、ｊ＝ｍとなったかを検出し、ｊ＝ｍとなっていなければステップＳ１５へ戻ってｊを１歩進させると、ステップＳ１６へと進む処理を繰り返す。このステップＳ１５からステップＳ１７の処理は、図７に示すｐ（ｊ｜ｉ）とｑ（ｉ，ｊ）との計算を示す図において、同じＸi内においてｊを歩進させて||Ｘj||のｊの値を歩進させて演算を行う処理であり、Ｘiは１つの区間に留まるものである。
【００５３】
このようにして処理を繰り返すうちにステップＳ１７を通過し、ステップＳ１８へ処理が進む。ステップＳ１８においては、ｉ＝ｍとなったかを検出し、ｉ＝ｍとなっていなければステップＳ１３へ戻ってｉを１歩進させると、ステップＳ１４へと進む処理を繰り返す。このステップＳ１３からステップＳ１８の処理は、図７に示すｐ（ｊ｜ｉ）とｑ（ｉ，ｊ）との計算を示す図において、||Ｘi||のｉの値を歩進させて演算を行う処理であって、Ｘiは現在の１区間から次の１区間へ移動して処理が行われる。１つＸiが区間を移動すると、移動した１区間においてステップＳ１５からステップＳ１７の処理が行われる。
【００５４】
上記のようにしてステップＳ１３からステップＳ１８の処理を繰り返すうちにステップＳ１８を通過すると、ステップＳ１９においては、既に求まっているｈを修正カオス尺度Ｈ´として結果出力部１０６から出力する。
【００５５】
＜Ｈ´の計算法２＞
＜Ｈ´の計算法２＞は、写像τと不変測度（確率密度関数）ρが未知であり、与えられた時系列データから修正カオス尺度Ｈ´を計算する場合に適用される。＜Ｈ´の計算法２＞では、写像τと不変測度（確率密度関数）ρが未知であることが前提であるから＜Ｈ´の計算法１＞において使用した式（３５）〜式（３７）ではなく、ｐ（ｉ）、ｐ（ｊ｜ｉ）、ｑ（ｉ，ｊ）を、与えられた時系列データから求めることになる。即ち、確率分布算出手段２０１は、データ長がｎ＋１の時系列データ｛ξ0，ξ1，ξ2，・・・，ξn｝から、確率分布ｐ（ｉ）、条件付確率ｐ（ｊ｜ｉ）、確率ｑ（ｉ，ｊ）を求める。
【００５６】
この計算法２による時系列データ解析装置或いは時系列データ解析用プログラムによって行われる処理を、図８、図９のフローチャートを参照して説明する。初期設定として、分割数をｍ、細分割数をＭ、時系列データの定義域をＩ＝［ａ，ｂ］とし、Δｘ＝（ｂ−ａ）／ｍとする（Ｓ２１）。更に、レジスタｈとレジスタｉ，ｊ，ｋの初期化を行い（Ｓ２２）、カウンタｃ1［］，ｃ2［］［］，ｃ3［］［］の初期化が行われる（Ｓ２３）。次に、データ長がｎ＋１の時系列データ｛ξ0，ξ1，ξ2，・・・，ξn｝の取り込みが行われる（Ｓ２４）。
【００５７】
ステップＳ２５からステップＳ２７までにおいてはｋをｎ−１まで１づつ歩進させて処理を行うステップである。ステップＳ２５においてｋを１歩進させ、ステップＳ２６ではｐ（ｉ）、ｐ（ｊ｜ｉ）、ｑ（ｉ，ｊ）を求めるために、カウンタｃ1［］，ｃ2［］［］，ｃ3［］［］の値を確定する。
【００５８】
具体的には、次の式（４０）、（４１）、（４２）からν1、ν2、ν3を求める。
【数４３】
0006869514.tif 000036
【００５９】
上記の通りの計算からν1、ν2、ν3が求められると、次の式（４３）、（４４）、（４５）からカウンタｃ1［］，ｃ2［］［］，ｃ3［］［］の値が求められる。
【数４４】
0006869514.tif 000037
【００６０】
上記ステップＳ２６の処理が終了すると、ステップＳ２７においてｋがｎ−１となったかを検出し、ｎ−１となっていなければ、ステップＳ２５へ戻りｋを１歩進してステップＳ２６の処理を行い、ステップＳ２７へ進みｋがｎ−１となったかを検出する処理を繰り返す。このようにして、処理を繰り返すとステップＳ２７においてｋ＝ｎ−１の成立が検出されると、図９のフローチャートに示す処理が行われる。即ち、ステップＳ３１からステップＳ３６までにおいてはｉをｍまで１づつ歩進させて処理を行うステップである。
【００６１】
ステップＳ３１においてはｉをｍまで１づつ歩進させると、ｐ（ｉ）へｃ1［ｉ］／ｎをセットし、ステップＳ３３へ進む。ステップＳ３３からステップＳ３５まではｊをｍまで１づつ歩進させて処理を行うステップである。ステップＳ３３においてｊを１歩進させ、ｐ（ｉ，ｊ）を求め、ｐ（ｊ｜ｉ）を求める。ステップＳ２６における処理によってｃ1［ｉ］，ｃ2［ｉ］［ｊ］は、それぞれ次の式（４６）、（４７）によって表されるように求まっている。
【数４５】
0006869514.tif 000038
【００６２】
上記からｐ（ｉ）は既に述べた通り、ｃ1［ｉ］／ｎがセットされた式（４８）により、ｐ（ｉ，ｊ）、ｐ（ｊ｜ｉ）は、それぞれ式（４９）、（５０）により表される式により計算される。
【数４６】
0006869514.tif 000039
【００６３】
このステップＳ３４では、時系列データからｑ（ｉ，ｊ）についても求める。ｑ（ｉ，ｊ）を求める意味は、分割区間Ｘi内にデータが存在する範囲の割合を求めることである。そこで、分割区間Ｘiを更にＭ等分した細分割区間Ｙiを用いる。即ち、区間Ｉをｍ×Ｍ等分した細分割区間をＹi（ｉ＝１，２，・・・，ｍ×Ｍ）とする。Ｉは、次の式（５１）により表すことができる。
【数４７】
0006869514.tif 000040
また、上記細分割区間をＹiと分割区間Ｘiの関係は以下の式（５２）の通りである。
【数４８】
0006869514.tif 000041
【００６４】
具体的な手法を以下に述べる。既にステップＳ２３において説明した通り、カウンタｃ3［ｉ］［ｊ］を用意し、初期化すると次の式（５３）のようになる。
【数４９】
0006869514.tif 000042
【００６５】
次に、ξk∈Ｘiの基で、ξk+1∈Ｙjとなる測度を次の式（５４）により示す。
【数５０】
0006869514.tif 000043
上記の式（５４）によりカウンタを更新する（ｃ3［ｉ］［ｊ］の最大値が１である点に注意する。）。従って、
【数５１】
0006869514.tif 000044
と、上記式（５３）により定義されるｔ（ｉ，ｊ）は、ξk∈Ｘiの基で、ξk+1∈Ｘjを満たすＸjにおいて、Ｘjに内にデータが存在する範囲の割合を意味する。なお、ｔ（ｉ，ｊ）は、式（３４）における第２項のｂi,jに相応する。ｔ（ｉ，ｊ）の最大値がＭであることから、Ｍは式（３４）における第２項のｂi,jに相応する。以上により、ｑ（ｉ，ｊ）は、上記ｔ（ｉ，ｊ）を用いて次の式（５６）に示されるように記載することができる。
【数５２】
0006869514.tif 000045
【００６６】
以上によって式（４９）、（５０）、（５６）を用いて得たｐ（ｉ，ｊ）、ｐ（ｊ｜ｉ）、ｑ（ｉ，ｊ）を式（３３）へ代入し修正カオス尺度Ｈ´を算出し、修正カオス尺度Ｈ´として結果出力部１０６から出力することができる。
【００６７】
以上のようにして＜Ｈ´の計算法２＞によりロジスティック写像から得られた修正カオス尺度と、従来のカオス尺度、更にリアプノフ指数の波形を図１０、図１１に示す。図１０では、分割数ｍが２０、細分割数Ｍが２０、データ長を定義するｎが１００００である場合を示す。また、図１１では、分割数ｍが１００、細分割数Ｍが１００、データ長を定義するｎが１００００００である場合を示す。いずれの場合も修正カオス尺度がリアプノフ指数に近似しており、時系列データ解析に極めて優れたものであることを示している。
【００６８】
＜Ｈ´の計算法３＞
本実施形態に係るデータ解析手法は、処理が所謂軽量にできているため、リアルタイムで得られるデータに対してリアプノフ指数による解析に極めて近似した解析を行うことが可能である。時系列データがリアルタイムに得られている場合に、修正カオス尺度Ｈ´を得る計算法を＜Ｈ´の計算法３＞と称し、以下に説明を行う。図１２のタイミングチャートに示されるように、時系列データは最上部の右横向き矢印により示されている通りリアルタイムで常時到来するものとする。
【００６９】
初めに時刻ｔ0から時刻ｔ1までの時間間隔Ｔ（＝ｔ1−ｔ0）における時系列データＤＡ１を抽出して第１番目の抽出データとする。この第１番目の抽出データは時刻ｔ1の直後に本実施形態に係る時系列データ解析用プログラムに渡され、処理が行われて第１番目の修正カオス尺度Ｈ1´が得られ、その後に結果出力部１０６から出力される。
【００７０】
次に、時刻ｔ0＋Δｔから時刻ｔ1＋Δｔまでの時間間隔Ｔ（＝ｔ2−ｔ1）における時系列データＤＡ２を抽出して第２番目の抽出データとする。この第２番目の抽出データは時刻ｔ1＋Δｔの直後に本実施形態に係る時系列データ解析用プログラムに渡され、処理が行われて第２番目の修正カオス尺度Ｈ2´が得られ、その後に結果出力部１０６から出力される。
【００７１】
次に、時刻ｔ0＋２Δｔから時刻ｔ1＋２Δｔまでの時間間隔Ｔ（＝ｔ3−ｔ2）における時系列データＤＡ３を抽出して第３番目の抽出データとする。この第３番目の抽出データは時刻ｔ1＋３Δｔの直後に本実施形態に係る時系列データ解析用プログラムに渡され、処理が行われて第３番目の修正カオス尺度Ｈ3´が得られ、その後に結果出力部１０６から出力される。
【００７２】
更に、同様に第４番目の抽出データから第４番目の修正カオス尺度Ｈ4´が得られ、以下同様に処理が続けられる。即ち、本実施形態では、データ長がｎ＋１の時系列データ｛ξ0，ξ1，ξ2，・・・，ξn｝をサンプリング期間をずらして複数回サンプリングし、上記修正カオス尺度算出手段２０４は、サンプリングした複数回分の時系列データ毎に修正カオス尺度算出を行う。斯くして、本実施形態では、次々に所要時間修正カオス尺度Ｈ4´が得られ、それぞれの時間単位でデータ解析を行うことができる。
【００７３】
なお、上記実施の形態では、ｋ番目の抽出データとｋ＋１番目の抽出データとがオーバーラップする部分があるが、ｋ番目とｋ＋１番目の間の抽出時間間隔Δｔにおいて大量のデータが到来するような場合には、オーバーラップ部分を設けることなく計算を行っても良いことは勿論である。
【００７４】
サンプリング期間をずらして複数回サンプリングして、サンプリングした複数回分の時系列データ毎に修正カオス尺度算出を行う（第１の処理態様という）のではなく、修正カオス尺度算出手段２０４は、データ長がｎ＋１の時系列データ｛ξ0，ξ1，ξ2，・・・，ξn｝を１度得て修正カオス尺度算出を行うよう（第２の処理態様という）にしても良い。また、第１の処理態様と第２の処理態様と双方を有し、データの長さや計算を行う時間帯に応じて、更には、オペレータの指示により、上記の２態様の一方を採用するようにしても良い。
【００７５】
本願の発明者らは、非特許文献１２において、心拍間隔データのカオス尺度と自律神経活動の関連について記述し、非特許文献１３において、心拍変動のカオス性と副交感神経活動との関係について記述した。このように心拍変動のカオス性と眠気や怠慢等の生理状態とは相関が見られることから、本実施形態の解析装置や解析用プログラムは、生体の生理状態について解析する場合に極めて有用な修正カオス尺度を時系列データから得ることができるものである。
【００７６】
電子機器等においては、電源電圧の不規則な揺らぎ（カオス的な変化）は、電圧供給先である内部の電子部品乃至素子等の抵抗値等に係る電気的特性が不規則に変動している可能性を示唆するものである。従って、電子機器等の電源電圧を時系列データとして本実施形態に係る解析装置や解析用プログラムを用いて解析を行うことによって電子部品乃至素子等の異常や故障の検出並びに予測へつながることが期待される。
【００７７】
更に、機械学習やディープラーニングは、入力データの特徴Ａに対し人間や自然界或いは機械・装置において生じる事象Ｂとの対応関係を過去のデータを用いて学習しておき、新たな入力データｎｅｗＡに対して最も生じる可能性が高い事象ＢＢを自動的に割り出すものと言える。そこで、本実施形態に係る解析装置や解析用プログラムを用いて入力データの特徴Ａを解析しておき、これに対する事象Ｂを対応付けて、機械学習やディープラーニングへ与えて学習させ、その後に到来する入力データについての、本実施形態に係る解析装置や解析用プログラムを用いて得た特徴から最も生じる可能性が高い事象の推定などに適用でき得るものと考えられる。
【００７８】
まとめるならば、機械学習を含む一般的なデータ解析に用いられる数学的手法は、データの規則性を見出すことを主目的として発展してきた背景がある。これに対し、本実施形態に係る解析装置や解析用プログラムのように規則性とは逆の事象である不規則性（カオス性）を把握して事象に対する特徴とする手法は一般的に行われているものではない。従って、本実施形態に係る解析装置や解析用プログラムを用いて得た不規則性の特徴を用いるものは、発生事象の予測や発生している事象の検出など幅広い応用が期待できる。
【符号の説明】
【００７９】
１００  時系列データ解析装置
１０２  主メモリ
１０３  バス
１０４  外部記憶装置
１０５  時系列データ供給部
１０６  結果出力部
２０１  確率分布算出手段
２０２  条件付確率算出手段
２０３  確率算出手段
２０４  修正カオス尺度算出手段
【図１】
0006869514.tif 000054
【図２】
0006869514.tif 000055
【図３】
0006869514.tif 000056
【図４】
0006869514.tif 000057
【図５】
0006869514.tif 000058
【図６】
0006869514.tif 000059
【図７】
0006869514.tif 000060
【図８】
0006869514.tif 000061
【図９】
0006869514.tif 000062
【図１０】
0006869514.tif 000063
【図１１】
0006869514.tif 000064
【図１２】
0006869514.tif 000065