JP5826892B1 - 変化点検出装置、変化点検出方法、及びコンピュータプログラム - Google Patents

Abstract

【課題】簡易な方法で外れ値の影響を低減させて変化点を検出することができる変化点検出装置、変化点検出方法、及びコンピュータプログラムを提供する。【解決手段】時系列的に連続する事象に対応する一連の値からなる系列を入力し、前記値を時系列に記録する記録手段と、前記系列における事象の変化点を検出する検出手段とを備える変化点検出装置において、前記検出手段は、前記記録手段に記録されている前記系列の記録期間中におけるトレンドをトレンドモデルに基づき算出するトレンド算出手段と、該トレンド算出手段が算出したトレンドの勾配を算出する勾配算出手段と、算出された勾配のマハラノビスの汎距離を算出する汎距離算出手段とを備え、前記検出手段は、算出されたマハラノビスの汎距離に基づき変化点を検出する。【選択図】図９

Description

本発明は、時系列的に変化する事象を分析して将来の事象を予測する方法に関し、特に、時系列的な事象に対応する値の系列から、外れ値の影響を低減させて事象の変化点を高精度に検出することができる変化点検出装置、変化点検出方法、及びコンピュータを変化点検出装置として機能させるコンピュータプログラムに関する。

経時的に変化する事象について、該事象に係る任意の時点までの時系列データを用いて以後の変化を予測することは様々な分野で有用である。例えば為替、株価の過去の推移の観測に基づいて一日後、一週間後、半年後の為替又は株価の予測を行なったり、製品の販売実績に基づいて次期の需要予測を求めたりすることは、人間の経験及び判断力に基づいて行なわれている。上述の為替、株価、製品需要などの予測は非常に困難であるが、それらの予測が人間の主観を除去した上で高精度に実現することが可能であるとすれば、生産、販売及び在庫管理、並びに物流、製品開発等の計画に非常に有用となる。

発明者はこれまでに、時系列的に変化する種々の事象に関する予測方法を提案してきた（特許文献１等）。既に提案した予測方法により、種々の事象の予測を実現できることを示されているが、更に予測精度を高めることが求められる。発明者が提案してきた予測方法の精度を高めるためには、過去の期間における時系列的な変化傾向（トレンド）を正確に捉えるべく、予測に用いる観測情報から適切に情報を選択することが有効である。情報を選択するための基準として、時系列的に変化する事象について変化点として定義される時点を特定することが有用である。変化点は言わば事象の潮目が変わる時点であり、観測される事象が実際に変化する時点、事象を発現させる状態が変化する時点のことである。変化点の検出は、データマイニング等の各種分野でも有用とされる。

変化点は、観測系によるノイズ、又は突発的な原因による偶然の変動（所謂外れ値）の影響を受ける。ノイズ、外れ値の除去、及び、外れ値の影響を受けない変化点の検出に関しては、統計学、機械学習、データマイニング等の種々の分野にて多くの提案がなされている。特に時系列的な変化を考慮した変化点の検出に関しては、時系列データの性質に対応する統計モデルを用いて外れ値及び変化点を統一的に検出する方法が特許文献２に開示されている。

特許第５０６８３８２号公報 特許第３８２１２２５号公報

特許文献２に開示された方法では、過去の事象の系列に依存する確率密度に基づき、忘却型の自己回帰モデル学習を行なう方法が利用されている。しかしながら、学習におけるパラメータの最適化等に課題が残される。

外れ値及び変化点の検出方法には、統計学的にマハラノビスの汎距離を利用した手法が用いられることがある。しかしながら、マハラノビスの汎距離を利用する場合、変化点を外れ値（ノイズ）と区別して検出することが困難な場合がある。

発明者は、マハラノビスの汎距離を利用しつつも、時系列的に変化する事象のトレンドを加味することによって容易に変化点の検出が可能となるとの知見を得た。

本発明は斯かる知見に基づいてなされたものであり、簡易な方法で外れ値の影響を低減させて変化点を検出することができる変化点検出装置、変化点検出方法、及びコンピュータを変化点検出装置として機能させるコンピュータプログラムを提供することを目的とする。

本発明に係る変化点検出装置は、時系列的に連続する事象に対応する一連の値からなる系列を１又は複数入力し、１又は複数の前記系列を時系列的に記録する記録手段と、前記事象の記録期間中における変化点を検出する検出手段とを備える変化点検出装置において、前記検出手段は、前記記録手段に記録されている前記系列に基づく時系列データの記録期間中におけるトレンドをトレンドモデルに基づき算出するトレンド算出手段と、該トレンド算出手段が算出したトレンドの勾配を算出する勾配算出手段と、算出された勾配のマハラノビスの汎距離を算出する汎距離算出手段とを備え、前記検出手段は、算出されたマハラノビスの汎距離に基づき変化点を検出するようにしてあることを特徴とする。

本発明に係る変化点検出装置は、前記トレンド算出手段は、前記記録手段に記録されている複数系列に対してマハラノビスの汎距離を算出して前記時系列データとする時系列データ化手段を更に備え、算出された時系列データの記録期間中におけるトレンドをトレンドモデルに基づき算出するようにしてあることを特徴とする。

