【発明の詳細な説明】

【０００１】

【発明の属する技術分野】本発明は再帰型ニューラルネットワークに係り、時間的に不連続に変化する測定量の変化のトレンドを推定する推定システムおよびその方法に関する。

【０００２】

【従来の技術とその問題点】従来より予測フィルタとして用いられているカルマンフィルタは、システム同定方法の古典であり、現在も様々な応用が考案されている。
一方、非線形性を持つ時系列データの予測には、静的な非線形システムの同定法として技術的に確立されつつあるニューラルネットワークを適用する様々な方法が考案されている。 しかし、カルマンフィルタの応用には現状の計算機の能力などからくる制限があり、また、ニューラルネットワークの時系列解析に対する応用は始まったばかりである。 このため、それぞれ次のような問題点が指摘される。 （１）カルマンフィルタを用いる方法 時系列データに不連続なトレンド変化が生じる場合やモデルの不確定さを表現するノイズにガウス型を仮定できない場合は、通常の線形ガウス型のカルマンフィルタでは予測や濾波がうまくいかない（北川源四郎．時系列解析プログラミング．岩波コンンピュータサイエンス．岩波書店，1993）。 もし、不連続な状態変化に線形ガウス型モデルで対処するなら、極端に次元の大きなモデルが必要になる。 この際、モデル次元を客観的に選択するための基準を設定することが難しくなる。

【０００３】そこで、近年、非線形非ガウス型の拡張カルマンフィルタが提案されている（Genshiro Kitagawa.
Non-Gaussian State-Space Modeling of Nonstationar
y Time Series. Journal of the American Statistical
Association, 82(400):1032-1041, 1987.）。 拡張カルマンフィルタは、不連続なトレンドや非ガウス型のノイズに対する予測や平滑に成功している。 しかし、拡張カルマンフィルタを実行するには、予測、濾波、平滑の各分布を直接計算しなければならない。 したがって、状態空間モデルが大きくなる場合は、拡張カルマンフィルタの適切なフィルタ係数を同定するには時間がかかる。 これに対して、線形ガウス型モデルでは、平均値と分散を推定するだけで各段階で計算に必要な確率分布を決定できるので、同定に関する計算量は低く抑えることができるものの、同定可能対象は制限される。 また、拡張カルマンフィルタを効果的に応用するには、異常値なども含めたノイズを適当に表現できる分布族に関する先見的な知識を必要とする。

【０００４】しかし、サンプルからブートスラップ法を用いてノイズの分布を推定するモンテカルロフィルタを用いれば、ノイズに関する先見的な知識を欠く場合でも適切な推定や濾波および平滑ガ遂行される（Genshiro K
itagawa. A Monte Carlo Filtering and Smoothing Met
hod for Non-Gaussian Nonlinear State Space Models.
Research Memorandum 462, The Institute of Statist
ical Mathematics, 121993.）。 このモンテカルロフィルタの手法により、より一般的な非線形非ガウス型の時系列に対する方法論が確立しつつある。 ただし、リサンプリングによる確率分布の算定に必要な時間は無視できないほど大きくなる。 （２）ニューラルネットワークを用いる方法 時系列データから時間窓によリデータを切り出し、時間的に少しづつずれている一連のパターンを生成し、フィードフォワード型のニューラルネットワークとパックプロパゲーションによって時系列データを学習する方法がある（A. Waibel. Modular Construction of Time-Dela
y Networks for Speech Recongnition.Neural Computat
ion, 1:382-399, 1989. ／ Jeng-Neng Hwang, Shyh-R
ong Lay, Martin Maechler, R. Douglas Martin, and J
ames Schimert. Regression Modeling in Back-Propaga
tion and Projection Pursuit Learning. IEEE Transac
tions on Neural Networks, 5(3):342-353, May 199
4.）。 この方法で時系列データを的確に学習するためには、ニューラルネットワークの規模が大きくなり、記憶領域が不足するという問題が残る。 これは時間軸のデータの変化を空間軸に展開し、入出力の相関をニューラルネットワークの重みで表現するためである。 これにより別の問題も生じる。 時系列データの確率構造を捉えるという視点による明確な記述がなされないことである。

【０００５】こうしたニューラルネットワークの規模の問題を解決する方法として、フィードバック構造を持った再帰型ニューラルネットワークが考案されている（Je
romeT. Connor, R. Douglas Martin, and LE Atlas.
Recurrent Neural Networksand Robust Time Series P
rediction. IEEE Transactions on Neural Networks, 5
(2):240-254, Mar 1994. ）。 再帰型ニューラルネットワークには、主に２つの型がある。 出力層を回帰させる方法（Jordan型）と中間層を回帰させる方法（Elman
型）である。 再帰型ニューラルネットワークには回帰する情報を蓄える層が特別に設けられていて、これを文脈層と呼ぶ。

【０００６】フィードバック構造により規模の問題は解決するが、文脈層がどれくらいの規模と密度で過去の情報履歴を回帰させれば適切な予測フィルタを構成できるのかが不明確である。 また、未知の時系列データ（パラメータ同定に用いるデータと同じ確率構造から生成される別のデータ）に対する動作や評価方法などが不明確である。 一般的なネットワークの結合を考える場合は空間的な計算コストが莫大になることは明らかであり、各種係数を探索する場合に、微分係数などの算定のために必要な情報が長大になる。

【０００７】そこで、ＡＲＭＡ（autoregressive movin
g average ）モデルに類似した制約構造の再帰型ニューラルネットワークを用いて、予測フィルタを構成する方法が考えられている（James Ting-Ho Lo. Synthetic Ap
proach to Optimal Filtering. IEEE Transactions on
Neural Networks, 5(5):803-811, Sep 1994. ／ GV
Puskorius and LA Feldkamp. Recurrent Neural Ne
tworks with the Decoupled Extended Kalman Filter A
lgorithm. Science of Artifitial Neural Networks, 1
710:461-473, 1992. ）。 この場合、通常のカルマンフィルタによってニューラルネットワークの内部状態の解釈が与えられる。 二乗誤差を最小にする評価基準とカルマンフィルタに準拠する計算法を用いて、与えられたパラメータにおける誤差を算定しながら、与えられたデータに適するパラメータを選択する方法も考案されている。 しかしながら、この手法は、先に挙げたカルマンフィルタの計算量に関する問題を内在している。 また、ネットワークの内部状態と時系列データとの関係付けにおいて不明確な点が多く、内部状態の解釈が難しい。

【０００８】本発明は、再帰型ニューラルネットワークの内部状態と時系列データの関係を明確にし、不連続に変化する時系列のトレンドを効率的に推定する時系列トレンド推定システムとその方法を提供することを目的とする。

【０００９】

【問題点を解決するための手段】図１は、本発明の時系列トレンド推定システムの原理図である。 図１の時系列トレンド推定システムは、入力手段１、ニューラルネットワーク手段２、予測値生成手段３、および出力手段４
を備える。

【００１０】入力手段１は、時系列データを入力する。
ニューラルネットワーク手段２は、文脈層を有するコラム構造の再帰型ニューラルネットワークを含み、過去の時系列データに関する情報を含む内部状態を出力する。

【００１１】予測値生成手段３は、上記内部状態を用いて予測値の候補の出現確率を求め、最も確率の高い候補を予測値として求める。 出力手段４は、上記予測値を未知データの推定結果として出力する。

【００１２】図２は、図１の時系列トレンド推定システムによる推定処理のフローチャートである。 図２のステップＳ１において、まず、入力手段１は現在の時刻のデータを入力する。

【００１３】ステップＳ２において、ニューラルネットワーク手段２は過去のデータ情報を持つ再帰型ニューラルネットワークの内部状態をコラム毎に設定する。 再帰型ニューラルネットワークは複数のコラムを備え、各コラムは入力データから出力を生成するニューラル素子と、ニューラル素子の過去の出力を保持する文脈層とを含む。 現在の時刻のデータが入力されると、ニューラル素子は入力データと同じコラムの文脈層からの回帰データとから新しい出力を生成し、文脈層のデータ情報は１
時刻分シフトされる。

【００１４】ステップＳ３において、予測値生成手段３
は内部状態により決められる予測方程式を解いて、予測値の候補を求める。 このとき、予測値生成手段３はニューラルネットワーク手段２から受け取った内部状態を用いて、予測方程式の解を求める。

【００１５】ステップＳ４において、予測値生成手段３
は内部状態を用いて予測値の候補の出現確率を求める。
ステップＳ５において、出力手段４は出現確率の最も大きなものを未知データの予測値として出力する。

【００１６】各コラムのニューラル素子に、時系列データの不連続なトレンドの離散性を表現する離散変数を割り付けることにより、離散変数と内部状態により一種のエネルギー関数を構成することができる。 このエネルギー関数から定義される確率分布関数を用いて、上記予測値の候補の出現確率を表す。 後述するように、こうして定義された確率分布関数に関する考察から、各コラムの内部状態は対応するニューラル素子の離散変数が１となる確率を与えていることが分かる。 また、ある離散変数が１となる確率が高ければ、そのコラムのパラメータにより決められる平均値を持つ確率分布の寄与が大きくなる。 したがって、内部状態は、データの特定の確率分布が選択される確率に関する情報を表しているといえる。
このようにして構成された時系列トレンド推定システムにおいては、ニューラルネットワークの内部状態を時系列データの離散値の出現確率と関連づけて解釈することが可能となる。 また、各ニューラル素子には同じコラム内のデータ情報のみが再帰的に入力されるので、各コラムの独立性が高い。 したがって、各コラムのパラメータを、不連続な変化トレンドの離散値に適合するように調整することにより、そのトレンドの推定が容易になる。

【００１７】図１のニューラルネットワーク手段２は、
実施形態の図３におけるコラム構造再帰型ニューラルネットワーク１９に対応し、予測値生成手段３は非線形方程式求解装置１８に対応する。 また、入力手段１および出力手段４は、表示・対話装置１２に対応する。

【００１８】

【発明の実施の形態】以下、図面を参照しながら本発明の実施の形態を詳細に説明する。 本実施形態においては、次のような方針でシステムを構築する。 （ａ）線形モデルでは追従できないトレンドの急激な変化に追従する能力を、ニューラルネットワークの隠れ素子に離散変数を割り付けることによって保証する。 この際、適切な追従能力を持ったフィルタを構成するために必要な同定法も与える。 （ｂ）線形モデルや単純な分布を想定することでは対処できない異常値に対して、ネットワークの非線形特性の一つである飽和特性を利用して対処する。 （ｃ）ガウス型の加重和によって、予測機構およびフィルタとしての測定量の確率密度関数を合成する。 確率密度関数は、ニューラルネットワークの内部状態に基づいて計算する。 これにより、ネットワークの内部状態と時系列との関係付けを明確にできる。 また、小規模のネットワークであれば効率的に実行できる計算法を用いる。 （ｄ）ネットワークの構造は Elman型の構成とし、文脈層にコラム構造を導入してその構造を制約する。 これにより、回帰情報の独立性を保証し、最適な係数の探索に必要となる内部状態の再計算を容易にする。

【００１９】図３は、本発明の時系列トレンド推定システムの構成図である。 図３の時系列トレンド推定システムは、制御装置１１、学習装置１４、および予測装置１
７を備える。 制御装置１１は表示・対話装置１２とシステム管理装置１３からなり、学習装置１４は観測装置１
５とネットワーク学習装置１６からなる。 また、予測装置１７は、非線形方程式求解装置１８とコラム構造再帰型ニューラルネットワーク１９からなる。 以後、コラム構造再帰型ニューラルネットワークをＣＳＳＲＮＮ（Co
lumn-Structured Simple Reccurent Neural Network ）
と記す。