本発明に係る変化点検出装置は、前記トレンド算出手段は、前記事象に対応する一連の値が、該値を発現させる状態が前記値に与える影響を表現する関数と白色雑音との和で記述される状態空間モデルに、各時点における前記一連の値がトレンドと白色雑音との和で記述されるトレンドモデル、及び前後する時点におけるトレンドは略等しいとするモデルを適用し、状態空間モデルの解法に基づきトレンドを算出することを特徴とする。

本発明に係る変化点検出装置は、前記勾配算出手段は、前後する時点の内、前の時点におけるトレンドから後の時点におけるトレンドへの変化率を算出するようにしてあることを特徴とする。

本発明に係る変化点検出装置は、前記検出手段は、前記汎距離算出手段により算出されたマハラノビスの汎距離を、前記記録期間中におけるマハラノビスの汎距離の最大値で除算した値を求める手段を更に備えることを特徴とする。

本発明に係る変化点検出方法は、時系列的に連続する事象に対応する１又は複数の一連の値からなる系列を入力し、１又は複数の系列を時系列的に記録する記録手段を備えるプロセッサが、前記事象の記録期間中における変化点を検出する変化点検出方法において、前記プロセッサは、前記記録手段に記録されている前記系列に基づく時系列データの記録期間中におけるトレンドをトレンドモデルに基づき算出し、算出されたトレンドの勾配を算出し、算出された勾配のマハラノビスの汎距離を算出し、算出されたマハラノビスの汎距離に基づき変化点を検出することを特徴とする。

本発明に係るコンピュータプログラムは、時系列的に連続する事象に対応する１又は複数の一連の値からなる系列を入力し、１又は複数の系列を時系列的に記録する記録手段を備えるコンピュータに、前記事象の記録期間中における変化点を検出させるコンピュータプログラムにおいて、前記コンピュータに、前記記録手段に記録されている前記系列に基づく時系列データの記録期間中におけるトレンドをトレンドモデルに基づき算出するステップと、算出されたトレンドの勾配を算出するステップと、算出された勾配のマハラノビスの汎距離を算出するステップと、算出されたマハラノビスの汎距離に基づき変化点を検出するステップとを実行させることを特徴とする。

本発明では、単一系列又は複数系列の時系列データにおけるトレンドが算出され、トレンドの勾配の周辺との相関強度（類似性）がマハラノビスの汎距離により算出される。変化点は、それまでの時系列データの推移の傾向と比して類似性が低いためにマハラノビスの汎距離が大きく算出されることから、算出結果により変化点が検出可能となる。

本発明では、複数系列の値のマハラノビスの汎距離が算出されて単一系列の時系列データ化され、トレンドが算出される。マハラノビスの汎距離を算出することにより、複数系列の値が影響する事象の総合的なトレンドが算出される。算出された総合的なトレンドの勾配の類似性が算出され、算出結果により複数系列の値についての変化点が検出可能となる。

本発明では、トレンド算出において、入力される系列を観測系列とする状態空間モデルに、各観測系列は各時点のトレンドと白色雑音との和で記述されるというトレンドモデル、及び前後する時点のトレンドは略等しいというモデルを適用して、状態空間モデルの解法に基づきトレンドが推定算出される。これにより、トレンドの勾配はトレンド変化に敏感に反応する値として算出することが可能となり、トレンドが変化していない場合は、変化の度合いが低く抑えられるように算出される。

本発明では、トレンドの勾配として、前後する時点におけるトレンド間において後の時点のトレンドへの変化率が算出される。時間的に後の時点のトレンドへのトレンド増の変化が正値として算出される。

本発明では、記録期間に亘り算出されたマハラノビスの汎距離が、最大値により正規化される。これにより、最大値との比較によって変化点が検出される。

本発明による場合、トレンドを考慮し、大きな変化を伴う外れ値についてはトレンドが変化したわけではないことを示す変化点スコアが算出され、これにより、外れ値の影響を低減させて簡易な方法で変化点を検出することができる。

実施の形態１における検出装置の構成を示すブロック図である。 実施の形態１の検出装置による検出処理手順の一例を示すフローチャートである。 トレンド勾配の算出概念を示すグラフである。 実施例１の時系列データＤ_t を示すグラフである。 実施例１の時系列データＤ_t に対して算出された信号トレンドを示すグラフである。 図５の一部を拡大した図である。 実施例１の時系列データＤ_t に対して算出された信号トレンドの勾配を示すグラフである。 実施例１の時系列データＤ_t の信号トレンドの勾配についてのＭＤを示すグラフである。 実施例１の時系列データＤ_t に対して導出された変化点スコアを示すグラフである。 実施例１の時系列データＤ_t に対して直接的にＭＤを算出した場合のグラフである。 実施例１の時系列データＤ_t の変化差分を示すグラフである。 図１１の変化差分に対して算出されたＭＤを示すグラフである。 図１２のＭＤを最大値で正規化して算出した変化点スコアを示すグラフである。 実施の形態１の検出装置により導出された変化点スコアと、対応する信号値の変化量との関係を示す散布図である。 比較例にて算出された変化点スコアと、対応する信号値の変化量との関係を示す散布図である。 実施の形態１の検出装置により導出されるＣＩ一致指数の変化点スコアを示すグラフである。 実施の形態１の検出装置により導出される大口電気使用量の変化点スコアを示すグラフである。 実施の形態１の検出装置により導出される日経平均株価の変化点スコアを示すグラフである。 実施の形態１の検出装置により導出される消費者物価指数の変化点スコアを示すグラフである。 実施の形態２の検出装置による検出処理手順の一例を示すフローチャートである。 複数系列の信号値の例を示す説明図である。 実施例３の複数系列の信号値に対して算出されたＭＤ及びＭＤのトレンドを示すグラフである。 実施例３の複数系列の信号値に対して算出されたＭＤトレンドの勾配を示すグラフである。 実施例３の複数系列の信号値から得られたＭＤトレンドの勾配のＭＤを示すグラフである。 実施例３の複数系列の信号値に対して導出された変化点スコアを示すグラフである。 実施例３にて導出された変化点スコアをＭＤトレンドと比較して示すグラフである。 実施例３にて導出された変化点スコアと、複数系列の信号値に係る指標値とを比較するグラフである。

以下、本発明の実施の形態を、図面に基づいて具体的に説明する。

（実施の形態１）
図１は、実施の形態１における検出装置１の構成を示すブロック図である。検出装置１は例えばパーソナルコンピュータ、サーバコンピュータ等のコンピュータを用いる。検出装置１は、制御部１０、記録部１１、一時記憶部１２、入力部１３及び出力部１４を備える。

制御部１０は、ＣＰＵ（Central Processing Unit ）を用いる。制御部１０は、以下に説明する検出プログラム１Ｐに基づき、パーソナルコンピュータを制御し、本実施の形態における検出装置１としての機能を発揮させる。

記録部１１は、ＲＯＭ（Read Only Memory）、ハードディスクドライブ等の不揮発性メモリを用いる。なお、記録部１１は、外付けのハードディスクドライブ、光学ディスクドライブ、通信網を介して接続される他の記録装置であってもよい。すなわち、記録部１１とは、制御部１０からアクセス可能な１又は複数の情報記録媒体の総称である。

記録部１１には、本実施の形態の変化点検出方法を実現するための各種手順を含む検出プログラム１Ｐが記録されている。また、記録部１１の記録領域の一部は、信号値を記録する領域として用いられる。制御部１０は、記録部１１に対して情報の読み書きが可能である。

一時記憶部１２は、例えば、ＤＲＡＭ（Dynamic Random Access Memory）、ＳＲＡＭ（Static Random Access Memory ）等の不揮発性メモリである。一時記憶部１２は、制御部１０の処理によって発生した情報を一時的に記憶する。

入力部１３は、キーボード、マウス等を用い、ユーザの操作に基づき情報を入力する。

出力部１４は、液晶モニタなどの表示部、又はプリンタ等の印刷部を用い、制御部１０による情報の処理結果を出力する。

このように構成される検出装置１にて、制御部１０が検出プログラム１Ｐに基づく処理を実行することにより、変化点を検出する。

図２は、実施の形態１の検出装置１による検出処理手順の一例を示すフローチャートである。

制御部１０は、入力部１３から時系列データである複数の信号値を入力する。制御部１０は、入力した複数の信号値を記録部１１に記録する（ステップＳ１）。なお、信号値は入力部１３から入力されるもののみならず、通信網を介して制御部１０が信号値を入力するようにしてもよいし、他の情報記録媒体から信号値を入力するようにしてもよい。

制御部１０は、記録した複数の信号値（観測値）に対し、隣り合う時点のトレンド成分は略等しいというトレンドモデルに基づいて記録期間全体における信号トレンドを算出する（ステップＳ２）。ステップＳ２における信号トレンドの算出には、状態空間モデルによる１次のトレンド成分を適用する。具体的には制御部１０は、トレンドｔ_n について、式（１）にて表現されるようにトレンドが時間軸において局所的に略一定の値をとり、前のタイミングにおけるトレンドｔ_n-1 に対してｔ_n ≒ｔ_n-1 であるというモデル（ランダムウォークモデル、参考：トレンド成分モデル−構造の確率的変化のモデル−、北川源四郎、１９９３、「ＦＯＲＴＲＡＮ７７時系列解析プログラミング」、pp.252-263）を、ステップＳ１で入力した信号値を観測値とする状態空間モデルに当てはめ、カルマンフィルタ及び平滑化処理等の処理を実行することによってトレンドを算出する。

制御部１０は、ステップＳ２で算出したトレンドに対し、風上差分からのトレンド勾配を算出する（ステップＳ３）。ステップＳ３において詳細には、制御部１０は、各時点について算出されたトレンドに対し、対象としている一の時点の次の時点におけるトレンドへの変化率を算出する。図３は、トレンド勾配の算出概念を示すグラフである。横軸は時間軸を示し、縦軸は算出されたトレンドを示す。ステップＳ３にて算出されるトレンド勾配は、図３中のＧに対応する。

制御部１０は、ステップＳ３で算出したトレンド勾配について、マハラノビスの汎距離（以下、ＭＤ：Maharanobis’s Distance と言う）を算出する（ステップＳ４）。これにより、トレンド勾配の類似性が求められる。なお、算出されたＭＤが近いほどトレンド（変化傾向）が類似し、遠いほどトレンドが非類似であることが示される。