【００２０】この時系列トレンド推定システムは、例えば図４に示すような計算機システムとして実施される。
図４の計算機システムは、ＣＰＵ（中央処理装置）２
１、メモリ２２、入出力端末２３、およびそれらを接続するバス２４を備える。 入出力端末２３は、例えばディスプレイ装置やキーボードを有する端末装置であり、表示・対話装置１２に対応する。 また、バス２４には、必要に応じてハードディスク等の外部記憶装置やプリンタ等が接続される。 システム管理装置１３、学習装置１
４、および予測装置１７の各機能は、ＣＰＵ２１がメモリ２２に格納されたプログラムを実行することにより実現される。

【００２１】制御装置１１は、時系列トレンド推定システムと他のシステムまたはユーザ（以下、クライアントと呼ぶ）との間のデータの授受に必要なインタフェースである。 また、必要な計算機資源の確保およびデータの流れの制御を行う。 図３において、実線の矢印はデータの流れを表し、破線の矢印は制御情報の経路を表している。 時系列トレンド推定システムの動作には、学習モードと予測モードがある。 これらの２つのモードは制御装置１１のシステム管理装置１３により制御される。 システム管理装置１３は、表示・対話装置１２からの情報に基づいて、モードの切替を行なう。 学習モードでは、学習装置１４が起動される。 このときデータの入力ラインは制御装置１１から学習装置１４に接続される。 必要な長さの時系列データが集まると、それをもとに学習が遂行される。 このとき、予測装置１７は、ＣＳＳＲＮＮ１
９の内部状態の更新に必要な入力を学習装置１４から得る。 予測モードでは、学習装置１４は起動待ち状態になり、入力ラインは直接予測装置１７に接続される。 予測装置１７は、最適なパラメータセットをもとに予測を遂行する。 その予測結果は、システム管理装置１３を経由して、表示・対話装置１２によりクライアントに提示される。

【００２２】学習装置１４の観測装置１５は、観測データｘ ^(t)を時間の経過に従って読み込む。 次に、学習装置１４はネットワーク学習装置１６を起動し、読み込んだデータを時系列データ｛ｘ ₁ ，． ． ． ，ｘ _N ｝としてネットワーク学習装置１６に渡す。 ネットワーク学習装置１６は最適なパラメータセットΦ _optを学習によって獲得し、予測装置１７のＣＳＳＲＮＮ１９に出力する。
その後、学習装置１４は起動待ち状態となる。

【００２３】ネットワーク学習装置１６は、時系列データ｛ｘ ₁ ，． ． ． ，ｘ _N-1 ｝とパラメータセットΦを予測装置１７に渡す。 ここで、予測装置１７に最初に渡すΦは、ランダムなパラメータセットである。 そして、予測装置１７から予測値 外１（以後、｛ｘ ₁ハット，．．．，ｘ _Nハット｝と記す。）を受け取り、時系列

【００２４】

【外１】

【００２５】データとＣＳＳＲＮＮ１９によって表現されたモデルとの食い違いを評価する評価関数を計算しながら、その値に基づいてパラメータセットΦを修正する。 この操作は、評価関数が最小になるような最適なパラメータセットが見つかるまで繰り返される。

【００２６】予測装置１７は、与えられたパラメータセットΦをもとにＣＳＳＲＮＮ１９を駆動して、ＣＳＳＲ
ＮＮ１９の内部状態を更新し、内部状態の情報をもとに非線形方程式求解装置１８を使って時系列データの予測を行なう。

【００２７】ＣＳＳＲＮＮ１９は、入力層、隠れ素子層、および文脈層からなる再帰型ニューラルネットワークである。 文脈層は、ＣＳＳＲＮＮ１９の各隠れ素子の過去の出力を保存するシフトレジスタからなる。 各隠れ素子の過去の出力履歴は、それ自身にだけ回帰するのであって、その他の素子に直接伝達されることはない。 このＣＳＳＲＮＮ１９の文脈層のシフトレジスタ内部の値全体を内部状態と呼ぶ。 ＣＳＳＲＮＮ１９は、早い時刻のデータから順に時間軸に沿って入力｛ｘ ₁ ，． ． ． ，
ｘ _N-1 ｝を受け取りながら、逐次内部状態を更新する。

【００２８】非線形方程式求解装置１８は、予測のための非線形方程式を解いて、データの予測値 外２ （以後、ｘ ^(t)ハットと記す）を計算する。 非線形方程式の係数

【００２９】

【外２】

【００３０】は、ＣＳＳＲＮＮ１９の内部状態から決められる。 次に、各装置の構成要素と動作、およびそれらの間のデータの流れを説明する。

【００３１】表示・対話装置１２は、クライアントに対して時系列データおよび推定したトレンドを提示する。
時系列データおよび堆定トレンドは、予測装置１７からシステム管理装置１３を通じて表示・対話装置１２に提供される。 クライアントは、提示された結果に不満があれば、予測装置１７のパラメータの再学習を要求できる。 その際、クライアントは、ＣＳＳＲＮＮ１９の規模を指定することができる。 また、学習装置１４の学習アルゴリズムも変更することができる。 学習アルゴリズムの変更の際には、学習基準の変更や学習アルゴリズムの各パラメータの指定・変更が、表示・対話装置１２を通じて行われる。 クライアントによる再計算の要求や各パラメータの変更指示は、システム管理装置１３に伝達される。

【００３２】システム管理装置１３は、学習装置１４および予測装置１７に必要な計算機資源を確保して、各装置を起動する。 システム管理装置１３は、図３の破線で示された制御経路を通じて、学習装置１４と予測装置１
７の動作を制御する。 学習モードでは、学習装置１４を呼び出して、予測装置１７のパラメータの同定を遂行させる。 この際、システム管理装置１３が直接予測装置１
７と通信することはない。 予測モードでは、システム管理装置１３は予測装置１７と直接通信して、時系列データを伝達し、トレンド予測値を受け取る。 トレンド予測値は、表示・対話装置１２を通じてクライアントに提供され、システム管理装置１３はクライアントからのフィードバックを受けとる。 クライアントが予測結果に満足しない場合は、その要求に応じて学習装置１４を起動し、予測装置１７のパラメータの再学習を開始させる。
この際、クライアントの要求事項として、学習基準の変更やＣＳＳＲＮＮ１９の規模の変更があれば、そのために必要な計算機資源を確保し、パラメータの初期値を設定して、学習装置１４を呼びだす。 学習装置１４によって更新されたパラメータにより、予測装置１７は、再びＣＳＳＲＮＮ１９を駆動して予測を行なう。 新しく得られた予測データは、表示・対話装置１２へ渡され、クライアントに提示される。

【００３３】学習装置１４の主要な機能は、ネットワーク学習装置１６によって提供される。 観測装置１５は、
必要な個数の標本を連続して収集し、パラメータの学習に必要な時系列データＳ＝｛ｘ ₁ ，． ． ． ，ｘ _N ｝を生成するとき以外は、起動待ち状態にある。

【００３４】図５は、ネットワーク学習装置１６の構成図である。 図５のネットワーク学習装置１６は、シミュレーティド・アニーリング制御装置３１、ランダムシンプレックスによるパラメータの精錬装置３２、および評価関数算定装置３３からなる。 ネットワーク学習装置１
６は、予測装置１７の適切な動作を保証するパラメータとして、例えば、評価関数

【００３５】

【数１】

【００３６】が最小となるΦ＝Φ _optを学習によって求める。 （１）式において、ｘ _iハットは予測装置１７から受け取る予測値であり、そのときにＣＳＳＲＮＮ１９
に与えられたΦに依存している。 ネットワーク学習装置１６は、評価関数算定装置３３により予測装置１７と通信しながら学習を遂行する。

【００３７】評価関数算定装置３３は、時系列データ｛ｘ ₁ ，． ． ． ，ｘ _N ｝と予測値列｛ｘ ₁ハット，． ． ． ，ｘ _Nハット｝とを入力として、（１）式の評価関数 merit（Ｓ｜Φ）の値を計算し、出力する。
｛ｘ ₁ハット，． ． ． ，ｘ _Nハット｝は、Φおよび｛ｘ
₁ ，． ． ． ，ｘ _N ｝を予測装置１７へ入力した時にその出力として得られる。

【００３８】アニーリング制御装置３１は、パラメータの選択基準、アニーリングの温度制御ルール、ランダム探索の回数の上限値Ｎ _r 、およびアニーリングの試行回数Ｎ _aを、制御データとしてシステム管理装置１３から受け取る。 そして、それらに基づいてシミュレーティド・アニーリングによる最適パラメータの探索処理を制御する。 アニーリング制御装置３１は、Ｎ _a個の初期値のパラメータセットΦ ^{(i )} _init （ｉ＝１，．．．，Ｎ _a ）
を生成し、ランダムシンプレックスによるパラメータの精錬装置３２に与える。

【００３９】パラメータ精錬装置３２は、必要に応じて、パラメータセットΦを入力として評価関数算定装置３３を起動し、 merit（Ｓ｜Φ）の値を受け取る。 パラメータ精錬装置３２は、公知の滑降シンプレックス法による局所最適化法を用いてパラメータセットΦを精練する。 各アニーリング試行において初期値のパラメータセットΦ ⁽ⁱ⁾ _initから、 dim（Φ）＋１個の点からなる初期値シンプレックスを生成して、探索を開始する。 ここで、 dim（Φ）はパラメータセットの独立変数の個数を表す。

【００４０】図６は、 dim（Φ）＝２の場合の初期値シンプレックスの例を示している。 図６において、パラメータセットΦは２つの独立変数φ１、φ２を用いてΦ＝
（φ１，φ２）と表され、φ１φ２平面上の任意の１点Φ ⁽ⁱ⁾ _initが初期値として与えられる。 パラメータ精錬装置３２は、点Φ ⁽ⁱ⁾ _initをもとに頂点Φ ₁ 、Φ ₂ 、Φ
₃を持つ三角形の初期値シンプレックス３４を生成する。 dim（Φ）＝３の場合には、初期値シンプレックスは４つの頂点を持つ３次元図形となる。

【００４１】各アニーリング試行において、パラメータ精錬装置３２は、シンプレックスの頂点の中で評価関数の値が最悪（最大）となる点、２番目に悪い点、および最良（最小）となる点を調べる。 次に、最悪値の点を除いた頂点の重心に対して、最悪値の点を対称移動して対称点（反点）を求め、その反点を最悪点の代わりに採用して新しいシンプレックスを生成する。 もし、反点において、評価関数の値が最良点の値よりも良ければ、重心から更にその方向に２倍の距離だけ進んだ点を新しい頂点とする。 また、反点の値が最悪値の次に悪い値よりも悪ければ、最悪点を重心方向へ１次元的に移動させ、新しい頂点を生成する。 それでも最良点よりも良い点が見つからなければ、最良点の方向へシンプレックス全体を収縮させる。 以上の操作を繰り返すことで、シンプレックスの形を変動させながら評価関数の値を小さくする方向を発見し、最終的にシンプレックスの大きさを縮小させることにより極小となる点を求める。

【００４２】このように、滑降シンプレックス法においては、シンプレックスが更新される毎に評価関数を計算する必要がある。 この計算は、既に説明したように、評価関数算定装置３３が予測装置１７と通信することによって実行される。

【００４３】図６においては、 merit（Ｓ｜Φ ₂ ）＞ m
erit（Ｓ｜Φ ₁ ）＞ merit（Ｓ｜Φ ₃ ）であるから、点Φ ₂が最悪点となる。 そこで、パラメータ精錬装置３２
は、初期値シンプレックス３４から点Φ ₂を除いた残りの頂点Φ ₁ 、Φ ₃の重心（この場合はΦ ₁とΦ ₃を結ぶ線分の中点）を求める。 そして、求めた重心に関して点Φ ₂と対称な点をφ１φ２平面内で求め、その点での m
erit（Ｓ｜Φ）を計算して、 merit（Ｓ｜Φ ₁ ）、 mer
it（Ｓ｜Φ ₃ ）と比較する。 もし、対称点での値が mer
it（Ｓ｜Φ ₃ ）より小さければ、重心からその方向にさらに離れた点を新しい頂点として、次のシンプレックスを生成する。

【００４４】また、パラメータ精錬装置３２は、与えられた制御温度Ｔ ^(j) （ｊ＝１，２，．．．，Ｍ）によって決まる偏差を使って、シンプレックスを探索空間においてランダムウォーク（酔歩）させ、評価関数の大域的最小値の近傍を探索する。 各制御湿度における酔歩の回数はＮ _rである。 この算法によって、各アニーリング試行から評価関数の大域的最小値を与えるパラメータセットの候補Φ ⁽ⁱ⁾ _res （ｉ＝１，．．．，Ｎ _a ）が求まる。

【００４５】アニーリング制御装置３１は、これらの候補の中で評価関数の値が最小となるパラメータセットを、最適なパラメータセットΦ _optとして選ぶ。 ネットワーク学習装置１６は、予測装置１７に最適なパラメータセットΦ _optを渡して、起動待ち状態になる。

【００４６】図７は、予測装置１７の構成図である。 図１７において、ＣＳＳＲＮＮ１９の隠れ素子層は、ｍ個の隠れ素子５１−１、５１−２、・・・、５１−ｍからなる。 各隠れ素子は、図８に示すニューロン（ニューラル素子）を有する。 任意の隠れ素子と隠れ素子の間に相互結合は存在しない。 各隠れ素子には、文脈層を構成するｓ個のレジスタがそれぞれ割り付けられている。 例えば、隠れ素子５１−１の出力側にはレジスタ５２−１−
１、５２−１−２、・・・、５２−１−ｓが設けられ、
隠れ素子５１−２の出力側にはレジスタ５２−２−１、
５２−２−２、・・・、５２−２−ｓが設けられ、隠れ素子５１−ｍの出力側にはレジスタ５２−ｍ−１、５２
−ｍ−２、・・・、５２−ｍ−ｓが設けられる。 １つの隠れ素子とそれに付随するレジスタによって、１つのコラムが形成される。 時刻ｔにおける隠れ素子に対する入力をｘ ^(t)とすると、同時刻におけるｊ番目の隠れ素子（ｊ＝１，．．．，ｍ）の出力は、

【００４７】

【数２】

【００４８】となる。 表記を簡潔に行なうため、幾つかの簡略表現を次のように導入する。

【００４９】

【数３】

【００５０】ここで、関数ＬＧ（ｘ）は、一般にロジスティック関数と呼ばれており、図８のニューロンの出力関数に相当する。 また、ｗ _jは入力ｘ ^(t)に対するｊ番目の隠れ素子の重みである。 Ｏ _j ^(ti)は時刻ｔ−ｉにおけるｊ番目の隠れ素子の出力であり、 外３ （以後、ベクトルＯ _j ^(t,s)と記す）は、ｊ番目のコラムの

【００５１】

【外３】

【００５２】ｓ個のレジスタの値を成分とするベクトルである。 また、ｕ _jiはｉ番目のレジスタからの再帰値に対するｊ番目の隠れ素子の入力重みであり、 外４
（以後、

【００５３】

【外４】

【００５４】ベクトルｕ _jと記す）は、それらの入力重みを成分とするベクトルである。 θ _jはｊ番目の隠れ素子のロジスティック関数のバイアスであり、一般的にはしきい値と呼ばれている。 以上の表記法を用いて、ＣＳ
ＳＲＮＮ１９のパラメータセットを表記するとΦ＝
｛（ｗ _j ，θ _j ，ベクトルｕ _j ）：（ｊ＝１，．．．，
ｍ）｝となる。 この場合、パラメータセットの次元 dim
（Φ）は（ｓ＋２）ｍとなる。 また、

【００５５】

【数４】

【００５６】で定義される 外５ （以後、Ｏ ^(t,s)バーと記す）を、時刻ｔにおけるＣＳＳ

【００５７】

【外５】

【００５８】ＲＮＮ１９の内部状態と呼ぶことにする。
内部状態Ｏ ^(t,s)バーの更新は、各コラム毎に行われる。 まず、（２）式により各隠れ素子の時刻ｔの出力Ｏ
_j ^(t) （ｊ＝１，．．．，ｍ）を計算する。 次に、レジスタ内の値Ｏ _j ^(t-1) ，． ． ． ，Ｏ _j ^(ts)をシフトして、Ｏ _j ^(t)を先頭のレジスタ５２−ｊ−１に保存する。 この際、レジスタ５２−ｊ−ｓ内にある最も古い時刻の出力値Ｏ _j ^(ts)が廃棄される。 そして、各コラムのレジスタ内の値はＯ _j ^(t-1) ，． ． ． ，Ｏ _j ^(ts)からＯ _j ^(t) ，． ． ． ，Ｏ _j ^(t-s+1)に更新される。 内部状態の更新は、各コラムにおいて時刻ｔ＋１までに終了する。

【００５９】このようなコラム構造を用いることにより、ＣＳＳＲＮＮ１９の各々の隠れ素子に、その隠れ素子自身の過去における出力履歴を時間の序列を保持したままで再帰させることができる。 したがって、各コラムの再帰情報の独立性が保証される。

【００６０】こうして、時刻ｔ＝１，２，． ． ． 毎に観測値ｘ ^(t)が新しく入力されると、ＣＳＳＲＮＮ１９は状態をＯ ^(t,s)バーからＯ ^(t+1,s)バーに更新する。 このネットワークの内部状態Ｏ ^(t,s)バーから、時刻ｔにおけるｘの予測密度関数および時刻ｔにおけるｘの予測値（時刻ｔ＋１の入力に対する予測値）が計算できる。
予測密度関数は、ネットワークの内部状態によって定義されたエネルギー関数から導くことができる。 このエネルギー関数は、各隠れ素子に対して隠れ変数ｈ _j ∈
｛０，１｝（ｊ＝１，．．．，ｍ）を導入し、隠れ素子の出力を隠れ変数の値が１となる確率と解釈することにより定式化できる。 ｈ _jは０または１の離散値をとるため、離散変数と呼ぶこともできる。 ｘ ^(t)に対する確率密度を表す予測確率密度関数を、隠れ変数ｈ _jを用いて表すと、

【００６１】

【数５】

【００６２】となる。 ここで、Ｚ ^(t)は、

【００６３】

【数６】

【００６４】で与えられる正規化係数であり時間と共に変化する。 （３）式の予測確率密度関数は、複数のガウス確率密度関数を混合したものである。 例えば、隠れ素子１個からなるＣＳＳＲＮＮ１９（ｍ＝１）を考えると、予測確率密度関数は、（３）式より、

【００６５】

【数７】

【００６６】となる。 （５）式の予測確率密度関数は、
２個のガウス密度関数の合成関数である。 この場合、ｘ
^(t) ＝０にピークを持つ１番目のガウス密度関数に対する重みは１／Ｚ ^(t)であり、ｘ ^(t) ＝ｗ ₁にピークを持つ２番目のガウス密度関数に対する重みは、

【００６７】

【数８】

【００６８】となる。 同様にして、ｍ個の隠れ素子からなるＣＳＳＲＮＮ１９が表現する確率密度関数は２ ^m個のガウス密度関数を合成したものになる。 各ガウス密度関数に対する合成の重みは、

【００６９】

【数９】

【００７０】の関数となる。 一方、Ｏ _j ^(tk)は、
（２）式およびベクトルＯ _j ^(t,s)の定義式より、入力ｘ ^(t-1) ，． ． ． ，ｘ ⁽¹⁾の値に依存しているので、合成の重みは入力ｘ ^(t-1) ，． ． ． ，ｘ ⁽¹⁾の変化によって変化する。 これは、事前の入力により予測確率密度関数の形状が変化することを意味しており、このことを用いて不連続トレンドの推定が可能になる。

【００７１】予測確率密度関数の微分から、トレンドの予測に必要な非線形方程式が得られ、

【００７２】

【数１０】

【００７３】となる。 （３）式および（６）式の導出方法については、後に詳しく説明する。 図７の非線形方程式求解装置１８は、予測確率密度分布および予測値の計算を遂行する装置である。 非線形方程式求解装置１８
は、目標関数の零点を求める零点算定装置４１、予測値選択装置４２、目標関数評価装置４３、正規化係数算定装置４４、確率密度算定装置４５からなる。

【００７４】（６）式から、目標関数ＴＧ（ｘ）を、

【００７５】

【数１１】

【００７６】とおくことができる。 ＴＧ（ｘ）＝０となるｘの値が（６）式の解となり、（３）式の予測確率密度の極大値を与える値に対応する。 零点算定装置４１
は、区間分割と単純な囲い込み法を繰り返して、目標関数の零点を与える変数ｘの値ｘ ¹ _peak ，． ． ． ，ｘ ^k
_peakを近似的に同定する。 このとき、必要に応じて変数値ｘを目標関数評価装置４３に与えて、（７）式のＴＧ
（ｘ）の値を計算させる。

【００７７】予測値選択装置４２は、予測確率密度関数の極大値を与える変数値ｘ ¹ _peak ，． ． ． ，ｘ ^k _peakの中で最大の極大値を与える変数値を選択する。 このとき、予測値選択装置４２は、確率密度算定装置４５と通信することで予測値の候補ｘ ¹ _peak ，． ． ． ，ｘ ^k _peak
に対する確率密度関数値ｖ ¹ _peak ，． ． ． ，ｖ ^k _peakを得る。 これらの値から高速なソーティング法を用いて、
確率密度関数の最大値を与えるｘ _maxを求める。 予測値の候補が少ない場合は、確率密度関数値の単純な比較によって選択することが可能である。 また、確率密度が２
つ以上の異なる変数値において最大値を取る場合は、それらの変数値の上で一様な確率分布を用いて、確率的にどれか１つを選択し、ｘ _maxとする。 こうして選択された変数値ｘ _maxが、次の入力の予測値ｘ ^(t)ハットとして出力される。

【００７８】目標関数評価装置４３、正規化係数算定装置４４、確率密度算定装置４５の３つの装置は、その構成および動作がＣＳＳＲＮＮ１９の内部状態に深く根ざしている。 これらの各装置について、図９から図１１までを参照しながら詳しく説明する。

【００７９】図９は、目標関数評価装置４３の構成を示している。 目標関数評価装置４３は、ＣＳＳＲＮＮ１９
からの｛ベクトルｕ _j ，ベクトルＯ _j ^(t,s) ｝（ｊ＝
１，． ． ． ，ｍ）と零点算定装置からの変数値ｘを入力として処理し、ＴＧ（ｘ）を計算する。 そのために、目標関数評価装置４３は、ＣＳＳＲＮＮ１９に対応するコラム構造を持つ内部状態情報変換器６１、および加算器（Σ）６４、６５を備える。 内部状態情報変換器６１
は、加算器６２−１、６２−２、・・・、６２−ｍとロジスティック関数演算器（Logistic）６３−１、６３−
２、・・・、６３−ｍを備え、各加算器６２−ｊとロジスティック関数演算器６３−ｊがｊ番目のコラムを形成している。 図９において、矢印で表した各データ経路に付加された変数および定数は、その値を経路上のデータに乗算することを意味している。 記号の付加されていない経路を伝播するデータの値は、伝播の途中で変化することはない。

【００８０】各加算器６２−ｊは、あらかじめセットされた重み係数ｗ _jを入力ｘに乗じ、その結果にしきい値θ _jとベクトルの内積（ベクトルｕ _j・ベクトルＯ _j
^(t,s) ）を加えて、ロジスティック関数演算器６３−ｊ
に入力する。 ロジスティック関数演算器６３−ｊは入力された値ｙ _jからＬＧ（ｙ _j ）を計算し、内部状態情報変換器６１は、さらにそれらの出力ＬＧ（ｙ _j ）（ｊ＝
１，． ． ． ，ｍ）にｗ _jを乗じた値を出力とする。 次に、加算器６４は、内部状態情報変換器６１のｍ個の出力の総和を求める。 加算器６５は、零点算定装置４１から入力された変数値ｘから加算器６４の出力を差し引いて、その結果を目標関数評価装置４３の出力とする。