制御部１０は、ステップＳ４で算出したＭＤに基づき変化点スコアを算出する（ステップＳ５）。詳細には、制御部１０は、入力した信号値の全期間に亘って算出したＭＤの最大値で、算出されたＭＤ全てを除算する。ステップＳ５の処理により、期間全体に亘って算出したＭＤの内の最大値に対応する時点の変化点スコアが１．０となるようにして各時点での変化点スコアが算出される。これにより、非類似性が高いトレンドの勾配に係る時点が変化点として検出される。

このようにして時系列データである信号値に対して導出される変化点スコアに基づき、信号値の変化点を検出することができる。検出装置１により導出される変化点スコアは、トレンドが考慮されているから、トレンドに従って真に状態が変化したことに応じて推定導出されるものである。このとき上述したように、実施の形態１における検出装置１では、算出した値が基準とする所定の値以上であるか否か等の判定条件は不要である。

次に、検出装置１による変化点検出方法を適用した具体的な実施例について説明する。

（実施例１）
実施例としてまず、人工的に発生させた時系列データＤ_t に本発明を適用し、外れ値の影響を低減させることができることを示す。

図４は、実施例１の時系列データＤ_t を示すグラフである。図４の横軸は時間軸に相当するデータの番号を示し、縦軸は信号値を示している。図４に示す時系列データＤ_t は大きく１０段階に変化し、各段階で細かく変動する信号値を人工的に発生させたものである。図４中のＡ〜Ｉの９つの時点が変化点に対応する。また、実施例１の時系列データＤ_t には、外れ値として扱われるべき局所的に変動する値が含まれている。図４中の符号ａ，ｂ，ｃにより外れ値を示している。

検出装置１の制御部１０によって、図４に示した時系列データＤ_t における変化点Ａ〜Ｉを変化点として検出し、その際に外れ値ａ，ｂ，ｃの影響を低減させることができることを以下に示す。

検出装置１の制御部１０は、時系列データＤ_t を入力し（Ｓ１）、入力した時系列データＤ_t に対して信号トレンドを算出する（Ｓ２）。

図５は、実施例１の時系列データＤ_t に対して算出された信号トレンドを示すグラフであり、図６は、図５の一部を拡大して示している。図５及び図６の横軸は時間軸に相当するデータの番号を示し、縦軸はトレンド値（信号値）を示している。図５及び図６中では、太線にて信号トレンドが示され、細線にて演算前の時系列データＤ_t の信号値が示されている。信号トレンドの算出により、図５に示すように、変化点Ａ〜Ｉにおけるトレンドの変化の幅は、元の信号値の変化の幅と同程度であるものの、外れ値ａ，ｂ，ｃに対応するトレンドの変化は、元の信号値の変化の幅の１／３程度に低減されている。

次に検出装置１の制御部１０は、算出した信号トレンドの勾配を算出する（Ｓ３）。

図７は、実施例１の時系列データＤ_t に対して算出された信号トレンドの勾配を示すグラフである。なお図７には比較のために図５のトレンド値（信号値）を併せて示している。図７の横軸は、時間軸に相当するデータの番号を示している。図７の縦軸は、左側にトレンド値（信号値）を示し、右側に算出されたトレンド勾配に対応する値を示している。図７では、変化点において値が増加する場合は勾配が正、値が減少する場合は勾配が負で示されていることがわかる。

検出装置１の制御部１０は、算出した勾配についてＭＤを算出する（Ｓ４）。

図８は、実施例１の時系列データＤ_t の信号トレンドの勾配についてのＭＤを示すグラフである。図８の横軸は時間軸に相当するデータの番号を示し、縦軸（左側）はＭＤを示す。なお図８には比較のために図７の信号トレンドの勾配を下部に示し、縦軸（右側）にトレンド勾配に対応する値を示している。図８に示すように、外れ値ａ，ｂ，ｃに対応するＭＤは、変化点Ａ〜Ｉに対応するＭＤと比較して値が低く算出され、特に外れ値ａに関しては、他の変化点Ａ〜Ｉと比較して無視できる程度までに値が小さくなっている。

検出装置１の制御部１０は、算出したＭＤを最大値によって除算することによって変化点スコアを算出する（Ｓ５）。

図９は、実施例１の時系列データＤ_t に対して導出された変化点スコアを示すグラフである。図９の横軸は時間軸に相当するデータの番号を示し、縦軸（右側）に変化点スコアを示す。最大値によって除算されているから、変化点スコアは最大値１．０である。なお図９には比較のために図４の時系列データＤ_t の信号値を示し、縦軸（左側）に信号値に対応する値を示している。図９に示すように、変化点Ａ〜Ｉ、及び外れ値ａ，ｂ，ｃ夫々に対応する変化点スコアが導出されていることが確認される。

このようにして、実施の形態１における検出装置１により、時系列データである信号値に対し、外れ値の影響を低減させて変化点を検出することが可能である。

比較として、時系列データＤ_t に対して信号トレンドを算出することなしに、ＭＤを算出した場合の算出結果を挙げる。

図１０は、実施例１の時系列データＤ_t に対して直接的にＭＤを算出した場合のグラフである。図１０の横軸は時間軸に相当するデータの番号を示し、縦軸（右側）はＭＤを示す。図１０には比較のために、図４の時系列データＤ_t の信号値を上部に示し、縦軸（左側）に信号値に対応する値を示している。図１０に示すように、全体として値は小さくなるが変化点Ａ〜Ｉと外れ値を同様にして検出することは困難である。