【００８１】図１０は、正規化係数算定装置４４の構成を示している。 正規化係数算定装置４４は、（４）式のＺ ^(t)を計算する装置であり、ＣＳＳＲＮＮ１９に対応するコラム構造を持つ内部状態情報変換器７１、２値ベクトル生成器７３、ベクトル乗算器７４、ベクトル加算器７５、ノルム算定器７６、加算器７７、指数関数計算器（ｅｘｐ）７８、累積型加算器７９、逆数演算器（ｒ
ｅｃｉｐ）８０を有する。 ２値ベクトル生成器７３、ベクトル乗算器７４、およびベクトル加算器７５は、ＣＳ
ＳＲＮＮ１９の隠れ素子に割り付けた隠れ変数ｈ ₁ ，． ． ． ，ｈ _mに関係する計算を遂行する。 ２値ベクトル生成器７３は、０または１の値をとるｍ個の隠れ変数ｈ ₁ ，． ． ． ，ｈ _mを１つの隠れ変数ベクトル 外６ （以後、ベクトルｈと記す）の成分として、２ ^m個の可能なベクトルｈ _i （ｉ＝１，．．．，

【００８２】

【外６】

【００８３】２ ^m ）を全て生成する。 ２値ベクトル生成器７３は、隠れ変数ベクトルを生成することにより、
（４）式を計算する一連の操作の流れを統括する。 内部状態情報変換器７１は、加算器７２−１、７２−２、・
・・、７２−ｍを備え、各加算器７２−ｊがｊ番目のコラムを形成している。 そして、各加算器７２−ｊがしきい値θ _jと内積（ベクトルｕ _j・ベクトルＯ _j ^(t,s) ）
を加えて出力する。 内部状態情報変換器７１の出力は、
ベクトル（ベクトルｕ ₁・ベクトルＯ ₁ ^(t,s) ＋
θ ₁ ，． ． ． ，ベクトルｕ _m・ベクトルＯ _m ^(t,s) ＋θ
_m ）となる。 正規化係数算定装置４４の１回の呼び出しについて、内部状態情報変換器７１は１度だけ起動され、ＣＳＳＲＮＮ１９の内部状態に関する情報を出力する。 その出力結果は、ベクトル乗算器７４の一方の入力として保存される。 ベクトル乗算器７４は、隠れ変数ベクトルｈ _i （ｉ＝１，．．．，２ ^m ）とＣＳＳＲＮＮ１
９の内部状態情報のベクトルを入力として、それらの内積

【００８４】

【数１２】

【００８５】を計算し、出力する。 ここで、各隠れ変数ベクトルをベクトルｈ _i ＝（ｈ _i1 ，．．．，
ｈ _ij ，． ． ． ，ｈ _im ）と表記する。 ベクトル加算器７５
は、ＣＳＳＲＮＮ１９の入力重み係数ｗ ₁ ，． ． ． ，ｗ
_mと各隠れ変数ベクトルｈ _i （ｉ＝１，．．．，２ ^m ）
から

【００８６】

【数１３】

【００８７】を計算する。 ノルム算定器７６は、（８）
式のノルム（絶対値）の２乗を算定し、その値に１／２
を乗じて、

【００８８】

【数１４】

【００８９】の値が計算される。 その結果は、加算器７
７によりベクトル乗算器７４の出力値に加算され、指数関数計算器７８に入力される。 指数関数計算器７８は、

【００９０】

【数１５】

【００９１】の値を計算し、その結果を累積型加算器７
９に渡す。 以上の操作を２ ^m個の各隠れ変数ベクトルｈ
_iについて繰り返すことで、（４）式の右辺のベクトルｈ ^(t)に関する総和が計算される。 累積型加算器７９
は、全ての隠れ変数ベクトルについての総和を求め、逆数演算器８０は、累積型加算器７９の出力に（２π）
^1/2を乗じて、その逆数を求める。 こうして、正規化係数算定装置４４からは（４）式のＺ ^(t)の逆数が最終的に出力される。

【００９２】図１１は、確率密度算定装置４５の構成を示している。 時刻ｔにおいて、確率密度算定装置４５
は、予測値選択装置４２から予測値の候補ｘ＝
ｘ ¹ _peak ，． ． ． ，ｘ ^k _peakを受け取り、これらの値に対する時刻ｔの予測確率密度関数の値ｖ ¹ _peak ，． ． ． ，ｖ ^k _peakを出力する。 時刻ｔにおける予測確率密度関数は（３）式で与えられ、確率密度算定装置４５は、ＣＳＳＲＮＮ１９に対応するコラム構造を持つ内部状態情報変換器８１、乗算器（П）８４、指数関数計算器８５、およびノルム算定器８６を有する。
内部状態情報変換器８１は、加算器８２−１、８２−
２、・・・、８２−ｍと指数関数計算器（１＋ｅｘｐ）
８３−１、８３−２、・・・、８３−ｍを備え、各加算器８２−ｊと指数関数計算器８３−ｊがｊ番目のコラムを形成している。

【００９３】各加算器８２−ｊは、予測値選択装置４２
からの入力値ｘ、ＣＳＳＲＮＮ１９の入力重み係数ｗ _j 、しきい値θ _j 、およびＣＳＳＲＮＮ１９のｊ番目のコラムの内部情報である内積値（ベクトルｕ _j・ベクトルＯ _j ^(t,s) ）から、ｗ _j ｘ＋ベクトルｕ _j・ベクトルＯ _j ^(t,s) ＋θ _jを求める。 また、各指数関数計算器８３−ｊは、加算器８２−ｊの出力から１＋ｅｘｐ（ｗ
_j ｘ＋ベクトルｕ _j・ベクトルＯ _j ^(t,s) ＋θ _j ）の値を計算する。 ノルム算定器８６は入力値ｘのノルムの２
乗を計算し、指数関数計算器８５は、その値を用いて入力ｘの分布Ｎ（０，１）の密度関数における値ｅｘｐ
（−｜ｘ｜ ² ／２）を求める。 乗算器８４は、内部状態情報変換器８１のｍ個の指数関数計算器の出力を積算し、さらに正規化係数算定装置から受け取る１／Ｚ ^(t)
と指数関数計算器８５から受け取る正規分布密度関数値ｅｘｐ（−｜ｘ｜ ² ／２）を積算する。 ただし、Ｎ
（０，１）は平均値０、分散１の正規分布を表す。 このようにして、（３）式のｘ ^(t)にｘ ¹ _{pe ak} ，． ． ． ，ｘ
^k _peakをそれぞれ代入した時の各ｐ（ｘ ^(t) ｜Ф，Ｏ
^(t,s)バー）の値ｖ ¹ _peak ，． ． ． ，ｖ ^k _peakが計算され、予測値選択装置４２に出力される。

【００９４】以上の実施形態において、評価関数算定装置３３が計算する評価関数 merit（Ｓ｜Φ）は、必ずしも（１）式の形式の平均２乗誤差である必要はない。 例えば、時系列データとＣＳＳＲＮＮの内部状態により決まる予測確率密度関数の対数尤度であってもかまわない。 その場合、評価関数は、

【００９５】

【数１６】

【００９６】と書ける。 ランダムシンプレックス法とシミュレーティド・アニーリング法を用いる図５のネットワーク学習装置１６において、評価関数を（９）式のように置き換えたとしても、学習アルゴリズムの基本動作は変わらない。

【００９７】また、時系列データｘをｙ＝（ｘ−ｍ）／
σと変換して、平均値の平行移動と分散の変更が可能なモデルを用いることもできる。 この場合、基底の個数と分散のトレードオフおよび分散の推定が必要となるが、
その方法は容易である。

【００９８】さらに、本実施形態においては、１次元入カデータに対する変化トレンドの推定を行っているが、
この推定処理は多次元にもそのまま拡張可能であり、モデルの解釈法などの変更はない。

【００９９】次に、（３）式および（６）式の導出方法について説明する。 一般に、隠れ素子数ｍ個、各隠れ素子に対するレジスタ数ｓ個を有するＣＳＳＲＮＮを考える。 また、入力データの次元はｄ次元とする。 このとき、入力重み変数もまたｄ次元となる。 （３）、（６）
式ではｄ＝１と置いた場合を考えている。

【０１００】ＣＳＳＲＮＮの各隠れ素子に確率変数ｈ _j
∈｛０，１｝（ｊ＝１，．．．，ｍ）を割り付ける。 これらを隠れ変数と呼び、ベクトルｈ＝（ｈ ₁ ，．．．，
ｈ _m ）で表現する。 時刻ｔのＣＳＳＲＮＮのｄ次元の入力 外７ （以後、ベクトル

【０１０１】

【外７】

【０１０２】ｘ ^(t)と記す）と隠れ変数ベクトルｈ ^(t)
とを合わせて、時系列解析系の状態（ベクトルｘ ^(t) ，
ベクトルｈ ^(t) ）と呼ぶことにする。 ＣＳＳＲＮＮの内部状態Ｏ ^(t,s)バーおよびパラメータセットФから状態（ベクトルｘ ^(t) ，ベクトルｈ ^(t) ）のエネルギー関数を次式で定義する。

【０１０３】

【数１７】

【０１０４】このエネルギー関数から、状態（ベクトルｘ ^(t) ，ベクトルｈ ^(t) ）の条件付き確率を次式で定義する。

【０１０５】

【数１８】

【０１０６】ここで、

【０１０７】

【数１９】

【０１０８】である。 （１１）式の条件付き確率から、
ベクトルｘ ^(t) 、ベクトルｈ ^(t)それぞれの条件付き確率を導くことができる。 まず、全ての隠れ変数ベクトルｈ ^(t)上の総和を計算することで、（１１）式からベクトルｈ ^(t)を消去すると、

【０１０９】

【数２０】

【０１１０】を得る。 （１２）式にｄ＝１を代入すると（３）式が得られる。 また、ベクトルｘ ^(t)に関して（１１）式を積分することで、ベクトルｘ ^(t)を消去すると、

【０１１１】

【数２１】

【０１１２】を得る。 ここで、ｈ _j ^(t)はベクトルｈ
^(t)のｊ番目の成分である。 次に、ベイスの法則により（１１）式と（１２）式から条件付き確率

【０１１３】

【数２２】

【０１１４】が求められる。 ここで、

【０１１５】

【数２３】

【０１１６】である。 （１４）式の右辺の帰結の式において、各隠れ変数に関する確率が積の形で含まれることより、明らかに各隠れ変数の独立性を仮定することができる。 そこで、（１５）式を用いて、

【０１１７】

【数２４】

【０１１８】を導くことができる。 （１６）式の右辺は入力ベクトルｘ ^(t)に対するｊ番目の隠れ素子の出力を表しており、同時にまた、隠れ変数ｈ _j ^(t)が１である条件付き確率を与えている。 実際、入力ベクトルｘ ^(t)
の次元が１次元である場合には、（１６）式の右辺は（２）式の右辺に一致する。 このとき、（２）式のＯ _j
^{(t )}は、ｈ _j ^(t) ＝１に対応する時系列データｘ ^(t)の確率密度を与える。 したがって、（１６）式から、ＣＳ
ＳＲＮＮの内部状態の解釈と時系列データとの関係が明らかになる。 また、（１４）式において各隠れ変数の独立性を仮定することは、ＣＳＳＲＮＮのコラム間で情報の交換をしないことと対応している。

【０１１９】次に、学習により最適なパラメータセットが見つかったと仮定する。 このとき、予測装置１７は、
予測確率密度関数ｐ（ベクトルｘ ^(t) ｜Φ，Ｏ ^(t,s)バー）を用いて時刻ｔの予測を行う。 具体的には、最も確からしいベクトルｘ ^(t)の予測値として、確率密度関数のピークに対応する値を選ぶ。 そこで、微分方程式

【０１２０】

【数２５】

【０１２１】を満たすベクトルｘ ^(t)を予測値とする。
（１７）式は、

【０１２２】

【数２６】

【０１２３】と等価である。 ここで、（１２）式のｐ
（ベクトルｘ ^(t) ｜Φ，Ｏ ^(t,s)バー）を（１８）式に代入して、予測方程式となる非線形方程式

【０１２４】

【数２７】

【０１２５】を得る。 入力ベクトルｘ ^(t)の次元が１次元である場合には、（１９）式は（６）式に一致する。
次に、図１２から１４までを参照しながら、図３の時系列トレンド推定システムの動作フローを説明する。

【０１２６】図１２は、時系列トレンド推定システムの全体処理のフローチャートである。 図１２においてシステムが立ち上げられると、まず、あらかじめ決められた初期設定に基づいて制御装置１１が起動される（ステップＳ２１）。 表示・対話装置１２およびシステム管理装置１３は入力待ち状態となり（ステップＳ２２）、クライアントからの要求が入力されると（ステップＳ２
３）、システム管理装置１３は、まず観測対象の変更要求があるかどうかを判定する（ステップＳ２４）。

【０１２７】変更要求があれば、観測対象を変更して時系列データのを収集を開始し（ステップＳ２５）、次に学習アルゴリズムの変更要求があるかどうかを判定する（ステップＳ２６）。 変更要求があれば、指定された項目を修正して必要な計算機資源を確保し（ステップＳ２
７）、次にＣＳＳＲＮＮ１９の規模の変更要求があるかどうかを判定する（ステップＳ２８）。 変更要求があれば、ＣＳＳＲＮＮ１９の規模を修正して必要な計算機資源を確保し（ステップＳ２９）、次に予測装置１７の設定変更要求があるかどうかを判定する（ステップＳ３
０）。 ステップＳ２４で変更要求がなければステップＳ
２６の処理に移り、ステップＳ２６で変更要求がなければステップＳ２８の処理に移り、ステップＳ２８で変更要求がなければステップＳ３０の処理に移る。

【０１２８】ステップＳ３０で変更要求があれば、予測装置１７の設定を変更して再起動し（ステップＳ３
３）、次に学習装置１４の設定変更要求があるかどうかを判定する（ステップＳ３４）。 ステップＳ３０で変更要求がなければ、予測装置１７が起動されているかどうかを調べ（ステップＳ３１）、起動されていればステップＳ３４の処理に移る。 予測装置１７が起動されていなければ、それを起動して（ステップＳ３２）、ステップＳ３４の処理に移る。

【０１２９】ステップＳ３４で変更要求があれば、学習装置１４の設定を変更して再起動し（ステップＳ３
７）、次に観測データの提示要求があるかどうかを判定する（ステップＳ３８）。 ステップＳ３４で変更要求がなければ、学習装置１４が起動されているかどうかを調べ（ステップＳ３５）、起動されていればステップＳ３
８の処理に移る。 学習装置１４が起動されていなければ、それを起動して（ステップＳ３６）、ステップＳ３
８の処理に移る。 ステップＳ３８で提示要求があれば、システム管理装置１３は、観測データを観測装置１
５から表示・対話装置１２に転送し（ステップＳ３
９）、次に予測データの提示要求があるかどうかを判定する（ステップＳ４０）。 提示要求があれば、予測装置１７を呼び出して予測データを受け取り、表示・対話装置１２に転送して（ステップＳ４１）、次に学習の開始要求があるかどうかを判定する（ステップＳ４２）。 要求があれば、学習装置１４を呼び出して学習を行わせ、
予測装置１７のパラメータセットを更新して（ステップＳ４３）、入力待ち状態となる（ステップＳ２２）。 ステップＳ３８で提示要求がなければステップＳ４０の処理に移り、ステップＳ４０で提示要求がなければステップＳ４２の処理に移り、ステップＳ４２で要求がなければステップＳ２２で入力待ち状態となる。

【０１３０】学習装置１４は、制御装置１１からの呼び出しに応じて処理を開始する。 図１３は、学習装置１４
の処理のフロ一チャートである。 図１３において学習装置１４が立ち上げられると、まず制御装置１１からの呼び出しがあるまで待機状態となる（ステップＳ５１）。
呼び出しがあると、学習装置１４は、学習アルゴリズムのパラメータの変更要求があるかどうかを判定する（ステップＳ５２）。 変更要求があれば、新しい設定パラメータを受け取って学習アルゴリズムを更新し（ステップＳ５３）、次に学習基準の変更要求があるかどうかを判定する（ステップＳ５４）。 変更要求があれば、新しい学習基準を受け取ってこれまでの学習基準を更新し（ステップＳ５５）、次に時系列データの学習要求があるかどうかを判定する（ステップＳ５６）。 ステップＳ５２
で変更要求がなければステップＳ５４の処理に移り、ステップＳ５４で変更要求がなければステップＳ５６の処理に移り、ステップＳ５６で学習要求がなければステップＳ５１で待機状態となる。

【０１３１】ステップＳ５６で学習要求があれば、ｉ＝
１とおいて（ステップＳ５７）、パラメータセットΦの次元である dim（Φ）次元の空間のランダムな１点Φ
⁽ⁱ⁾ _{in it}を生成する（ステップＳ５８）。 次に、Φ ⁽ⁱ⁾
_initから初期値シンプレックスＳＬ ⁽ⁱ⁾ ₍₀₎を生成する（ステップＳ５９）。 次に、予測装置１７を呼び出してＳＬ ⁽ⁱ⁾ ₍₀₎の各頂点のΦの値と時系列データ｛ｘ ₁ ，． ． ． ，ｘ _N ｝を与え、対応する予測値列｛ｘ
₁ハット，． ． ． ，ｘ _Nハット｝を受け取る（ステップＳ６０）。 そして、ＳＬ ⁽ⁱ⁾ ₍₀₎の各頂点に対する評価関数 merit（Ｓ｜Φ）の値を計算する（ステップＳ６
１）。

【０１３２】次に、ｊ＝１とおいて（ステップＳ６
２）、制御温度Ｔ ^(j)において評価関数の極小値を与えるパラメータセットを、シンプレックスＳＬ ⁽ⁱ⁾ _(j-1)
から出発して、滑降シンプレックス法とＮ _r回の酔歩の試行により探索する（ステップＳ６３）。 このとき、シンプレックスの収縮過程で新しく生成する頂点の評価関数の値を、予測装置１７と通信しながら計算する（ステップＳ６４）。 そして、新しく得られたシンプレックスをＳＬ ⁽ⁱ⁾ _(j)とする。

【０１３３】次に、Ｔ ^(j)とアニーリングの最終到達温度Ｔ _minとを比較する（ステップＳ６５）。 Ｔ ^(j)がＴ
_minより高ければ、ｊ＝ｊ＋１とおき（ステップＳ６
６）、Ｔ ^(j) ＝κＴ ^(j-1)とおいて温度を下げる（ステップＳ６７）。 ここで、０＜κ＜１である。 そして、ステップＳ６３以降の処理を繰り返す。 ステップＳ６５でＴ ^(j)がＴ _min以下になれば、そのときのシンプレックスＳＬ ⁽ⁱ⁾ _(j)の各頂点の内で評価関数が最小となる点を求め、それを局所最適パラメータセットΦ ⁽ⁱ⁾ _resとして保存する（ステップＳ６８）。

【０１３４】次に、ｉをアニーリングの試行回数Ｎ _aと比較し（ステップＳ６９）、ｉがＮ _aに達していなければ、ｉ＝ｉ＋１とおいて（ステップＳ７０）、ステップＳ５８以降の処理を繰り返す。 ステップＳ６９でｉがＮ
_aに達すると、得られたＮ _a個のΦ ⁽ⁱ⁾ _res （ｉ＝
１，． ． ． ，Ｎ _a ）の中から、評価関数が最小となるものを最適パラメータセットとして選択する（ステップＳ
７１）。 そして、その最適パラメータセットをＣＳＳＲ
ＮＮ１９にセットして（ステップＳ７２）、待機状態に戻る（ステップＳ５１）。

【０１３５】予測装置１７は、学習モードにおいては学習装置１４から呼び出され、予測モードにおいては制御装置１１から呼び出される。 両モードに置ける予測装置１７の基本的な動作は同じであり、観測データの供給元と予測結果の出力先が異なるだけである。 図１４は、予測装置１７の処理のフロ一チャートである。

【０１３６】図１４において予測装置１７が立ち上げられると、まず制御装置１１または学習装置１４からの呼び出しがあるまで待機状態となる（ステップＳ８１）。
呼び出しがあると、予測装置１７は、ＣＳＳＲＮＮ１９
のパラメータセットの変更要求があるかどうかを判定する（ステップＳ８２）。 変更要求があれば、新しいパラメータセットを受け取り、これまでのパラメータセットを更新して（ステップＳ８３）、次に予測の要求があるかどうかを判定する（ステップＳ８４）。 ステップＳ８
２で変更要求がなければステップＳ８４の処理に移り、
ステップＳ８４で予測要求がなければステップＳ８１で待機状態となる。

【０１３７】ステップＳ８４で予測要求があれば、ＣＳ
ＳＲＮＮ１９の内部状態Ｏ ^(t,s)バーをセットし（ステップＳ８５）、ｔ＝１とおいて予測を開始する（ステップＳ８６）。 まず、非線形方程式（６）を解いて時系列データの予測値ｘ ^(t)ハットを求め（ステップＳ８
７）、予測結果として出力する（ステップＳ８８）。 次に、時刻ｔにおける真値ｘ ^(t)を読み込み（ステップＳ
８９）、ＣＳＳＲＮＮ１９を駆動して内部状態Ｏ ^(t,s)
バーをＯ ^(t+1,s)バーに更新する（ステップＳ９０）。
次に、ｔが時刻の上限値Ｎを越えたかどうかを判定し（ステップＳ９１）、Ｎを越えていなければ、ｔ＝ｔ＋
１とおいて（ステップＳ９２）、ステップＳ８７以降の処理を繰り返す。 そして、ステップＳ９１でｔがＮを越えれば、予測を終了し（ステップＳ９３）、待機状態に戻る（ステップＳ８１）。

【０１３８】次に、図３の時系列トレンド推定システムによるデータ変化のトレンドの推定例について説明する。 以下では、記述の簡便さを図って、ｍ個の隠れ素子を有し、各隠れ素子にそれぞれｓ個のレジスタが割り付けられているＣＳＳＲＮＮをｈｍｒｓ−ＮＮと表記する。 隠れ素子に割り付けられたレジスタの数ｓは、文脈層の深さを表す。 例えば、隠れ素子２個、各隠れ素子毎のレジスタ数が１個のＣＳＳＲＮＮは、ｈ２ｒ１−ＮＮ
と表される。

【０１３９】まず、区分的に定常となるデータにおけるトレンド推定結果について説明する。 区分的に定常な時系列生成装置（不図示）から生成された次のような時系列を、予測の対象として選ぶ。

【０１４０】

【数２８】

【０１４１】ただし、Ｎ（μ，σ ² ）は、平均μ、分散σ ²のガウス密度関数を表す。 （２０）式により生成される時系列の例は図１５に示されている。 図１５の時系列において、時刻ｔの４つの区間１≦ｔ≦５０、５１≦
ｔ≦１００、１０１≦ｔ≦１５０、１５１≦ｔ≦２００
の境界における平均値の不連続な切り替わりがデータの跳躍トレンドとして現れていることが伺える。 このような不連続な跳躍トレンドを時系列トレンド推定システムにより予測する。 このとき、最適パラメータセットの決定に用いられる学習データと、予測結果と比較するテストデータは、（２０）式に従ってそれぞれ個別に生成されるものとする。

【０１４２】ｈ２ｒ１−ＮＮまたはｈ２ｒ１０−ＮＮを備えるシステムに、図１５の時系列を学習データとして与えて学習させ、未知のテストデータに対するそれぞれのシステムによる予測結果を求める。 ｈ２ｒ１−ＮＮ、
ｈ２ｒ１０−ＮＮによる予測結果は、それぞれ図１６、
１７のようになる。 図１６、１７において、実線はＣＳ
ＳＲＮＮが予測した結果を表し、破線は時系列生成装置が実際に生成したテストデータを表す。