次に他の比較例として、ＭＤを算出する前に信号値の変化差分（風上差分、次の信号値への変化分）を算出した場合の算出結果を挙げる。

図１１は、実施例１の時系列データＤ_t の変化差分を示すグラフである。図１１の横軸は時間軸に相当するデータの番号を示し、縦軸（右側）は変化差分を示す。図１１には比較のために図４の時系列データＤ_t の信号値を上部に示している。図１２は、図１１の変化差分に対して算出されたＭＤを示すグラフである。図１２の横軸は時間軸に相当するデータの番号を示し、縦軸（右側）はＭＤを示す。図１２においても比較のために図４の時系列データＤ_t の信号値を上部に示し、縦軸（左側）に信号値に対応する値を示している。図１３は、図１２のＭＤを最大値で正規化して算出した変化点スコアを示すグラフである。図１３の横軸は時間軸に相当するデータの番号を示し、縦軸（右側）は変化点スコアを示す。図１３においても比較のために図４の時系列データＤ_t の信号値を上部に示し、縦軸（左側）に信号値に対応する値を示している。図８のグラフと比較した場合に、特に外れ値ｃについて他の変化点Ｈ，Ｉ等と区別することが困難であることがわかる。

変化点と外れ値とを区別できるか否かについて、実施の形態１の検出装置１により導出された変化点スコアと、比較例にて算出された変化点スコアとを比較して説明する。図１４は、実施の形態１の検出装置１により導出された変化点スコアと、対応する信号値の変化量との関係を示す散布図である。これに対し図１５は、比較例にて算出された変化点スコアと、対応する信号値の変化量との関係を示す散布図である。図１４及び図１５では、横軸に信号値の変化量を示し、縦軸に変化点スコアを示している。図１４及び図１５では、変化点Ａ〜Ｉを黒丸にて示し、外れ値ａ，ｂ，ｃを白抜きの菱形にて示している。変化量は、各時点における前又は後の値からの振れ幅である。例えば図４を参照した場合、変化点Ａの変化量は２．５、変化点Ｂの変化量は３より少し大きい程度、変化点Ｃの変化量は略ゼロ、変化点Ｄの変化量は７より少し小さい程度、変化点Ｅの変化量は１０より大きく、変化点Ｆの変化量は１０より少し小さい。変化点Ｇの変化量は７．５程度、変化点Ｈの変化量も７前後、変化量Ｉの変化量も７．５程度である。また、外れ値ａの変化量は２．５程度、外れ値ｂの変化量は６程度、外れ値ｃの変化量は１０より少し小さく見える。定量的には変化量は、各段階における信号値の平均値に対する差分の絶対値として求めた。

図１４に示す散布図では、変化点Ａ〜Ｉ及び外れ値ａ，ｂ，ｃの変化点スコアを変化量に対してプロットした場合、図１４に示すように、変化点Ａ〜Ｉの変化点スコアに対する実線の近似曲線と、外れ値ａ〜ｃの変化点スコアに対する破線の近似曲線とは、有意に分離されている。これに対し、図１５に示す散布図では、変化点Ａ〜Ｉの変化点スコアに対する実線の近似曲線と、外れ値ａ〜ｃの変化点スコアに対する破線の近似曲線とが重なっており、変化点Ａ〜Ｉと外れ値ａ〜ｃとを区別することが困難であることが示されている。つまり、図１５に示すように、トレンドを考慮することなしにＭＤにより変化点の検出を行なう場合（図１１〜図１３に示した方法）では、変化量が大きい外れ値については特に変化点スコアも大きく算出されてしまう。一方で本実施の形態１の検出装置１による場合、図１４に示すように、トレンドを考慮して変化点の検出を行なうときには、変化量の大きさに惑わされない変化点の検出が可能であると言える。これは、実施の形態１の検出装置１により、大きな変化量を持つ外れ値ｃに対してトレンドが変化したわけではないことを示す変化点スコアが導出されているということを示している。

ランダムウォークモデルというものは本来、ある時点の状態から他の状態への変化傾向は、任意の時点で不定であって次の時点での状態確率は現在の状態に対して正規分布するものであり、時系列データは状態確率の連鎖であるから、任意の時点で次の状態を予測することは不可能であることを記述していると言える。これに対し、本願発明者は、時系列的なデータに関しては、時間的に隣り合う信号値間の関係というのは、局所的にはほぼ等しいながらも時系列的には、ある程度の変化傾向を持ちながら変化していくものであることを表現することにも適しているとして、時系列データにおける外れ値（ノイズ）を除いた変化点を検出する方法に適用可能であるとの知見を得た。発明者はこのような知見に基づき上述したような方法を発明したものである。本発明による場合、本当にトレンドが変化したことに応じて観測された信号値に対応する変化点が推定検出される。