【０１４３】これらの結果を見ると、確かに跳躍型のトレンドがうまく予測されていることが分かる。 しかし、
支脈層の規模が小さいｈ２ｒ１−ＮＮの場合は、交替現象とノイズの振幅の揺らぎとがうまく分割できていない。 図１６において、３５１≦ｔ≦４００の区間にある予測値Ｄ１は、周期的なトレンドから逸脱している。 これに対して、図１７のｈ２ｒ１０−ＮＮの場合はこうした予測値の逸脱がなく、より正確にトレンドを予測しているといえる。 このように、交替現象の追従にある程度のロバスト性（ノイズ等に引き摺られない強靱さ）を持たせるためには、文脈層の規模をある程度の大きさに設定する必要がある。

【０１４４】次に、異常値を含む時系列データに対するトレンドの推定結果について説明する。 ノイズの摺らぎでは説明できない極端な偏差として、異常値の問題がある。 図１８は、異常値が混入している跳躍型トレンドの時系列を示している。 図１８の時系列は、基本的なデータを（２０）式を用いて生成し、その上に異常値Ｅ１、
Ｅ２、Ｅ３、Ｅ４を加算することにより得られる。 これらの異常値は、発生間隔がポアソン分布に従うように生成される。 またそれらの値は平均１０、分散０．５の正規分布に従って生成され、確率１／２で正値として、また確率１／２で負値として跳躍型トレンドの時系列に加重されるものとする。 こうした異常値を生成する確率に関する情報は、時系列推定システムには全く与えられない。 ここで、異常値とは、例えば正規分布を仮定するなら、その平均値からの逸脱が大きく出現確率がほぼ０であるといえるような値のことである。

【０１４５】ＣＳＳＲＮＮとしてｈ２ｒ１−ＮＮ、ｈ２
ｒ１０−ＮＮ、およびｈ１０ｒ１０−ＮＮを選び、それらを備えるシステムにそれぞれ図１８に示した時系列を学習させる。 そして、それらのシステムによる予測結果を、図１８と同様の機構で生成された別のデータでテストした結果をそれぞれ図１９、２０、２１に示す。 図１
９、２０、２１において、データＥ５、Ｅ６、Ｅ７、Ｅ
８、Ｅ９、Ｅ１０、Ｅ１１、Ｅ１２、Ｅ１３、Ｅ１４、
Ｅ１５、Ｅ１６、Ｅ１７、Ｅ１８は、テストデータの異常値を表す。

【０１４６】図１９のｈ２ｒ１−ＮＮの場合は、予測値における異常値はＤ２、Ｄ３、Ｄ４、Ｄ５、Ｄ６、Ｄ
７、Ｄ８、Ｄ９、Ｄ１０、Ｄ１１の１０個である。 これらの異常な予測値は、それぞれ直前のテストデータの異常値Ｅ６、Ｅ７、Ｅ１０、Ｅ１２、Ｅ１３、Ｅ１４、Ｅ
１５、Ｅ１６、Ｅ１７、Ｅ１８の影響を受けて生成されたものと考えられる。 これに対して、図２０のｈ２ｒ１
０−ＮＮの場合は、予測値における異常値はＤ１２、Ｄ
１３、Ｄ１４、Ｄ１５、Ｄ１６の５個だけである。 さらに、図２１のｈ１０ｒ１０−ＮＮの場合は、予測値における異常値はＤ１７のみとなっている。 図２１において、予測値のトレンドが小刻みに上下しているのは、多数の隠れ素子により離散値がより細かく表現されることに起因する。 これらの結果より、ＣＳＳＲＮＮにおいて回路の規模を大きくすると、より大きなロバスト性を実現できることがわかる。 したがって、本発明のシステムを用いれば、時系列の異常値が存在していても不連続トレンドを的確に推定することができる。

【０１４７】次に、時系列ＳとパラメータセットΦにより決まるモデルとの適合度を測る関数 merit（Ｓ｜Φ）
（最小２乗法や最尤推定法等における評価関数）を用いて、ＣＳＳＲＮＮを予測フィルタとして構成した時に、
時系列の性質がどのように内部表現として獲得されるのかを考える。

【０１４８】まず、隠れ変数ベクトルとガウス密度関数の関係について考察する。 ＣＳＳＲＮＮの隠れ変数ベクトルｈ ^(t)と入力ベクトルｘ ^(t)上において、同時密度関数ｐ（ベクトルｘ ^(t) ，ベクトルｈ ^(t) ｜Φ，Ｏ
^(t,s)バー）が（１１）式により定義される。 この同時密度関数から導かれるベクトルｘ ^(t)の周辺密度ｐ（ベクトルｘ ^(t) ｜Φ，Ｏ ^(t,s)バー）とベクトルｈ ^(t)の周辺分布関数Ｐ（ベクトルｈ ^(t) ｜Φ，Ｏ ^(t,s)バー）
の間には、ガウス密度関数を基底とする線形関係がある。 このことは、以下の考察から容易に分かる。

【０１４９】時刻ｔにおけるＣＳＳＲＮＮの隠れ変数ベクトルの予測分布は、（１３）式のＰ（ベクトルｈ ^(t)
｜Φ，Ｏ ^(t,s)バー）によって決まる。 ＣＳＳＲＮＮが表現する事象は、隠れ変数ベクトルが確率的に重なった状態である。 以後、表現上の簡約のため、予測における隠れ変数ベクトルの確率的重ね合わせのことを様相と呼ぶことにする。 （１３）式は、可能な全ての隠れ変数ベクトルの集合上で定義され、それらの重率を予測する式である。

【０１５０】ここで、隠れ変数ベクトルを要素表示するとともに、各要素に非負整数を指標として割り当てる。
隠れ素子の個数をｍとすると、２ ^m個の隠れ変数ベクトルは、

【０１５１】

【数２９】

【０１５２】のように表記される。 例えば、ｍ＝２のときは、隠れ変数ベクトルはベクトルｈ ₍₀₎ ＝（０，
０）、ベクトルｈ ₍₁₎ ＝（１，０）、ベクトルｈ ₍₂₎ ＝
（０，１）、ベクトルｈ ₍₃₎ ＝（１，１）の４つである。

【０１５３】時刻ｔにおける入力ベクトルの予測分布は、（１２）式より、

【０１５４】

【数３０】

【０１５５】と書き換えられる。 ここで、 外８ （以後、ベクトルｗ _jと記す）はｊ番目の

【０１５６】

【外８】

【０１５７】隠れ素子の入力重みベクトルであり、

【０１５８】

【数３１】

【０１５９】である。 これは、平均値

【０１６０】

【数３２】

【０１６１】、分散１のガウス密度関数である（以後、
形式的に

【０１６２】

【数３３】

【０１６３】と表記する）。 （２２）式は、（１３）式を用いて次のように書き換えることができる。

【０１６４】

【数３４】

【０１６５】ここで、

【０１６６】

【数３５】

【０１６７】である。 （２３）式は、入力ベクトルの予測密度関数と隠れ変数ベクトルにより指定されるガウス密度関数の線形結合で表現されている。 その結合係数は、対応する隠れ変数ベクトルの予測分布である。 したがって、これらの結合係数はＣＳＳＲＮＮの内部状態の非線形関数になる。 各隠れ変数ベクトルｈ _(i) ＝
（ｈ _i1 ，．．．，ｈ _im ）を入力ベクトル空間上のガウス密度関数に対応させる写像

【０１６８】

【数３６】

【０１６９】は、各隠れ素子の重みの集合｛ベクトルｗ
₁ ，． ． ． ，ベクトルｗ _m ｝によって決まる。 ＣＳＳ
ＲＮＮとして、隠れ素子２個で入力ベクトルの次元が１
の場合を考える。 隠れ素子の入力重みをそれぞれｗ ₁ 、
ｗ ₂とすると、各隠れ変数ベクトルに対して基底となるガウス密度関数は、それぞれ、

【０１７０】

【数３７】

【０１７１】となる。 今、仮に隠れ変数ベクトルの時刻ｔにおける予測分布が次のように求められたとする。

【０１７２】

【数３８】

【０１７３】このとき、（２３）式により入力ｘ ^(t)の予測分布は、

【０１７４】

【数３９】

【０１７５】となる。 この場合、およそＮ（ｗ ₁ ，１）
に近い形状の分布となることが分かる。 一般に、時刻ｔ
における入力ｘ ^(t)の予測密度の形状は、入力重み間の距離や基底となるガウス密度関数に対する重みに応じて変わる。 予測密度関数のすそが広がることもあれば、いずれの基底関数のピーク（平均値）に対しても予測密度関数のピークとの間に偏差があることもある。 また、予測密度関数のピークが多数個あることもあり、その他にも様々に変化し得る。 隠れ変数ベクトルの予測分布が、
特定の隠れ変数ベクトルｈ _(i)だけに極端に偏る場合は、入力に対する予測密度関数はＮ（Σｈ _ij・ベクトルｗ _j ，１）にほぼ一致する。 このとき、その平均値Σｈ
_ij・ベクトルｗ _jは、ｈ _ij ＝１となる隠れ素子の入力重みベクトルｗ _jにより決められることになる。

【０１７６】上述のような関係を急激なトレンド変化を伴う時系列の予測実験に利用してみる。 まず、次式に示す平均値移動型の時変ガウス分布により時系列Ｓを生成する。

【０１７７】

【数４０】

【０１７８】そして、隠れ素子が２個、各隠れ素子に割り付けられたレジスタが１個のｈ２ｒ１−ＮＮを時系列Ｓに適合させる。 そのために、負値対数尤度により定義された（９）式の merit（Ｓ｜Φ）を最小にするパラメータΦ _opt ＝｛ｗ ₁ ，θ ₁ ，ｕ ₁ ，ｗ ₂ ，θ ₂ ，ｕ ₂ ｝
を、数値的最適化法により見つける。 適当なパラメータが選択されると、それをΦ _optとしてｈ２ｒ１−ＮＮを構成し、図２２に示すテスト用の時系列の予測を行う。
このときの予測の対象は、刻々と入力されるテストデータの次の時刻の入力値である。 ｈ２ｒ１−ＮＮの内部状態の更新は、実際に入力される観測値を用いて行う。 このような設定は、オープンループと呼ばれる。 図２３
は、ｈ２ｒ１−ＮＮによる予測結果を示している。 図２
３の予測値は、入力ｘ ^(t)の予測分布のピークに対応している。 図２２と図２３を比較すると、ｈ２ｒ１−ＮＮ
はほぼ正確に時系列のトレンドの不連続な切り替えを追跡し、各トレンドをうまく予測していることが分かる。
予測密度全体の時間的変化（予測密度関数列）は図２４
に示されている。 図２４において、入力値ｘの予測確率密度Ｐが時間ｔの経過とともに変化する様子が示されている。 各時刻の予測確率密度Ｐのピークに対応するｘの値を予測値としてプロットしたものが図２３である。

【０１７９】ここで、ＣＳＳＲＮＮが時系列Ｓのどのような性質をその内部表現として獲得しているかを考察する。 まず、ｈ２ｒ１−ＮＮの各隠れ変数ベクトルの予測分布と入力の予測分布との関係を定性的に説明するために必要な情報として、予測に用いたｈ２ｒ１−ＮＮの最適パラメータセットを図２５に示す。 図２５の各パラメータの値から、隠れ変数ベクトルとガウス密度関数との対応関係が、

【０１８０】

【数４１】

【０１８１】となることが分かる。 次に、ｈ２ｒ１−Ｎ
Ｎが予測する隠れ変数ベクトルｈ ₍₀₎ ＝（０，０）、ベクトルｈ ₍₁₎ ＝（１，０）、ベクトルｈ ₍₂₎ ＝（０，
１）、ベクトルｈ ₍₃₎ ＝（１，１）に対する確率密度を、それぞれ図２６、２７、２８、２９に示す。 図２３
の時系列のトレンド予測と図２４の予測密度全体の変化、および図２６から図２９までの各隠れ変数ベクトルの予測分布を比較検討すると、次のようなことが言える。

【０１８２】Ｎ（１，１）から生成される時系列の区間では、隠れ変数ベクトルｈ ₍₁₎の予測分布が係数の中で主項となり、予測密度の形状は主にＮ（１．０２６９
６，１．０）によって形作られる。 また、時系列Ｓに対してｈ２ｒ１−ＮＮが隠れ変数ベクトルの予測として適切な様相を実現するために、無視できない割合でベクトルｈ ₍₀₎の分布Ｎ（０．０，１．０）が重なり合うようなパラメータが選ばれている。 このため、予測密度関数の負方向のすそが広がり、そのピークがＮ（１．０２６
９６，１．０）のピークからｘの負方向にずれていることが分かる。 ここで、適切な様相を実現することは、入力に対する予測密度から計算される負値対数尤度を最小にすることを意味する。

【０１８３】Ｎ（−１，１）から生成される時系列の区間では、隠れ変数ベクトルｈ ₍₂₎の予測分布が係数の中で主項となり、予測密度の形状は主にＮ（−１．１０８
５３，１．０）によって形作られる。 Ｎ（−０．０８１
５７，１．０）の重みが無視できるほど小さくないので、予測密度関数の正方向のすそが広がり、そのピークがＮ（−１．１０８５３，１．０）のピークからｘの正方向にずれていることが分かる。

【０１８４】以上の議論は、（２３）式により理論的に裏付けられている。 この実験結果から次のような解釈が得られる。 時系列Ｓに対してＣＳＳＲＮＮの最適なパラメータとは、ＣＳＳＲＮＮによって表現可能な予測密度関数列の中から尤もらしい関数列を生成するものである。 最適なパラメータを設定することにより、ＣＳＳＲ
ＮＮは時系列Ｓの背後にある法則性（例えば、トレンドの切り替えなど）を相異なる様相から様相への変化として表現する。 時系列Ｓ自体は、様相の列によって指定された予測密度関数列により記述される。

【０１８５】実験結果においては、ＣＳＳＲＮＮの予測密度関数のモードによってトレンドが表現されている。
ここで、モードとは、ある時刻における予測密度関数の最大ピークに対応するデータ値（最頻値）を指す。 図２
６、２７、２８、２９から明らかなように、様相としては大きく２つに大別される。 大別された各様相を１つのガウス密度関数基底Ｎ（Σｈ _ij・ベクトルｗ _j ，１）によって高い精度で近似できるならば、そのガウス密度関数が１つのトレンドを近似的に表現する。 この場合、そのトレンドはＣＳＳＲＮＮの重みベクトルとして明示的に分散表現されていることが分かる。

【０１８６】また、時系列のトレンドが、様相の主項となる隠れ変数ベクトルに対応するガウス密度関数基底Ｎ
（Σｈ _ij・ベクトルｗ _j ，１）のピークと数値的にほぼ一致する場合、そのトレンドのデータ値はＣＳＳＲＮＮ
の入力荷重Σｈ _ij・ベクトルｗ _jにより分散表現される。

【０１８７】次に、予測フィルタとしてＣＳＳＲＮＮを用いる場合、文脈層がどのような役割を果たすかについて考える。 （２５）式により生成される時系列の予測フィルタとして要求されることは、時系列のトレンドを予測することとトレンドが切り替わったときにそれを的確に追跡することである。

【０１８８】（２５）式から導かれる性質として、次のようなものを挙げることができる。 文脈層は時間的に変化し、予測密度関数の形状を変える。 文脈層の内容Ｏ
^(t,s)バーが変化することにより様相が変化し、予測密度関数の形状が変わる。

【０１８９】

【数４２】

【０１９０】は、ｊ番目の隠れ素子の出力履歴が様相に関与する度合を示す。 （１３）式において、Θ _j ^(t)が大きければｊ番目の成分が１である隠れベクトルの重率が大きくなる。

【０１９１】では、具体的な時系列とＣＳＳＲＮＮの文脈層の回帰係数およびしきい値との関係を考える。 ここでは、ｈ２ｒ１−ＮＮとの比較のためにｈ２ｒ２−ＮＮ
による予測値を図３０に示し、その予測密度関数列を図３１に示す。 図２３と図３０とを比較するとｈ２ｒ２−
ＮＮのトレンド予測には乱れが少なく、安定な予測を行っていることが分かる。 また、図３１のどの時刻における予測密度関数の形状もガウス密度関数に近い。 ｈ２ｒ
２−ＮＮで使用した最適なパラメータセットは図３２に示されている。 ｈ２ｒ２−ＮＮの場合も、ｈ２ｒ１−Ｎ
Ｎと同様にして、予測フィルタとして適当なパラメータを見つけることができる。

【０１９２】次に、ｈ２ｒ２−ＮＮが予測する隠れ変数ベクトルｈ ₍₀₎ ＝（０，０）、ベクトルｈ ₍₁₎ ＝（１，
０）、ベクトルｈ ₍₂₎ ＝（０，１）、ベクトルｈ ₍₃₎ ＝
（１，１）に対する確率密度を、それぞれ図３３、３
４、３５、３６に示す。 これらの図から、各区間の様相はそれぞれ１つの隠れ変数ベクトルでほぼ完全に表現できていることが分かる。 Ｎ（１，１）が支配する区間に対してはベクトルｈ ₍₂₎が対応し、また、Ｎ（−１，
１）の区間に対してはベクトルｈ ₍₃₎が対応する。 したがって、トレンドが入力荷重により分散表現されていることが分かる。 以上のことから、ｈ２ｒ１−ＮＮに比べて、ｈ２ｒ２−ＮＮの方が予測フィルタとしての性質をより多く満足していると考えられる。

【０１９３】図３７、３８は、特別なテストデータとそれに対するｈ２ｒ１−ＮＮとｈ２ｒ２−ＮＮとによる予測結果を示している。 各ＣＳＳＲＮＮのパラメータとしては、それぞれ図２５、３２に示したものを用いた。 このテストデータは、｛−１．０，−０．５，−０．２，
−０．１，０．０，０．１，０．２，０．５，１．０｝
のいずれかの値がある程度続いた後、不連続に切り替わるパターンを表している。 そして、それらの値の組み合わせを適宜決めることにより生成されている。 この実験の目的は、文脈層の回帰係数（レジスタ係数）としきい値の役割を明らかにすることである。 すなわち、最適パラメータのＣＳＳＲＮＮにおいて、入力ｘ ^(t)に対する様相の変化およびトレンドの切り替わりを定性的に把握することが狙いである。 図３７、３８において、黒塗りの四角形がテストデータを表し、白抜きの四角形がｈ２
ｒ１−ＮＮによる予測値を表し、白抜きの丸がｈ２ｒ２
−ＮＮによる予測値を表す。

【０１９４】テストデータの全区間（時間軸をｔとして、ｔ∈［１，２，．．．，２３０］）において、全体的な傾向として、ｈ２ｒ１−ＮＮに比べてｈ２ｒ２−Ｎ
Ｎの方がトレンドからの変位（変差）の大きなデータに対して鈍感であり、トレンドの切り替わりに対する追従の速さに遅れがある。

【０１９５】図３７のｔ∈［１，． ． ． ，３０］においては、Ｎ（１，１）またはＮ（−１，１）の平均値をトレンドとして推定している状態で、分散１を越える変位のデータが突然入った場合、以後の予測がどのように変化するかを調べるためにテストデータを構成した。 例えば、Ｎ（１，１）の平均値１をトレンドと推定している状態で、その予測から外れてＮ（−１，１）の平均値−
１の値に近い入力ｘ ^{(t )} ＝−０．５が突然入った場合、
ｈ２ｒ１−ＮＮの予測値もｈ２ｒ２−ＮＮの予測値も僅かに揺らぐだけである。 その予測値の乱れは、若干、ｈ
２ｒ１−ＮＮの方が大きい。 Ｎ（−１，１）の平均値−
１をトレンドとして推定している状態で、Ｎ（１，１）
の平均値１の値に近い入力ｘ ^(t) ＝０．５が突然入った場合についても同様である。 ｈ２ｒ１−ＮＮ、ｈ２ｒ２
−ＮＮの内部状態Ｏ ^(t,1)バー、Ｏ ^(t,2)バーから決まるそれぞれの様相は、トレンドからの変差が分散を越えるような入力が３個程度含まれても現状を保持する。

【０１９６】ｔ∈［３１，． ． ． ，６０］においては、
現在の予測がトレンドと一致しているものとして、トレンドそのものではないが、その周辺の値、つまり分散よりも小さな変位の値が連続して入力された場合、予測がどのように変化したかを調べるためにデータを構成した。 この場合、ｈ２ｒ１−ＮＮ、ｈ２ｒ２−ＮＮともに予測値に変動は見られない。 この区間では、内部状態Ｏ
^(t,1)バー、Ｏ ^(t,2)バーの変化による様相の変化はほとんど生じないと考えられる。

【０１９７】ｔ∈［６１，． ． ． ，１００］においては、現在の予測がトレンドと一致しているものとして、
平均値から分散以内の変差を持った入力が連続している状態で、分散の２倍の変差を持ったデータが突然入力された場合を考える。 これは、揺らぎや異常値を模したデータが入力された後の予測の変動を調べるためである。
予測値にはある程度の乱れが生じ、入力の変化の影響が僅かながら持続する。 ｈ２ｒ１−ＮＮに比べてｈ２ｒ２
−ＮＮの方が、影響は残らないことが読みとれる。 内部状態Ｏ ^(t,1)バー、Ｏ ^(t,2)バーが変化し、それらの変化は様相の変化として現れている。 このとき、予測密度関数のピークの移動は小さいが、ガウス密度関数に比べると予測密度関数のすそは広がる傾向にある。

【０１９８】図３８のｔ∈［１０１，． ． ． ，１５０］
においては、現在の予測がトレンドと一致しているものとして、変差が分散よりも大きく、分散の２倍以下であるようなデータが連続して入力される場合の予測の変化を調べた。 予測の変動は、ｈ２ｒ１−ＮＮとｈ２ｒ２−
ＮＮではかなり異なっている。 現在のトレンドを−１として入力ｘ ^(t) ＝０．５を選び、その値を連続的に入力し続けた場合、予測密度関数は徐々にＮ（１，１）に近い形に変化する。 ｈ２ｒ１−ＮＮでは連続的かつ直線的に変化している。 一方、ｈ２ｒ２−ＮＮでは指数関数的に急激に変化してＮ（−１，１）に近い形からＮ（１，
１）に変化する。 ｈ２ｒ１−ＮＮに比べてｈ２ｒ２−Ｎ
Ｎでは、様相の保持が強固であることが分かる。

【０１９９】ｔ∈［１５１，． ． ． ，２３０］においては、内部状態の変化が様相の変化としてほとんど現れないという意味で保守的であることを示している。 ここでの入力はトレンドに対して分散を越える変差を持つデータではあるが、先の区間で用いられたデータよりも変差の小さいものを用いた。 ｈ２ｒ１−ＮＮでは予測密度関数のピークが連続的に移動するが、ｈ２ｒ２−ＮＮでは全く動かないといえる。 ｈ２ｒ２−ＮＮにおいて、入力の変差が分散を越える大きさであっても、更新された内部状態の変化が様相の変化としてほとんど現れていない。

【０２００】以上、簡単にまとめると、（２５）式で生成されたデータに対して適切なパラメータを選択するとき、平均値が一定に保たれている各区間において、ＣＳ
ＳＲＮＮのその区間における様相が、入力の変動に関わらず、できるかぎり保存されるようなパラメータが選ばれている。 このような様相の保存という要求がある一方で、様相の切り替えを的確に行うという要求にも答えなければならない。 したがって、文脈層のレジスタに割り付ける重みベクトルｕ _jとしきい値θ _jとして、ＣＳＳ
ＲＮＮの様相が保守的になり、かつ、様相から様相への遷移では中間的な様相を経由しないような値を選択することができればよい。 それらの値は、適当な個数のレジスタを持つＣＳＳＲＮＮにおいて、負値対数尤度を最小にするパラメータを探索することによって得られる。