（実施例２）
次に実際の経済データに対して実施の形態１の検出装置１による変化点検出を実行した場合の例を示す。第１に、内閣府発表の景気動向指数の内のＣＩ（Composite index ）一致指数に対して変化点スコアを導出した。図１６は、実施の形態１の検出装置１により導出されるＣＩ一致指数の変化点スコアを示すグラフである。図１６では、変化点スコアにＣＩ一致指数を併せて示している。図１６は、横軸に時間軸を示し、縦軸（左側）にＣＩ一致指数を示し、縦軸（右側）に変化点スコアの大きさを示している。用いられたＣＩ一致指数は、１９８５年１月から２０１３年５月までの期間における各月の値である。図１６に示すように、２００８年後半の景気後退に対応する変化点と２００９年後半の景気回復に対応する変化点とが検出されていることがわかる。なお、導出された変化点スコアは、１９８５年１月から２０１３年１２月までの期間中におけるトレンド勾配のＭＤの最大値で正規化された値であるから、記録期間中における変化点の変化点スコアは、２００８年後半の最大の変化点との比較で算出される。

第２に、検出装置１により、内閣府発表の景気動向指数の内の大口電気使用量に対して変化点スコアを導出した。図１７は、実施の形態１の検出装置１により導出される大口電気使用量の変化点スコアを示すグラフである。図１７においても図１６同様に、変化点スコアに大口電気使用量を併せて示している。図１７は、横軸に時間軸を示し、縦軸（左側）に大口電気使用量（ｋＷ単位）を示し、縦軸（右側）に変化点スコアの大きさを示している。用いられた大口電気使用量は、２００３年１月から２０１３年５月までの期間における各月の値である。図１７に示すように、大口電気使用量の変化点スコアの導出により、２００８年後半の景気後退に対応する変化点と２００９年後半の景気回復に対応する変化点とが検出されていることがわかる。

第３に、検出装置１により、日経平均株価に対して変化点スコアを導出した。図１８は、実施の形態１の検出装置１により導出される日経平均株価の変化点スコアを示すグラフである。図１８においても図１６同様に、変化点スコアに日経平均株価を併せて示している。図１８は、横軸に時間軸を示し、縦軸（左側）に日経平均株価を示し、縦軸（右側）に変化点スコアの大きさを示している。用いられた日経平均株価は、２００３年１月から２０１３年５月までの期間における各月の平均値である。図１８に示すように、日経平均株価の変化点スコアの導出により、２００８年後半の景気後退に対応する変化点と２００９年後半の景気回復に対応する変化点と、更には２０１３年初頭の株価上昇に対応する変化点とが検出されていることがわかる。

第４に、検出装置１により、消費者物価指数に対して変化点スコアを導出した。図１９は、実施の形態１の検出装置１により導出される消費者物価指数の変化点スコアを示すグラフである。図１９においても図１６同様に、変化点スコアに消費者物価指数を併せて示している。図１９は、横軸に時間軸を示し、縦軸（左側）に消費者物価指数を示し、縦軸（右側）に変化点スコアの大きさを示している。用いられた消費者物価指数は、２０００年１月から２０１４年１月までの期間における各月の平均値である。図１９に示すように、消費者物価指数の変化点スコアの導出により、２００８年前半及び後半の景気後退に対応する変化点が検出されていることがわかる。

このようにして、時系列的に隣り合う時点におけるトレンドは局所的には略等しいというモデルに基づいてトレンドを算出し、トレンドの勾配、勾配の類似性をＭＤにより算出することによって、トレンドが変化したと判断される変化点の検出が可能となる。

（実施の形態２）
実施の形態１では、検出装置１は一系列の時系列データＤ_t に対して変化点スコアを導出する構成とした。これに対し実施の形態２では、検出装置１は複数系列の時系列データに対して総合的な変化点スコアを導出する。実施の形態２における検出装置１の構成は、以下に示す処理手順以外は実施の形態１における構成と同一であるから、共通する構成には同一の符号を付して詳細な説明を省略する。

図２０は、実施の形態２の検出装置１による検出処理手順の一例を示すフローチャートである。なお図２０のフローチャートに示す処理手順の内、実施の形態１の図２のフローチャートに示す処理手順と共通するステップについては同一のステップ番号を付して詳細な説明を省略する。

検出装置１の制御部１０は、入力部１３から同一期間における時系列データである複数系列の信号値を各入力し、入力した複数系列の信号値を系列ごとに記録部１１に記録する（ステップＳ２１）。なお、信号値は入力部１３から入力されるもののみならず、通信網を介して制御部１０が信号値を入力するようにしてもよいし、他の情報記録媒体から信号値を入力するようにしてもよい。

制御部１０は、記録した複数系列の信号値についてのＭＤを算出する（ステップＳ２２）。具体的に制御部１０は、複数系列を夫々次元とする多変量ベクトルとし、対応する時点における各系列の信号値をベクトルの各成分として演算を行なう。

次に制御部１０は、算出される時系列のＭＤに対し、隣り合う時点のＭＤのトレンド成分は略等しいというトレンドモデルに基づいて記録期間全体におけるＭＤのトレンドを算出する（ステップＳ２３）。

ステップＳ２３におけるＭＤトレンドの算出は、実施の形態１の信号トレンドのモデルにおける信号トレンドに代替して行なう。つまり、ステップＳ２３におけるＭＤトレンドの算出には、状態空間モデルによる１次のトレンド成分を適用する。具体的には制御部１０は、ＭＤトレンドＭＤＴ_n について、式（２）にて表現されるようにトレンドが時間軸において局所的に略一定の値をとり、前のタイミングにおけるトレンドＭＤＴ_n-1 に対してＭＤＴ_n ≒ＭＤＴ_n-1 であるというモデル（ランダムウォークモデルを、ステップＳ２２で算出された時系列のＭＤを観測値とする状態空間モデルに当てはめ、カルマンフィルタ及び平滑化処理等の処理を実行することによってトレンドを算出する。