【０２０１】ＣＳＳＲＮＮの予測分布は、入力および隠れ素子の過去の出力に応じて、２つのトレンドに対応する主な２つの様相の間を移り合う。 様相間の切り替えは揺らぎに対しては保守的であると同時に、実際の切り替えに対しては不連続的に対処するように、しきい値およびレジスタの重み係数が選ばれている。

【０２０２】次に、レジスタの個数について考えてみる。 様相の切り替えを的確に行うという要求に加えて、
計算効率という点も考えて文脈層の大きさを考える必要がある。 多数のレジスタを持つＣＳＳＲＮＮを予測フィルタとして用いる場合、適切なパラメータを選択する時の困難さやトレンドの切り替え時の動作遅れなど好ましくない側面が浮かび上がる。 レジスタの個数を増やせば負値対数尤度を小さくすることはできるが、予測値が時系列に引きずられるため予測フィルタとしては適切ではない。 したがって、時系列に対してレジスタの個数を適切に算定する必要がある。 そこで、単純なヒューリスティクスとして、形式的に情報量基準を援用し、レジスタの個数（文脈層の深さ）の適切さを計算する方法を導入する。

【０２０３】ここでは、（２５）式によって生成された時系列Ｓに対して、いくつかのｈ２ｒｓ−ＮＮ（１≦ｓ
≦５）の形式的情報量基準 ２×（負値対数尤度）＋２×（隠れ素子１つ当たりのレジスタの個数） を計算する。 その結果は図３９に示されている。

【０２０４】図３９から、分布のすそが重なる領域に出現するデータや異常値に対して各区間に対応する様相を保存するには、予測レジスタの個数を２個として過去２
単位時間の出力を保存すれば十分であることが分かる。
レジスタが２個以上のＣＳＳＲＮＮでは、レジスタ数の増加に伴う負値対数尤度の変化は小さい。 また、トレンドの予測値の形状に大きな差が見られないことや、レジスタの個数が増えると逆にＣＳＳＲＮＮのトレンドの予測がテストデータに引きずられる傾向があることなどを実験的に確かめることができる。 このような形式的情報量基準を用いる評価方法は、（２５）式のような構造の情報源に対する予測フィルタを構成するという目的に合致したヒューリスティクスと考えられる。

【０２０５】以上の実験および考察より、ＣＳＳＲＮＮ
による予測フィルタの内部表現と時系列データの関係について次のような結論が得られた。 （ｅ）ＣＳＳＲＮＮにおいて、予測密度関数ｐ（ベクトルｘ ^(t) ｜Φ，Ｏ ^(t,s)バー）が隠れ変数ベクトルｈ
_(i)の予測分布Ｐ（ベクトルｈ _(i) ｜Φ，Ｏ ^(t,s)バー）を重みとするガウス密度関数Ｎ（Σｈ _ij・ベクトルｗ _j ，１）の線形結合で表現される。 基底となるガウス密度関数の平均値は、隠れ変数ベクトルと隠れ素子の入力重みベクトルｗ _jとにより決まり、Σｈ _ij・ベクトルｗ _jとなる。 したがって、予測密度関数の形状は、隠れ変数ベクトルが確率的に重合した状態（様相）と入力重みベクトルとにより決まる。 （ｆ）時系列Ｓに対してＣＳＳＲＮＮの最適なパラメータとは、ＣＳＳＲＮＮによって表現可能な予測密度関数列の中から尤もらしい関数列を生成するものである。 ゆえに、最適なパラメータのＣＳＳＲＮＮでは、時系列Ｓ
の持っている性質（確率的な構造）が様相の重率に反映される。 （ｇ）時系列Ｓが平均値移動型のガウス密度関数によって生成される場合、ＣＳＳＲＮＮは時系列Ｓの背後にあるトレンドの切り替えなどの法則性を相異なる様相から様相への遷移として表現する。 時系列Ｓそれ自体は、様相の列によって指定された予測密度関数列に従うものと見倣される。 （ｈ）時系列のトレンドが、様相の主項となる隠れ変数ベクトルに対応するガウス密度関数基底Ｎ（Σｈ _ij・ベクトルｗ _j ，１）のピークと数値的にほぼ一致し、かつ、主項の重率が１．０に近ければ、トレンドの値はＣ
ＳＳＲＮＮの入力荷重Σｈ _ij・ベクトルｗ _jにより分散表現される。 （ｉ）作動中のＣＳＳＲＮＮでは、予測密度関数は、入力および隠れ素子の過去の出力に応じて、複数のトレンドに対応する主な複数の様相間を移り合う。 平均値移動型ガウス密度関数に対する予測フィルタとしては、様相間の切り替えが揺らぎに対して保守的であると同時に、
それが不連続的に実施されることが要求される。 したがって、内部状態の変動が様相の変動に可能な限り影響しないという条件で、可能な限り内部状態の変動の許容範囲を広げるように、しきい値およびレジスタの重み係数が選ばれなければならない。 このとき、評価関数 merit
（Ｓ｜Φ）を数値的に最適化することで、要求を満たすパラメータが見つかる。 （ｊ）予測フィルタを実現するために、ＣＳＳＲＮＮの大きさ（隠れ素子の個数とレジスタの個数）を適切に選択する必要がある。 平均値移動型ガウス密度関数に対して、形式的情報量基準を用いて文脈層の深さ（レジスタの個数）および各パラメータの値を定めると、適切な予測フィルタを構成することができる。

【０２０６】上述した実施形態から分かるように、本発明の時系列トレンド推定システムは、不連続に急変するトレンドを的確にトレースすることができる。 したがって、次のような非定常非線形の時系列解析を必要とする分野において、高速で的確なトレンドの推定に利用できる。 ♯１：コンピュータ・ネットワークの資源管理 ネットワーク・トラフィックやネットワーク全体のＣＰ
Ｕ負荷を、時系列データとして推定することができる。
また、利用可能なメモリ量を用いてネットワークの利用度をモデル化し、利用度の変化を推定することができる。

【０２０７】ネットワーク・トラフィックなどの時間的に変化する測定量からネットワークの状態を予測することにより、ネットワークの資源を適切に制御することができる。 特に、測定量が区分的に定常であって、各区間毎に不連続に変化するトレンドを含む場合でも、ネットワークの状態が的確に予測される。 ♯２：生体信号処理 脳波による睡眠ステージの個人別のモデル化を行い、睡眠障害の予測を行うことができる。 また、心電による個人別の負荷モデルを作成し、身体の異常の予測を行うことができる。

【０２０８】

【発明の効果】本発明によれば、再帰型ニューラルネットワークの内部状態と時系列データとの関係が明確になり、その関係を用いて時間的に不連続に変化する時系列のトレンドを効率的に推定することが可能になる。

【図面の簡単な説明】

【図１】本発明の原理図である。

【図２】時系列トレンド推定処理のフローチャートである。

【図３】実施形態のシステム構成図である。

【図４】実施形態の計算機システムの構成図である。

【図５】ネットワーク学習装置の構成図である。

【図６】２次元のパラメータセットの初期値シンプレックスを示す図である。

【図７】予測装置の構成図である。

【図８】ニューラル素子を示す図である。

【図９】目標関数評価装置の構成図である。

【図１０】正規化係数算定装置の構成図である。

【図１１】確率密度算定装置の構成図である。

【図１２】トレンド推定システムの処理フローチャートである。

【図１３】学習装置の処理フローチャートである。

【図１４】予測装置の処理フローチャートである。

【図１５】跳躍型トレンドの時系列を示す図である。

【図１６】ｈ２ｒ１のニューラルネットワークによる予測結果を示す図である。

【図１７】ｈ２ｒ１０のニューラルネットワークによる予測結果を示す図である。

【図１８】異常値を含む時系列を示す図である。

【図１９】異常値を含む場合のｈ２ｒ１のニューラルネットワークによる予測結果を示す図である。

【図２０】異常値を含む場合のｈ２ｒ１０のニューラルネットワークによる予測結果を示す図である。

【図２１】異常値を含む場合のｈ１０ｒ１０のニューラルネットワークによる予測結果を示す図である。

【図２２】テスト用の時系列を示す図である。

【図２３】ｈ２ｒ１のニューラルネットワークによる予測値を示す図である。

【図２４】ｈ２ｒ１のニューラルネットワークによる予測分布を示す図である。

【図２５】ｈ２ｒ１のニューラルネットワークの最適パラメータセットの例を示す図である。

【図２６】ｈ２ｒ１のニューラルネットワークにおける隠れ変数ベクトル（０，０）の確率密度を示す図である。

【図２７】ｈ２ｒ１のニューラルネットワークにおける隠れ変数ベクトル（１，０）の確率密度を示す図である。

【図２８】ｈ２ｒ１のニューラルネットワークにおける隠れ変数ベクトル（０，１）の確率密度を示す図である。

【図２９】ｈ２ｒ１のニューラルネットワークにおける隠れ変数ベクトル（１，１）の確率密度を示す図である。

【図３０】ｈ２ｒ２のニューラルネットワークによる予測値を示す図である。

【図３１】ｈ２ｒ２のニューラルネットワークによる予測分布を示す図である。

【図３２】ｈ２ｒ１のニューラルネットワークの最適パラメータセットの例を示す図である。

【図３３】ｈ２ｒ２のニューラルネットワークにおける隠れ変数ベクトル（０，０）の確率密度を示す図である。

【図３４】ｈ２ｒ２のニューラルネットワークにおける隠れ変数ベクトル（１，０）の確率密度を示す図である。

【図３５】ｈ２ｒ２のニューラルネットワークにおける隠れ変数ベクトル（０，１）の確率密度を示す図である。

【図３６】ｈ２ｒ２のニューラルネットワークにおける隠れ変数ベクトル（１，１）の確率密度を示す図である。

【図３７】２つのニューラルネットワークの予測結果を示す図（その１）である。

【図３８】２つのニューラルネットワークの予測結果を示す図（その２）である。

【図３９】レジスタの個数と形式的情報基準の例を示す図である。

【符号の説明】

１ 入力手段 ２ ニューラルネットワーク手段 ３ 予測値生成手段 ４ 出力手段 １１ 制御装置 １２ 表示・対話装置 １３ システム管理装置 １４ 学習装置 １５ 観測装置 １６ ネットワーク学習装置 １７ 予測装置 １８ 非線形方程式求解装置 １９ コラム構造再帰型ニューラルネットワーク ２１ ＣＰＵ ２２ メモリ ２３ 入出力端末 ２４ バス ３１ アニーリング制御装置 ３２ パラメータ精錬装置 ３３ 評価関数算定装置 ４１ 零点算定装置 ４２ 予測値選択装置 ４３ 目標関数評価装置 ４４ 正規化係数算定装置 ４５ 確率密度算定装置 ５１−１，５１−２，５１−ｍ，５１−ｊ ニューラル素子（隠れ素子） ５２−１−１，５２−１−２，５２−１−ｓ，５２−２
−１，５２−２−２，５２−２−ｓ，５２−ｍ−１，５
２−ｍ−２，５２−ｍ−ｓ レジスタ ６１，７１，８１ 内部状態情報変換器 ６２−１，６２−２，６２−ｍ，６４，６５，７２−
１，７２−２，７２−ｍ，７７，８２−１，８２−２，
８２−ｍ 加算器 ６３−１，６３−２，６３−ｍ ロジスティック関数演算器 ７３ ２値ベクトル生成器 ７４ ベクトル乗算器 ７５ ベクトル加算器 ７６，８６ ノルム算定器 ７８，８３−１，８３−２，８３−ｍ，８５ 指数関数計算器 ７９ 累積型加算器 ８０ 逆数演算器 ８４ 乗算器