制御部１０は、ステップＳ２３にて算出されたＭＤトレンドに対してＭＤトレンドの勾配を算出する（ステップＳ２４）。ステップＳ２４における勾配の算出は、実施の形態１の図２のフローチャートに示したステップＳ３の処理同様に、各時点について、次の時点におけるトレンドへの変化率を算出することである。

以後制御部１０は、算出した勾配のＭＤを算出し（Ｓ４）、算出したＭＤに基づき変化点スコアを算出する（Ｓ５）。

このようにして実施の形態２の検出装置１は、複数系列の信号値について各系列で変化点を各検出するのではなく、複数の項目から総合的に判断される変化点を検出するために、複数系列の信号値から変化点スコアを導出する。

以下、実施の形態２の検出装置１による変化点検出方法を適用した具体的な実施例について説明する。

（実施例３）
図２１は、複数系列の信号値の例を示す説明図である。図２１には、消費者物価指数を構成する６８項目が示されている。実施例３では、６８項目の月次データ（６８系列の時系列データ）を複数系列の信号値とし、複数系列の信号値の総合的な変化点スコアを検出装置１により導出する。

検出装置１の制御部１０は、複数系列の信号値を入力し（Ｓ２１）、入力した複数系列の信号値のＭＤを算出し（Ｓ２２）、ＭＤのトレンドを算出する（Ｓ２３）。

図２２は、実施例３の複数系列の信号値に対して算出されたＭＤ及びＭＤのトレンドを示すグラフである。図２２の横軸は時間軸を示し、縦軸（左側）にＭＤに対応する数値を示している。このように、６８系列の信号値（時系列データ）から定義される各時点における“場”の観測値が単一系列の信号値として求められる。ＭＤトレンドの算出により、複数系列の信号値のＭＤから、各系列の信号値を発現させる“場”のトレンドが得られる。

単一系列の信号値として求められた後は、実施の形態１における演算過程と同様である。検出装置１の制御部１０は、算出されたＭＤトレンドに対して勾配を算出する（Ｓ２４）。

図２３は、実施例３の複数系列の信号値に対して算出されたＭＤトレンドの勾配を示すグラフである。図２３の横軸は時間軸を示し、縦軸（右側）にＭＤトレンドの勾配に対応する数値を示している。なお図２３には、比較のために図２２のＭＤ及びＭＤトレンドを上部に示し、縦軸（左側）にＭＤに対応する数値を示している。ＭＤトレンドの勾配の算出により、複数系列の信号値を発現させる“場”のトレンドの局所的な変化が得られる。

次に検出装置１の制御部１０は、算出されたＭＤトレンドの勾配のＭＤを算出する（Ｓ４）。

図２４は、実施例３の複数系列の信号値から得られたＭＤトレンドの勾配のＭＤを示すグラフである。図２４の横軸は時間軸を示し、縦軸（左側）にＭＤトレンドの勾配のＭＤに対応する数値を示している。なお図２４には、比較のために図２３のＭＤトレンドの勾配を下部に示し、縦軸（右側）にＭＤトレンドの勾配に対応する値を示している。ＭＤトレンドの勾配のＭＤを算出することにより、“場”のトレンドの局所的な変化の類似性を示す指標（値が大きいほど遠く、非類似を示す）が得られる。

検出装置１の制御部１０は、算出されたＭＤトレンドの勾配のＭＤを、最大値によって除算することによって変化点スコアを算出する（Ｓ５）。

図２５は、実施例３の複数系列の信号値に対して導出された変化点スコアを示すグラフである。図２５の横軸は時間軸を示し、縦軸（右側）に変化点スコアを示している。変化点スコアは最大値で除算されて算出されているから、変化点スコアの最大値は１．０である。図２５には比較のために図２４のＭＤトレンドの勾配のＭＤを上部に示し、縦軸（左側）にＭＤトレンドの勾配のＭＤに対応する値を示している。なお図２５に示す変化点スコアは、算出された変化点スコアの内、上位３０点のみ表示している。このように変化点スコアの最大値と比較して変化点を検出することが可能となる。図２５に示すように、２００７年５月から２００８年３月までの期間において、２０００年１月以降最大の変化点の群が検出されている。

図２６は、実施例３にて導出された変化点スコアをＭＤトレンドと比較して示すグラフである。図２６は、複数系列の信号値に対して算出されたＭＤ及びＭＤのトレンド（図２２）を上部に示し、変化点スコア（図２５）を下部に示して比較したものである。図２６の横軸は時間軸を示し、縦軸（右側）に変化点スコアを示し、縦軸（左側）にＭＤに対応する値を示している。図２６に示すように、２０００年１月以降最大の変化点の群が検出された後は、ＭＤのトレンドが１．０以上を示しており、明らかに“場”が変化していることを読み取ることが可能である。

図２７は、実施例３にて導出された変化点スコアと、複数系列の信号値に係る指標値とを比較するグラフである。図２７には、図２１に示した複数項目の月次データを基に導出された変化点スコア（図２２）を下部に示し、図２１の複数項目に基づく総合指標値である消費者物価指数を上部に示している。また、図２７には、実施の形態２の検出装置１によって導出された変化点スコア（ハッチングで示す棒グラフ）に加え、消費者物価指数を単一系列の時系列データＤ_t として実施の形態１の検出装置１により導出した変化点スコア（白抜きの棒グラフ）を共に示している。図２７の横軸は時間軸を示し、縦軸（左側）は消費者物価指数に対応する数値を示し、縦軸（右側）は変化点スコアに対応する値を示している。

図２７に示すように、２００８年前半及び後半の景気後退に対応する消費者物価指数の変化点の前に、２００７年５月から、消費者物価指数の構成項目である６８系列の信号値の変化点群が検出されている。これにより、２００７年の６８系列の項目の変化が、翌年の消費者物価指数の変化点に繋がっているという構造の変化が観測される。これにより、消費者物価指数の変化と、該指数に関連するデータの変化との関係性、時間特性等の把握が可能となり、例えば消費者物価指数の変化の予測に、複数項目の月次データにおける変化点を用いるなどの各種産業的な利用が可能となる。

このようにして、着目する対象に係る単一系列の時系列データに対する変化点スコアのみならず、前記着目対象に影響する複数系列の時系列データから導出される変化点スコアを加味することにより、着目対象自体のトレンド及びトレンドに影響する“場”を考慮した変化点を、外れ値の影響を低減させて検出することができる。着目対象に関するシステム全体の特性変化の検知、時間特性の定量化などが可能となる。

なお、上述のように開示された本実施の形態はすべての点で例示であって、制限的なものではないと考えられるべきである。本発明の範囲は、上記した意味ではなく、特許請求の範囲によって示され、特許請求の範囲と均等の意味および範囲内でのすべての変更が含まれることが意図される。

１ 検出装置
１０ 制御部
１１ 記録部
１３ 入力部
１Ｐ 検出プログラム

Claims (7)

1. 時系列的に連続する事象に対応する一連の値からなる系列を１又は複数入力し、１又は複数の前記系列を時系列的に記録する記録手段と、前記事象の記録期間中における変化点を検出する検出手段とを備える変化点検出装置において、
   前記検出手段は、
   前記記録手段に記録されている前記系列に基づく時系列データの記録期間中におけるトレンドをトレンドモデルに基づき算出するトレンド算出手段と、
   該トレンド算出手段が算出したトレンドの勾配を算出する勾配算出手段と、
   算出された勾配のマハラノビスの汎距離を算出する汎距離算出手段と
   を備え、
   前記検出手段は、算出されたマハラノビスの汎距離に基づき変化点を検出するようにしてある
   ことを特徴とする変化点検出装置。
2. 前記トレンド算出手段は、
   前記記録手段に記録されている複数系列に対してマハラノビスの汎距離を算出して前記時系列データとする時系列データ化手段を更に備え、
   算出された時系列データの記録期間中におけるトレンドをトレンドモデルに基づき算出するようにしてある
   ことを特徴とする請求項１に記載の変化点検出装置。
3. 前記トレンド算出手段は、
   前記事象に対応する一連の値が、該値を発現させる状態が前記値に与える影響を表現する関数と白色雑音との和で記述される状態空間モデルに、各時点における前記一連の値がトレンドと白色雑音との和で記述されるトレンドモデル、及び前後する時点におけるトレンドは略等しいとするモデルを適用し、状態空間モデルの解法に基づきトレンドを算出する
   ことを特徴とする請求項１又は２に記載の変化点検出装置。
4. 前記勾配算出手段は、
   前後する時点の内、前の時点におけるトレンドから後の時点におけるトレンドへの変化率を算出するようにしてある
   ことを特徴とする請求項１乃至３のいずれか１つに記載の変化点検出装置。
5. 前記検出手段は、
   前記汎距離算出手段により算出されたマハラノビスの汎距離を、前記記録期間中におけるマハラノビスの汎距離の最大値で除算した値を求める手段
   を更に備えることを特徴とする請求項１乃至４のいずれかに記載の変化点検出装置。
6. 時系列的に連続する事象に対応する１又は複数の一連の値からなる系列を入力し、１又は複数の系列を時系列的に記録する記録手段を備えるプロセッサが、前記事象の記録期間中における変化点を検出する変化点検出方法において、
   前記プロセッサは、
   前記記録手段に記録されている前記系列に基づく時系列データの記録期間中におけるトレンドをトレンドモデルに基づき算出し、
   算出されたトレンドの勾配を算出し、
   算出された勾配のマハラノビスの汎距離を算出し、
   算出されたマハラノビスの汎距離に基づき変化点を検出する
   ことを特徴とする変化点検出方法。
7. 時系列的に連続する事象に対応する１又は複数の一連の値からなる系列を入力し、１又は複数の系列を時系列的に記録する記録手段を備えるコンピュータに、前記事象の記録期間中における変化点を検出させるコンピュータプログラムにおいて、
   前記コンピュータに、
   前記記録手段に記録されている前記系列に基づく時系列データの記録期間中におけるトレンドをトレンドモデルに基づき算出するステップと、
   算出されたトレンドの勾配を算出するステップと、
   算出された勾配のマハラノビスの汎距離を算出するステップと、
   算出されたマハラノビスの汎距離に基づき変化点を検出するステップと
   を実行させることを特徴とするコンピュータプログラム。