[go: up one dir, main page]

JP7163977B2 - 推定装置、学習装置、それらの方法、およびプログラム - Google Patents

推定装置、学習装置、それらの方法、およびプログラム Download PDF

Info

Publication number
JP7163977B2
JP7163977B2 JP2020569354A JP2020569354A JP7163977B2 JP 7163977 B2 JP7163977 B2 JP 7163977B2 JP 2020569354 A JP2020569354 A JP 2020569354A JP 2020569354 A JP2020569354 A JP 2020569354A JP 7163977 B2 JP7163977 B2 JP 7163977B2
Authority
JP
Japan
Prior art keywords
state
function
parameters
value
observable
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
JP2020569354A
Other languages
English (en)
Other versions
JPWO2020158032A1 (ja
Inventor
伸 村田
悠馬 小泉
登 原田
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
NTT Inc
NTT Inc USA
Original Assignee
Nippon Telegraph and Telephone Corp
NTT Inc USA
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nippon Telegraph and Telephone Corp, NTT Inc USA filed Critical Nippon Telegraph and Telephone Corp
Publication of JPWO2020158032A1 publication Critical patent/JPWO2020158032A1/ja
Application granted granted Critical
Publication of JP7163977B2 publication Critical patent/JP7163977B2/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/04—Architecture, e.g. interconnection topology
    • G06N3/047—Probabilistic or stochastic networks
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/04—Architecture, e.g. interconnection topology
    • G06N3/045—Combinations of networks
    • G06N3/0455—Auto-encoder networks; Encoder-decoder networks
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/04—Architecture, e.g. interconnection topology
    • G06N3/0475—Generative networks
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/08—Learning methods
    • G06N3/084—Backpropagation, e.g. using gradient descent
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/08—Learning methods
    • G06N3/09—Supervised learning

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Health & Medical Sciences (AREA)
  • Computing Systems (AREA)
  • Biomedical Technology (AREA)
  • Biophysics (AREA)
  • Computational Linguistics (AREA)
  • Data Mining & Analysis (AREA)
  • Evolutionary Computation (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Molecular Biology (AREA)
  • Artificial Intelligence (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Software Systems (AREA)
  • Health & Medical Sciences (AREA)
  • Probability & Statistics with Applications (AREA)
  • Testing And Monitoring For Control Systems (AREA)

Description

本発明は、状態空間モデルを用いて観測量から状態を推定する推定装置、学習装置、それらの方法、およびプログラムに関する。
系列データから対象の性質を解析するために、状態空間モデルと呼ばれる枠組みが広く用いられている。状態空間モデルは、観測できない隠れた「状態モデル」と観測した結果である「観測モデル」からなり、状態と呼ばれる量が時間発展し、それらの状態から観測過程を通じて観測量（たとえば電流・音圧・画像などの「量」）の系列データが生成されると考えるモデルである。
状態は非線形で時間推移（時間発展）するものであり、観測量は状態の時間推移に伴い非線形で時間推移する観測処理を行うことで得る(観測過程を通じて得る)ことができるものである。観測過程と時間発展の非線形性から、事前の仮定なしに状態空間モデルを観測量だけからすべて学習することは困難である。一方、近年研究が進んでいるKoopmanモード分解と呼ばれる手法は、状態空間モデルを別の領域（関数空間）で考えることにより、上述の非線形性を回避することが可能である(非特許文献１、２参照)。
Matthew O. Williams, Clarence W. Rowley, and Ioannis G. Kevrekidis., "A Kernel-Based Approach to Data-Driven Koopman Spectral Analysis", Journal of Computational Dynamic, 2:247-265, 2015. arXiv: 1411.2260. Matthew O. Williams, Ioannis G. Kevrekidis, and Clarence W. Rowley., "A Data-Driven Approximation of the Koopman Operator: Extending Dynamic Mode Decomposition", Journal of Nonlinear Science, 25(6):1307-1346, 2015.
しかしながら、状態が未知の場合にはKoopmanモード分解を適用できないため、依然として観測量だけから観測過程・時間発展・状態の全てを学習することはできない。
本発明は、観測量だけから観測過程・時間発展・状態を学習する学習装置、観測量だけから学習した観測過程・時間発展・状態を用いて、観測量から状態を推定する推定装置、それらの方法、およびプログラムを提供することを目的とする。
上記の課題を解決するために、本発明の一態様によれば、推定装置は、エンコーダを用いて、観測量から状態を推定する状態推定部と、デコーダとを用いて、状態から観測量を推定する観測量推定部と、時間展開を表すパラメータKを用いて、観測量が時間推移により変動した値である未来観測量を推定する未来観測量推定部と、を有し、エンコーダのパラメータと、デコーダのパラメータと、パラメータKとは、同時に最適化されたものである。
本発明によれば、観測量だけから観測過程・時間発展・状態を学習することができる。また、観測量だけから学習した観測過程・時間発展・状態を用いて、観測量から状態を推定することができるという効果を奏する。
第一実施形態の概要を示す図。 図２Ａは従来技術の空間状態モデルを説明するための図、図２Ｂは第一実施形態のオートエンコーダの枠組みを説明するための図。 第一実施形態に係る推定システムの構成例を示す図。 第一実施形態に係る学習装置の機能ブロック図。 第一実施形態に係る学習装置の処理フローの例を示す図。 学習段階のアルゴリズムを説明するための図。 第一実施形態に係る推定装置の機能ブロック図。 第一実施形態に係る推定装置の処理フローの例を示す図。 推定段階のアルゴリズムを説明するための図。 実際に画像データに基づくデータの系列データを入力としてパラメータを学習し、系列データの生成を行なった例を示す図。 状態に対応する観測量を予測するアルゴリズムを説明するための図。 第二実施形態に係る異常検知装置の機能ブロック図。 第二実施形態に係る異常検知装置の前処理の処理フローの例を示す図。 第二実施形態に係る異常検知装置の異常検知処理の処理フローの例を示す図。 第三実施形態の概要を示す図。 第三実施形態に係る学習装置の機能ブロック図。 第三実施形態に係る推定装置の推定部の機能ブロック図。 第三実施形態に係る推定装置の推定部の処理フローの例を示す図。
以下、本発明の実施形態について、説明する。なお、以下の説明に用いる図面では、同じ機能を持つ構成部や同じ処理を行うステップには同一の符号を記し、重複説明を省略する。以下の説明において、テキスト中で使用する記号「^」等は、本来直後の文字の真上に記載されるべきものであるが、テキスト記法の制限により、当該文字の直前に記載する。式中においてはこれらの記号は本来の位置に記述している。また、ベクトルや行列の各要素単位で行われる処理は、特に断りが無い限り、そのベクトルやその行列の全ての要素に対して適用されるものとする。
＜第一実施形態のポイント＞
本実施形態は、時間発展と観測過程だけでなく、同時に観測過程の「逆関数」を学習することで、観測量だけから観測過程・時間発展・状態の全てを学習することを可能とする。具体的にはオートエンコーダネットワーク(参考文献１参照)を用いる。
(参考文献１)G. E. Hinton. "Reducing the Dimensionality of Data with Neural Networks", Science, 313(5786):504-507, July 2006.
オートエンコーダはある入力をエンコーダと呼ばれるネットワークで低次元に変換し、デコーダと呼ばれるネットワークで復元する。ここで、エンコーダを「逆関数」、デコーダを観測過程とみなすことで、両方を同時に学習可能としたことが本実施形態の基礎技術である。上述したように、状態空間モデルは、時間経過で変化するかつ観測不能な状態から観測量を推定する観測モデルと、時間経過で変化する状態を推定する状態モデルからなる。つまり、状態空間モデルは状態から観測量を推定することが前提となる。しかしながら、前述したように状態が未知である場合にはそのようなモデルを構築することができないため、オートエンコーダの枠組みを利用し、エンコーダを観測量から状態を推定するモデルとし、デコーダを状態から観測量を推定するモデルとして学習することで、状態から観測量を推定するモデルの構築を可能とした。入出力を逆転させたモデルを学習するように構築したと言い換えてもよい。
以下で本実施形態の基礎技術の説明をする。
状態xtを考える。このとき状態は抽象・具体を問わない。状態は次のように時間発展する。
xt+1 =f(xt) (1)
この状態xtから、観測量ytが観測過程、
yt=g(xt) (2)
を通じて得られる。このとき、観測量は何らかの方法で数値化されるもの（たとえば電流電圧、温度、音圧、画像など）であり、マイクロホンアレイで得られるような多次元データでもよい。
状態空間モデルの枠組みにおける目標は、観測量の系列データ{y1,…,yT}が得られたときに、
・状態の系列データ{x1,…,xT}
・時間発展の関数f(xt)
・観測過程の関数g(xt)
を決定することである。
ただし、観測量だけから状態空間モデルをすべて決めるのは一般には難しい。これは、時間発展f(xt)と観測過程g(xt)が一般に非線形であり、事前知識なしに定めることが困難なことが一因である。
Koopmanモード分解は、前述の非線形性を回避可能な手法である。
（Koopmanモード分解）
基底関数
Figure 0007163977000001
で観測過程の関数を
Figure 0007163977000002
のように展開する。このとき、Koopman作用素の考え方を用いると、基底関数の時間発展は
Ψ(xt+1)=Ψ(f(xt))= KΨ(xt) (4)
のように書き直すことができる。
まとめると適切な変換zt=Ψ(xt)を用いて（無限次元の）関数空間上で議論をすることで、元の状態空間モデル
Figure 0007163977000003
が、線形の状態空間モデル
Figure 0007163977000004
と書き直せる。
前述の通り、状態空間モデルでは時間発展と観測過程の非線形性から、事前情報なしに生成モデルを学習することは困難である。一方、Koopmanモード分解は、上記の非線形性を、関数空間上で議論することで回避することが可能である。しかし、Koopmanモード分解は、状態xtが既知な場合に適用可能である。そのため、Koopmanモード分解を用いても、観測量の系列データだけから、時間発展・観測過程・状態を学習することは困難である。
そこで本実施形態では、Koopmanモード分解の枠組みに基づく、新たな状態推定手法を考案する。本実施形態を用いることで、観測量の系列データだけから、時間発展・観測過程・状態のすべてを同時に学習することが可能となる。
（基本的なアイデア）
図１に本実施形態の概要を記載している。観測量の系列データ{yt}を入力とし、状態xt、観測過程(^Ψ,B)、時間発展Kを学習する。学習後はある観測量ytから予測した観測量の系列データ^y1 (t),^y2 (t),…や状態xtを出力する。
状態xtから観測量^ytは
Figure 0007163977000005
と生成されるので、現在の観測量ytから状態xtを推定するには上記の逆問題を解く必要がある。すなわち、
Figure 0007163977000006
である。B-1(・)はBの擬似逆行列や、Ridge回帰で解析的に求めることができる。一方で、基底関数Ψ(・)の逆関数 Ψ-1(・)は一般に求めることは困難である。ここで、逆関数・基底関数・状態を決めるためにオートエンコーダネットワークの枠組みを利用した実装例を説明する。オートエンコーダネットワークはデータの次元削減に使われるニューラルネットワークであり、エンコーダを通して入力が中間層に写され、デコーダを通じて復元される。すなわち、
Figure 0007163977000007
である。ここで、エンコーダΦ-1を基底関数の逆関数Ψ-1(・)、デコーダΦを基底関数Ψ(・)とみなす。
別の言い方をすると、従来技術では、状態xtを基底で^ztに変換し、^ztと係数Bから予測した観測量^ytを導くような空間状態モデルになるよう、時間経過に応じて線形変換されるような基底を得ていた(図１、図２Ａ参照)。一方、本実施形態では、観測量ytを入力し状態xtを出力するエンコーダ、状態xtを入力し観測量^ytを出力するデコーダ、となるようオートエンコーダの枠組みを用いて学習する(図１、図２Ｂ参照)。
＜第一実施形態＞
本実施形態は二段階に分かれている。一つは観測量の時系列データから時間発展と観測過程ならびにその逆関数を学習する段階（以下、学習段階ともいう）である。一つは観測量から状態を取得する段階（以下、状態取得段階ともいう）である。
本実施形態に係る状態取得システムは、学習装置１００と状態取得装置２００とを含む（図３参照）。なお、状態取得装置２００は、観測量から状態、または、状態から観測量を推定するので、推定装置２００ともいう。同様に、状態取得段階を推定段階、状態取得システムを推定システムともいう。
学習装置１００により前述の学習段階を実行し、推定装置２００により前述の推定段階を実行する。まず、学習段階について説明する。
（学習段階）
図４は学習装置１００の機能ブロック図を、図５はその処理フローの例を示す。学習装置１００は、初期化部１１０、推定部１２０、目的関数計算部１３０、パラメータ更新部１４０とを含む。
学習装置１００は、学習用の観測量の系列データ{yt (L)}を入力とし、学習結果であるパラメータ(wenc,wdec,K,B)を出力する。wencはエンコーダで用いられるパラメータ(逆関数Ψ-1のパラメータ)を、wdecはデコーダで用いられるパラメータ(基底関数Ψのパラメータ)を、Kは時間発展を表すパラメータ、Bは展開係数を表す。
学習段階は図６のアルゴリズム1のように行われる。学習用の観測量が画像データに基づくデータである例を示す。ここでは、入力データとして画像データに基づくデータからなる系列データ{yt}（e.g.動画データ）を考える。画像データに基づくデータytは、例えば、各画素の画素値からなる、画素数だけの次元を持つデータでもよいし、画像データから得られる特徴量や動画データから得られる特徴量等でもよい。画像データに基づくデータを入力データとする場合、状態は、観測量ytに対応する抽象的な状態を表し、この抽象的な状態が物理的な意味を持つ値(例えば、観測量に比べ低次元の値)となり得る。例えば、状態は、画像データに表れる周期的な模様や物体の大きさや位置に類する量、画像データからなる動画データ内での動体の位相に類する量(周期的に動作する場合)を表す。
＜初期化部１１０＞
初期化部１１０は、学習に先立ち、推論に用いるパラメータwenc (k),wdec (k),K(k),B(k)を初期化し(Ｓ１１０)、初期化したパラメータwenc (0),wdec (0),K(0),B(0)を推定部１２０に出力する。また、更新回数を示すインデックスkをk=0とする。
＜推定部１２０＞
推定部１２０は、学習用の観測量の系列データ{yt (L)}と初期化または更新したパラメータwenc (k),wdec (k),K(k),B(k)を入力とし、(1)基底関数の値の推定(推定値zt)、(2)状態の推定(推定値xt)、(3)再構成した基底関数の値の推定(推定値^zt)、(4)基底関数の予測(予測値^zτ (t))、(5)観測量の予測(予測値^yτ (t))を行い(Ｓ１２０)、推定値zt、xt、^zt、予測値^zτ (t)、^yτ (t)を出力する。以下、推定および予測を合わせて推定ともいう。以下、推定の詳細を説明する。
(1)まず、推定部１２０は、現在の展開係数B(k)を用いて基底関数の値ztを推定する。すなわち、zt=B(k)+(yt (L))である。ここで、B(k)+(・)は回帰問題を解くことを表す。本実施形態ではRidge回帰の問題としてB(k)+(・)=(B(k)TB(k)+σI)-1B(k)T・とするが、擬似逆行列を用いた通常の線形回帰やLASSOなどのスパース推定アルゴリズムを用いてもよい。ただし、σはRidge回帰における所定の重みパラメータである。
(2)次に、式(9)のニューラルネットワークを用いて、状態の推定をxt=Φ-1(zt;wenc (k))として行う。
(3)更に、式(9)のニューラルネットワークを用いて、再構成された基底関数の値の推定を^zt=Φ(xt;wdec (k))として行う。
(4)次にτ=0,…,Τとして、基底関数の値の線形予測を行い、基底関数の値の予測値^zτ (t)=K(k)τ^ztを得る。
(5)予測値^zτ (t)をもとに、観測量の予測値^yτ (t) = B(k)^zτ (t)を得る。
＜目的関数計算部１３０＞
目的関数計算部１３０は、学習用の観測量の系列データ{yt+τ (L)}と、推定値の系列データ{zt+τ}、{xt+τ}、{^zt+τ}、予測値の系列データ{^zτ (t)}、{^yτ (t)}、パラメータwenc (k)、wdec (k)を入力とし、目的関数の値J(Θ)を計算し(Ｓ１３０)、出力する。なお、Θはパラメータwenc (k),wdec (k),K(k),B(k)の集合である。
(i)まず、観測量の予測誤差を次式により求める。
Figure 0007163977000008
ただし、δτは0<δτ<1であり、現時点tから近い時刻t+τの誤差ほど大きく評価するように設定するための重みパラメータであり、E[・]は期待値計算を表す。
(ii)また、基底関数の予測誤差を次式により求める。
Figure 0007163977000009
(iii)また、正則化項として、ニューラルネットワークの重みに対する正則化項Ω1を導入する。なお、Ω1=||wenc (k)||2 2+||wdec (k)||2 2である。
(iv)また、状態に構造を導入する。本実施形態では滑らかさΩ2=E[||xt+1-xt||2 2]と非ガウス性Ω3=E[log cosh(xt)]を導入する。ただし、cosh(・)は双曲線関数を表す。
(v)以上の項を重みづけした目的関数J(Θ)=aJ1+bJ2+p1Ω1+p2Ω2+p3Ω3の値を求める。ただし、a,b,p1,p2,p3は、それぞれJ1,J2,Ω1,Ω2,Ω3をどの程度重視するかを決定するためのパラメータであり、実験結果やシミュレーション結果等を用いて適宜設定する。
＜パラメータ更新部１４０＞
パラメータ更新部１４０は、目的関数J(Θ)を受け取り、各パラメータwenc (k),wdec (k),K(k),B(k)を更新する（Ｓ１４０）。例えば、Back propagationを用いて各パラメータに対する勾配ΔΘJを計算し、各パラメータwenc (k),wdec (k),K(k),B(k)を更新する(Θ(k+1)=Θ(k)+ΔΘJ)。
所定の条件（例えば、所定の回数のループを終えるか、目的関数が変化しなくなる等）を満たさない場合には(Ｓ１４１のnoの場合)、更新したパラメータwenc (k+1),wdec (k+1),K(k+1),B(k+1)を推定部１２０に、正則化項Ω1算出用に更新したパラメータwenc (k+1),wdec (k+1)を目的関数計算部１３０に出力し、k←k+1とし、Ｓ１２０～Ｓ１４０を繰り返す。
所定の条件を満たす場合には(Ｓ１４１のyesの場合)、パラメータの更新を打ち切り、モデルの学習を終える。最新の各パラメータwenc (k),wdec (k),K(k),B(k)を学習結果であるパラメータ(wenc,wdec,K,B)として、推定装置２００に出力する。
(推定段階)
図７は推定装置２００の機能ブロック図を、図８はその処理フローの例を示す。推定装置２００は、推定部２２０を含む。
なお、推定装置２００は、推定、予測に先立ち、学習装置１００から受け取ったパラメータ(wenc,wdec,K,B)を推定部２２０に設定する。
(アルゴリズム2)
図９のアルゴリズム2を実行する場合、推定装置２００は、観測量ytを入力とし、観測量ytに対応する状態を推定し、yt以降の観測量の系列データを予測し、推定値xt、予測した系列データ{^yτ (t)}を出力する。図９のアルゴリズム2では、推定部２２０に適当な画像データに基づくデータを観測量ytとして、Τステップ先の画像データまで予測することができる。例えば、観測量の系列データyt,yt+1,…,yt+Nを入力とし、推定値の系列データxt,xt+1,…,xt+Nと、予測したN個の系列データ{^yτ (t)},{^yτ (t+1)},…,{^yτ (t+N)}を出力してもよい。
＜推定部２２０＞
推定部２２０は、観測量ytを入力とし、観測量ytに対して所定の処理を行うことで状態を推定する（Ｓ２２０）。本実施形態では、所定の処理は、(1)基底関数の値の推定(推定値zt)、(2)状態の推定(推定値xt)である。さらに、推定部２２０は、(3)再構成した基底関数の値の推定(推定値^zt)、(4)基底関数の予測(予測値^zτ (t))、(5)観測量の予測(予測値^yτ (t))を行い(Ｓ２２０)、推定値xt、予測値の系列データ{^yτ (t)}を出力する。
(1)まず、推定部２２０は、展開係数Bを用いて基底関数の値ztを推定する。すなわち、zt=B+(yt)である。
(2)次に、式(9)のニューラルネットワークを用いて、状態の推定をxt=Φ-1(zt;wenc)として行う。
(3)更に、式(9)のニューラルネットワークを用いて、再構成された基底関数の値の推定を^zt=Φ(xt;wdec)として行う。
(4)次にτ=0,…,Τとして、基底関数の値の線形予測を行い基底関数の値の予測値^zτ (t)=Kτ^ztを得る。
(5)予測値^zτ (t)をもとに、観測量の予測値^yτ (t) = B^zτ (t)を得る。
なお、上述の(3)～(5)の処理は、状態から観測量を得る処理であり、上述の(1),(2)の処理の逆処理と言える。
図１０は、実際に画像データに基づくデータの系列データを入力としてパラメータを学習し、系列データの生成を行なった例を示す。
図１０の上段は実際の画像データに基づくデータの系列データ(学習用の観測量の系列データ)である。一方で下段は、上段左端"Input"として与えられた画像データをアルゴリズム2のytとして入力し、τ=1,…,10だけ生成した系列データ(予測した系列データ)である。
（アルゴリズム3）
図１１のアルゴリズム3を実行する場合、推定装置２００は、状態xtを入力とし、状態xtに対応する観測量を予測し、予測した系列データ{^yτ (t)}を出力する。図１１のアルゴリズム3では、推定部２２０に適当な状態xtを与え、Τステップ先の画像まで予測することができる。
推定部２２０は、ある状態xtを入力とし、(1)再構成した基底関数の値の推定(推定値^zt)、(2)基底関数の予測(予測値^zτ (t))、(3)観測量の予測(予測値^yτ (t))を行い(Ｓ２２０)、予測値の系列データ{^yτ (t)}を出力する。
(1)式(9)のニューラルネットワークを用いて、再構成された基底関数の値の推定を^zt=Φ(xt;wdec)として行う。
(2)次にτ=0,…,Τとして、基底関数の値の線形予測を行い、基底関数の値の予測値^zτ (t)=Kτ^ztを得る。
(3)予測値^zτ (t)をもとに、観測量の予測値^yτ (t) = B^zτ (t)を得る。
＜効果＞
学習段階では、観測量だけから観測過程・時間発展・状態を学習することができる。
推定段階では、観測量だけから学習した観測過程・時間発展・状態（学習したモデル）を用いて、観測量から状態を推定することができる。また、推定した状態や与えられた状態から観測量を予測することができる。すなわち、現在の観測量から状態を推定し、時間発展をシミュレートし、その状態を観測することで、系列データを予測できる。また、観測量や状態を与えることで、人工的にデータ（状態、観測量）を生成することができる。
センサなどで観測される系列データから状態を推定することができ、観測量の分析等に用いることができる。また、視覚的に把握しずらい(例えば次元数が多い)系列データから、(例えば次元数が少ない)状態を推定し、推定した状態を提示することで観測量を視覚的に把握しやすいものに変換する（可視化）ことができる。
さらに予測した値と実際の観測量との距離を適切に定義することで、系列データの異常検知へと応用することも可能である。
＜変形例＞
本実施形態では、観測量が画像データに基づくデータの場合について説明したが、他のデータであってもよい。例えば、音響データに基づくデータ、振動データに基づくデータ、音響データに基づくデータと振動データに基づくデータとの組合せなどが考えられる。以下に、より詳しく説明する。
（音響データ）
観測量を音響データに基づくデータとした場合、マイクから取れた音圧の波形データまたはその特徴量（STFT やlog-Mel power など）を入力ytとする。なお、マイクロホンアレイで収音する場合、音圧の波形データまたはその特徴量を素子数分だけ結合したベクトルを入力ytとする。
この場合、状態は、観測量ytに対応する抽象的な状態を表し、この抽象的な状態が物理的な意味を持つ値(例えば、観測量に比べ低次元の値)となり得る。例えば、状態は、音源の波形に類する量、音源位置に類する量（音源が動く場合）、位相に類する量（周期的な音の場合）を表す。
（振動データ）
観測量を振動データに基づくデータとした場合、振動ピックアップから取れた加速度の波形データまたはその特徴量（STFT やlog-Melpower など）を入力ytとする。複数の振動ピックアップでの収録する場合、波形データや特徴量を素子数分だけ結合したベクトルを入力とする。
この場合、状態は、観測量ytに対応する抽象的な状態を表し、この抽象的な状態が物理的な意味を持つ値(例えば、観測量に比べ低次元の値)となり得る。例えば、状態は、振動のモードに類する量、位相に類する量（準周期的に動く物体の振動の場合）などを表す。
（音響・振動データ）
観測量を上述の音響データに基づくデータと振動データに基づくデータとの組合せとした場合、振動ピックアップから取れた加速度の波形データまたはその特徴量とマイクから取れた音圧の波形データまたはその特徴量を素子数分だけ結合したベクトルを入力ytとする。
この場合、状態xtは、観測量ytに対応する抽象的な状態を表し、この抽象的な状態が物理的な意味を持つ値(例えば、観測量に比べ低次元の値)となり得る。例えば、状態は、音源（振動源）の波形に類する量、振動のモードに類する量などを表す。
＜第二実施形態＞
第一実施形態と異なる部分を中心に説明する。
本実施形態では、本発明を異常検知に適用する。
図１２は本実施形態に係る異常検知装置の機能ブロック図をを示す。
まず、第一実施形態で説明した方法により、学習装置１００は、学習用の観測量の系列データ{yt (L)}を入力とし、パラメータ(wenc,wdec,K,B)を学習し、出力する。
第一実施形態で説明した方法により、推定装置２００は、推定、予測に先立ち、学習装置１００から受け取ったパラメータ(wenc,wdec,K,B)を推定部２２０に設定する。
＜異常検知装置３００＞
異常検知装置３００は、誤差ベクトル計算部３１０と平均分散共分散行列計算部３２０と検知部３３０とを含む（図１２参照）。
異常検知装置３００は、異常検知処理と、異常検知処理の前に予めパラメータを求める前処理とを実行する。まず、前処理について説明する。
＜前処理＞
まず、正常時の観測量のデータセットDnormal={y1,y2,…,yT_1}を用意する。ただし、下付き添え字A_BはABを意味する。ここから，系列の長さLの部分系列Dt={yt+1,…,yt+L}、t=1,2,…,T1-Lを取り出す。
次に、第一実施形態で説明した方法により、推定装置２００は、観測量ytを入力とし、yt以降の観測量の系列データを予測し、長さLの予測した系列データPt={^y(t) 1,…,^y(t) L}を出力する。
図１３は前処理の処理フローの例を示す。
異常検知処理に先立ち、異常検知装置３００は、正常時の観測量のデータセットDnormalから得られるT1-L個の部分系列DtとT1-L個の予測した系列データPtとを入力とし、後述する平均μと分散共分散行列Sとを計算しておく。
誤差ベクトル計算部３１０は、T1-L個の部分系列Dt={yt+1,…,yt+L}とT1-L個の系列データPt={^y(t) 1,…,^y(t) L}とを入力とし、部分系列Dt={yt+1,…,yt+L}と予測した系列データPt={^y(t) 1,…,^y(t) L}の各要素で誤差を計算し（Ｓ３１０－Ａ）、誤差ベクトルet=[e(t) 1…e(t) L]Tを出力する。なお、e(t) i=yt+i-^y(t) i、i=1,…,Lである。観測量ytを入力としたとき(ここでは、観測量ytから予測した系列データPt={^y(t) 1,…,^y(t) L}と系列データPtに対応する部分系列Dt={yt+1,…,yt+L}を入力としたとき)の誤差ベクトルをet=[e(t) 1…e(t) L]Tと定める。このとき、誤差ベクトルetは観測量ytの次元がDのとき(D×L)の長さのベクトルである。上記の誤差ベクトルの計算をT1-L個の部分系列Dt及び系列データPtすべてに対して計算する。
平均分散共分散行列計算部３２０は、T1-L個の誤差ベクトルet=[e(t) 1…e(t) L]Tを入力とし、平均μと分散共分散行列Sを次式により計算し（Ｓ３２０）、検知部３３０に出力する。
Figure 0007163977000010
Figure 0007163977000011
＜異常検知処理＞
異常検知の対象となる観測量のデータセットDnew={y'1,y'2,…,yT_2}が得られたとき、前処理のときと同様に、系列の長さLの部分系列D't'={yt'+1,…,yt'+L}、t'=1,2,…,T2-Lを取り出す。
次に、第一実施形態で説明した方法により、推定装置２００は、観測量y't'を入力とし、y't'以降の観測量の系列データを予測し、長さLの予測した系列データP't'={^y'(t') 1,…,^y'(t') L}を出力する。
図１４は異常検知処理の処理フローの例を示す。
異常検知処理時には、異常検知装置３００は、異常検知の対象となる観測量のデータセットDnewから得られるT2-L個の部分系列D't'とT2-L個の予測した系列データP't'={^y'(t') 1,…,^y'(t') L}とを入力とし、部分系列D't'と系列データP't'に対する検知結果を出力する。なお、系列データP't'は観測量y't'から予測したものであり、部分系列D't'は系列データP't'に対応するものなので、部分系列D't'と系列データP't'に対する検知結果は、観測量y't'に対する検知結果といってもよい。
誤差ベクトル計算部３１０は、T2-L個の部分系列D't'={y't'+1,…,yt'+L}とT2-L個の系列データPt'={^y'(t') 1,…,^y'(t') L}とを入力とし、部分系列D't'={y't'+1,…,yt'+L}と予測した系列データPt'={^y'(t') 1,…,^y'(t') L}の各要素で誤差を計算し（Ｓ３１０－Ｂ）、誤差ベクトルe't=[e'(t') 1…e'(t') L]Tを出力する。なお、e'(t') i=y't'+i-^y'(t') i、i=1,…,Lである。
検知部３３０は、異常検知に先立ち平均μと分散共分散行列Sを受け取っておく。
検知部３３０は、T2-L個の誤差ベクトルe't=[e'(t') 1…e'(t') L]Tを入力とし、平均μと分散共分散行列Sと誤差ベクトルe't'とを用いて、各時刻t'=1,…,T2-Lにおける以下の異常度合いLt'を計算する（Ｓ３３０－１）。
Lt'=logdet(S)+(e't'-μ)TS-1(e't'-μ)T
この異常度合いLt'は誤差ベクトルを正規分布でフィッティングしたときの負の対数尤度に比例する量である。
次に、検知部３３０は、異常度合いLt'に対応する値と閾値pとの大小関係に基づき、異常か否かを判定し、異常を検知し（Ｓ３３０－２）、検知結果を出力する。例えば、異常度合いLt'>pのとき異常と判定し、Lt'≦pのとき異常ではないと判定する。閾値pは、実験やシミュレーション等により予め適切にさだめておく。
以上の構成により、本発明を異常検知に適用することができる。
＜第三実施形態＞
第一実施形態と異なる部分を中心に説明する。
第一実施形態では、状態xtから観測量の予測(予測値^yτ (t))を出力する関数に逆関数が存在する場合を検討している。しかし、例えば、状態xtが観測量の予測(予測値^yτ (t))より低次元な場合、観測量ytから状態xtを推定する逆関数が存在するかは定かではない。
また、基底関数の値ztの次元より観測量ytの次元が少ない場合、観測量ytから基底関数の値ztを定める問題は劣決定問題であり、適切に正則化項を導入する必要がある。
そこで、本実施形態では、生成モデルを考えることで、逆関数を考えることなく観測量から状態を推定する。図１５に提案するモデルを示す。図１５の右(破線部分)がKoopmanモード分解の数値計算手法である拡張動的モード分解（EDMD）部分を表す。
まず、観測量ytとyt+1からニューラルネットワークを用いて平均と分散パラメータを推定する((σt,μt)←~Ψ(yt,yt+1))。この処理が変分オートエンコーダのエンコードに相当し、この処理を実現する部分をエンコーダと呼ぶ。
次に、得られた平均の推定値μtと分散パラメータの推定値σtに従う多変量正規分布から潜在変数(状態xt)をサンプリングする。なお、図中のetは平均0,分散1の正規分布から得られる乱数である。
以降の処理は第一実施形態と同様である。以下に概略を説明する。
(3)再構成した基底関数の値の推定(推定値^zt←Ψ(xt))、(4)基底関数の予測(予測値^zτ (t)=Kτ^zt)、(5)観測量の予測(予測値^yτ (t) = B^zτ (t)))を行い、推定値xt、予測値の系列データ{^yτ (t)}を出力する。状態xtから観測量の予測値^yτ (t)を得る処理が変分オートエンコーダのデコードに相当し、この処理を実現する部分をデコーダと呼ぶ。
通常の変分オートエンコーダは、ニューラルネットワークの重みパラメータθを次の目的関数を最小化するように学習する。
Figure 0007163977000012
ただし、KL[A|B]は分布A、BのKullback-Leibler divergenceを表し、N(μ,σ)は平均μ、分散σの分布を表し、μθ(yt)、Σθ(yt)はそれぞれ重みパラメータθのニューラルネットワークに観測量yt+τを与えて推定した平均、分散パラメータを表す。
一方、本実施形態の変分オートエンコーダでは展開係数B、パラメータKも重みパラメータθと同時に学習し、最適化する。すなわち、
Figure 0007163977000013
を最小化する。ただし、^yτ (t) = BKτ^ztであり、μθ(yt+τ,yt+τ+1)、Σθ(yt+τ,yt+τ+1)はそれぞれ重みパラメータθのニューラルネットワークに観測量yt+τ,yt+τ+1を与えて推定した平均、分散パラメータを表す。第一実施形態と本実施形態との差分としては次の二点が挙げられる。
第一実施形態は逆関数を仮定していたが、本実施形態は確率的な生成モデルを仮定している。
第一実施形態では再構成誤差を目的関数として最小化していたが、本実施形態は再構成誤差に分布の近さを測るKullback-Leibler divergence 項が追加されている。
以下、本実施形態を実現する推定システムについて説明する。
本実施形態に係る推定システムは、学習装置１００と推定装置２００とを含む（図３参照）。
学習装置１００により学習段階を実行し、推定装置２００により推定段階を実行する。まず、学習段階について説明する。
（学習段階）
図１６は学習装置１００の機能ブロック図を、図５はその処理フローの例を示す。学習装置１００は、初期化部１１０、推定部１２０、目的関数計算部１３０、パラメータ更新部１４０とを含む。
学習装置１００は、学習用の観測量の系列データ{yt}を入力とし、学習結果であるパラメータ(wenc,wdec,K,B)を出力する。wencはエンコーダで用いられるパラメータを、wdecはデコーダで用いられるパラメータ(基底関数Ψのパラメータ)を、Kは時間発展を表すパラメータ、Bは展開係数を表す。
学習段階は以下のように行われる。
初期化部１１０およびパラメータ更新部１４０の処理内容は第一実施形態と同様なので、説明を省略する。ただし、パラメータ更新部１４０では、目的関数J(Θ(k))に代えて、目的関数L(B(k),K(k),θ(k))を受け取り、処理を行う。
＜推定部１２０＞
推定部１２０は、学習用の観測量の系列データ{yt}と初期化または更新したパラメータwenc (k),wdec (k),K(k),B(k)を入力とし、(1)状態の平均と分散パラメータの推定(推定値σt,μt)、(2)状態の推定(推定値xt)、(3)再構成した基底関数の値の推定(推定値^zt)、(4)基底関数の予測(予測値^zτ (t))、(5)観測量の予測(予測値^yτ (t))を行い(Ｓ１２０)、予測値^yτ (t)を出力する。(3)～(5)は第一実施形態と同様なので、(1),(2)について説明する。
(1)まず、推定部１２０は、学習用の観測量の系列データ{yt}と現在のパラメータwenc (k)とからニューラルネットワークを用いて、状態の平均と分散パラメータとを推定する。すなわち、(σt,μt)=~Ψ(yt;wenc (k))である。なお、本実施形態では、ニューラルネットワークを用いるため、入力を2つ以上の学習用の観測量ytとしてもよい。例えば、(σt,μt)=~Ψ(yt,yt+1;wenc (k))として、2つの学習用の観測量yt,yt+1から状態の平均と分散パラメータとを推定してもよい。このように、2つ以上の学習用の観測量を用いて状態の平均と分散パラメータとを推定することで、状態の特徴をより捉えやすくなると考えられる。
(2)次に、(σt,μt)に従う多変量正規分布から状態xtをサンプリングする。
＜目的関数計算部１３０＞
目的関数計算部１３０は、学習用の観測量の系列データ{yt+τ}と、予測値の系列データ{^yτ (t)}、パラメータwenc (k)を入力とし、目的関数の値L(B(k),K(k),θ(k))を計算し(Ｓ１３０)、出力する。なお、θ(k)はパラメータwenc (k),wdec (k)の集合である。
Figure 0007163977000014
ただし、^yτ (t) = B(k)K(k)τ^ztである。なお、ここでは、K(k)を同時に更新するために、すくなくとも2つ以上の時系列の観測量yτとyτ+1とを必要とし、Τは１以上の整数の何れかとする。c,dは、各項をどの程度重視するかを決定するためのパラメータであり、実験結果やシミュレーション結果等を用いて適宜設定する。
(推定段階)
図７は推定装置２００の機能ブロック図を、図８はその処理フローの例を示す。推定装置２００は、推定部２２０を含む。
なお、推定装置２００は、推定、予測に先立ち、学習装置１００から受け取ったパラメータ(wenc,wdec,K,B)を推定部２２０に設定する。
(観測量から状態を推定し、観測量を予測)
観測量から状態を推定し、観測量を予測する場合(第一実施形態のアルゴリズム2の場合)、推定装置２００は、観測量ytを入力とし、観測量ytに対応する状態を推定し、yt以降の観測量の系列データを予測し、推定値xt、予測した系列データ{^yτ (t)}を出力する。例えば、観測量の系列データyt,yt+1,…,yt+Nを入力とし、推定値の系列データxt,xt+1,…,xt+Nと、予測したN個の系列データ{^yτ (t)},{^yτ (t+1)},…,{^yτ (t+N)}を出力してもよい。
＜推定部２２０＞
推定部２２０は、観測量ytを入力とし、観測量ytに対して所定の処理を行うことで状態を推定する（Ｓ２２０）。本実施形態では、所定の処理は、状態の平均と分散パラメータの推定(推定値σt,μt)、(2)状態の推定(推定値xt)である。さらに、推定部２２０は、(3)再構成した基底関数の値の推定(推定値^zt)、(4)基底関数の予測(予測値^zτ (t))、(5)観測量の予測(予測値^yτ (t))を行い(Ｓ２２０)、推定値xt、予測値の系列データ{^yτ (t)}を出力する。(3)～(5)は第一実施形態と同様なので、(1),(2)について説明する。
(1)まず、推定部２２０は、観測量ytとパラメータwencとからニューラルネットワークを用いて、状態の平均と分散パラメータとを推定する。すなわち、(σt,μt)=~Ψ(yt;wenc)である。なお、本実施形態では、入力を2つ以上の観測量ytとしてもよく、学習装置１００で学習したニューラルネットワークに対応する個数の観測量ytを入力とする。
(2)次に、(σt,μt)に従う多変量正規分布から状態xtをサンプリングする。
（観測量から観測量を予測）
観測量から観測量を予測する場合(第一実施形態のアルゴリズム3の場合)、第一実施形態と同様である。
＜効果＞
第一実施形態と同様の効果を得ることができる。生成モデルを考えることで、逆関数を考えることなく状態を推定する。なお、本実施形態と第一実施形態の変形例や第二実施形態を組合せてもよい。
＜その他＞
第一実施形態および第三実施形態の推定部２２０を図１７の機能ブロック図で表現することもできる。図１８は推定部２２０の処理フローの例を示す。
推定部２２０は、状態推定部２２１と、観測量推定部２２２と、未来観測量推定部２２３とを有する。さらに観測量推定部２２２は、中間値推定部２２２Ａと、中間観測値推定部２２２Ｂとを含む。状態推定部２２１では、第一実施形態と第三実施形態とにおける処理内容が異なる。観測量推定部２２２および未来観測量推定部２２３では、第一実施形態と第三実施形態とにおける処理内容が同じである。
（状態推定部２２１）
状態推定部２２１は、オートエンコーダのエンコーダを用いて、観測量から状態を推定し（Ｓ２２１）、出力する。
第一実施形態では、状態推定部２２１は、推定処理に先立ちニューラルネットワークのパラメータwenc,展開係数Bを受け取る。状態推定部２２１は、観測量ytを入力とし、展開係数Bを用いて基底関数の値ztを推定する。すなわち、zt=B+(yt)である。さらに、式(9)のニューラルネットワークを用いて、状態の推定をxt=Φ-1(zt;wenc)として行う。
第三実施形態では、状態推定部２２１は、推定処理に先立ちニューラルネットワークのパラメータwencを受け取る。状態推定部２２１は、１つ以上の観測量ytを入力とし、１つ以上の観測量ytとパラメータwencとからニューラルネットワークを用いて、状態の平均と分散パラメータとを推定する。例えば、2つの観測量yt、yt+1から(σt,μt)=~Ψ(yt,yt+1;wenc)として、状態の平均と分散パラメータとを推定する。次に、状態推定部２２１は、(σt,μt)に従う多変量正規分布からサンプリングすることで状態xtを推定する。
（観測量推定部２２２）
観測量推定部２２２は、オートエンコーダのデコーダとを用いて、状態から観測量を推定し（Ｓ２２２）、出力する。
なお、第一実施形態の場合、状態推定部２２１が行う処理を第一の関数で、観測量推定部２２２が行う処理を第二の関数で定義され、第一の関数は第二の関数の逆関数である。
（中間値推定部２２２Ａ）
中間値推定部２２２Ａは、推定処理に先立ちニューラルネットワークのパラメータwdecを受け取る。
中間値推定部２２２Ａは、状態xtを入力とし、式(9)のニューラルネットワークを用いて、再構成された基底関数の値の推定を^zt=Φ(xt;wdec)として行い（Ｓ２２２Ａ）、推定値^ztを出力する。なお、再構成された基底関数の値の推定値を中間値ともいう。
（中間観測値推定部２２２Ｂ）
中間観測値推定部２２２Ｂは、推定処理に先立ち展開係数Bを受け取る。
中間観測値推定部２２２Ｂは、推定値^ztを入力とし、推定値^ztから、観測値を推定し（Ｓ２２２Ｂ）、推定値^ytを出力する。以下の式において、τ=0に対応する。
^zτ (t)=Kτ^zt
^yτ (t) = B^zτ (t)
つまり、^yt = B^ztである。
（未来観測量推定部２２３）
未来観測量推定部２２３は、推定処理に先立ちK、Bを受け取る。
未来観測量推定部２２３は、時間展開を表すパラメータKを用いて、観測量が時間推移により変動した値である未来観測量を推定し（Ｓ２２３）、出力する。
まず、未来観測量推定部２２３は、パラメータKを用いて、推定値^ztが時間推移により変動した値である未来中間値^zτ (t)を得る。つまり、^zτ (t)=Kτ^ztである。ただし、τ=1,…,Τである。
さらに、未来観測量推定部２２３は、展開係数Bを用いて、未来中間値^zτ (t)から未来観測量^yτ (t)を推定する。つまり、^yτ (t) = B^zτ (t)である。
なお、第一実施形態および第三実施形態の学習装置の推定部１２０についても同様に表現することができる。ただし、学習用の観測量と、学習対象のパラメータを用いて処理を行う。
＜ハードウェア構成＞
学習装置および推定装置は、例えば、中央演算処理装置（CPU: Central Processing Unit）、主記憶装置（RAM: Random Access Memory）などを有する公知又は専用のコンピュータに特別なプログラムが読み込まれて構成された特別な装置である。学習装置および推定装置は、例えば、中央演算処理装置の制御のもとで各処理を実行する。学習装置および推定装置に入力されたデータや各処理で得られたデータは、例えば、主記憶装置に格納され、主記憶装置に格納されたデータは必要に応じて中央演算処理装置へ読み出されて他の処理に利用される。学習装置および推定装置の各処理部は、少なくとも一部が集積回路等のハードウェアによって構成されていてもよい。学習装置および推定装置が備える各記憶部は、例えば、RAM（Random Access Memory）などの主記憶装置、またはリレーショナルデータベースやキーバリューストアなどのミドルウェアにより構成することができる。ただし、各記憶部は、必ずしも学習装置および推定装置がその内部に備える必要はなく、ハードディスクや光ディスクもしくはフラッシュメモリ（Flash Memory）のような半導体メモリ素子により構成される補助記憶装置により構成し、学習装置および推定装置の外部に備える構成としてもよい。
＜その他の変形例＞
本発明は上記の実施形態及び変形例に限定されるものではない。例えば、上述の各種の処理は、記載に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。その他、本発明の趣旨を逸脱しない範囲で適宜変更が可能である。
＜プログラム及び記録媒体＞
また、上記の実施形態及び変形例で説明した各装置における各種の処理機能をコンピュータによって実現してもよい。その場合、各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記各装置における各種の処理機能がコンピュータ上で実現される。
この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。
また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ－ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させてもよい。
このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶部に格納する。そして、処理の実行時、このコンピュータは、自己の記憶部に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。また、このプログラムの別の実施形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよい。さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。なお、プログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。
また、コンピュータ上で所定のプログラムを実行させることにより、各装置を構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。

Claims (6)

  1. エンコーダを用いて、観測量から状態を推定する状態推定部と、
    デコーダとを用いて、前記状態推定部で推定した前記状態から観測量を推定する観測量推定部と、
    時間発展を表すパラメータKを用いて、前記観測量が時間推移により変動した値である未来観測量を推定する未来観測量推定部と、を有し、
    前記エンコーダのパラメータと、前記デコーダのパラメータと、前記パラメータKとは、同時に最適化されたものであり、
    前記状態推定部が行う処理は第一の関数で、前記観測量推定部が行う処理は第二の関数で定義され、
    前記第一の関数は前記第二の関数の逆関数であり、
    前記第二の関数はオートエンコーダのデコーダのパラメータである、基底関数のパラメータを用い、前記第一の関数は前記オートエンコーダのエンコーダのパラメータである、前記基底関数の逆関数のパラメータを用い、
    前記未来観測量推定部は、推定した前記状態と、前記基底関数のパラメータ、前記パラメータK、展開係数とを用いて、(1)再構成した基底関数の値の推定、(2)基底関数の予測、(3)観測量の予測を行い、
    前記状態は非線形で時間推移するものであり、
    前記観測量は状態の時間推移に伴い非線形で時間推移する観測処理を行うことで得ることができるものである、
    推定装置。
  2. 請求項１の推定装置であって、
    前記エンコーダのパラメータと、前記デコーダのパラメータと、前記パラメータKとは、前記状態推定部をエンコーダとし、前記観測量推定部をデコータとする変分オートエンコーダにより最適化されたものであり、
    前記観測量は時系列観測量である、
    推定装置。
  3. 請求項１から請求項２の何れかの推定装置であって、
    前記観測量推定部は、
    前記状態から中間値を推定する中間値推定部と、
    推定した前記中間値から、観測値を推定する中間観測値推定部と、を含み、
    前記未来観測量推定部は、前記パラメータKを用いて得られる、前記中間値が時間推移により変動した値である未来中間値から前記未来観測量を推定する、
    推定装置。
  4. 請求項１の推定装置で用いるパラメータを学習する学習装置であって、
    学習用の観測量の系列データと、前記基底関数のパラメータ、前記逆関数のパラメータ、前記パラメータK、展開係数とを用いて、(1)基底関数の値の推定、(2)状態の推定、(3)再構成した基底関数の値の推定、(4)基底関数の予測、(5)学習用の観測量の予測を行う推定部と、
    前記学習用の観測量の系列データと、基底関数の推定値の系列データ、状態の推定値の系列データ、再構成した基底関数の推定値の系列データ、基底関数の予測値の系列データ、観測量の予測値の系列データを用いて、(i)観測量の予測誤差、(ii)基底関数の予測誤差、(iii)前記基底関数のパラメータ、前記逆関数のパラメータに基づくニューラルネットワークの重みに対する正則化項、(iv)前記状態の滑らかさと非ガウス性を求め、求めた値から目的関数の値を求める目的関数計算部と、
    前記目的関数に基づき、前記基底関数のパラメータ、前記逆関数のパラメータ、前記パラメータK、前記展開係数とを更新する更新部とを含む、
    学習装置。
  5. エンコーダを用いて、観測量から状態を推定する状態推定ステップと、
    デコーダとを用いて、前記状態推定ステップで推定した前記状態から観測量を推定する観測量推定ステップと、
    時間発展を表すパラメータKを用いて、前記観測量が時間推移により変動した値である未来観測量を推定する未来観測量推定ステップと、を有し、
    前記エンコーダのパラメータと、前記デコーダのパラメータと、前記パラメータKとは、同時に最適化されたものであり、
    前記状態推定ステップにおいて行う処理は第一の関数で、前記観測量推定ステップにおいて行う処理は第二の関数で定義され、
    前記第一の関数は前記第二の関数の逆関数であり、
    前記第二の関数はオートエンコーダのデコーダのパラメータである、基底関数のパラメータを用い、前記第一の関数は前記オートエンコーダのエンコーダのパラメータである、前記基底関数の逆関数のパラメータを用い、
    前記未来観測量推定ステップは、推定した前記状態と、前記基底関数のパラメータ、前記パラメータK、展開係数とを用いて、(1)再構成した基底関数の値の推定、(2)基底関数の予測、(3)観測量の予測を行い、
    前記状態は非線形で時間推移するものであり、
    前記観測量は状態の時間推移に伴い非線形で時間推移する観測処理を行うことで得ることができるものである、
    推定方法。
  6. 請求項１から請求項３の何れかの推定装置、または、請求項４の学習装置としてコンピュータを機能させるためのプログラム。
JP2020569354A 2019-01-28 2019-08-30 推定装置、学習装置、それらの方法、およびプログラム Active JP7163977B2 (ja)

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
JP2019011981 2019-01-28
JP2019011981 2019-01-28
PCT/JP2019/034216 WO2020158032A1 (ja) 2019-01-28 2019-08-30 推定装置、学習装置、それらの方法、およびプログラム

Publications (2)

Publication Number Publication Date
JPWO2020158032A1 JPWO2020158032A1 (ja) 2021-09-30
JP7163977B2 true JP7163977B2 (ja) 2022-11-01

Family

ID=71840041

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2020569354A Active JP7163977B2 (ja) 2019-01-28 2019-08-30 推定装置、学習装置、それらの方法、およびプログラム

Country Status (3)

Country Link
US (1) US20240028872A1 (ja)
JP (1) JP7163977B2 (ja)
WO (1) WO2020158032A1 (ja)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP7542821B2 (ja) * 2020-10-27 2024-09-02 日本電信電話株式会社 学習方法、系列解析方法、学習装置、系列解析装置、及びプログラム

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2018052875A1 (en) 2016-09-15 2018-03-22 Google Llc Image depth prediction neural networks
JP2019502212A (ja) 2016-01-14 2019-01-24 株式会社Ｐｒｅｆｅｒｒｅｄ Ｎｅｔｗｏｒｋｓ 時系列データ適合およびセンサ融合のシステム、方法、および装置

Family Cites Families (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8095210B2 (en) * 2007-01-19 2012-01-10 California Institute Of Technology Prosthetic devices and methods and systems related thereto
US10643153B2 (en) * 2017-04-24 2020-05-05 Virginia Tech Intellectual Properties, Inc. Radio signal identification, identification system learning, and identifier deployment
US11334807B1 (en) * 2017-06-23 2022-05-17 Virginia Tech Intellectual Properties, Inc. Learning approximate estimation networks for communication channel state information
US20190228110A1 (en) * 2018-01-19 2019-07-25 General Electric Company System and method for abstracting characteristics of cyber-physical systems
US12033080B2 (en) * 2018-06-15 2024-07-09 Google Llc Sparse recovery autoencoder
US20200082272A1 (en) * 2018-09-11 2020-03-12 International Business Machines Corporation Enhancing Data Privacy in Remote Deep Learning Services
US10671889B2 (en) * 2018-09-27 2020-06-02 Deepmind Technologies Limited Committed information rate variational autoencoders
EP3667439A1 (en) * 2018-12-13 2020-06-17 ABB Schweiz AG Predictions for a process in an industrial plant
US11172195B2 (en) * 2019-01-26 2021-11-09 Qualcomm Incorporated Excluding intra coded reference samples from local illumination compensation parameter derivation

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2019502212A (ja) 2016-01-14 2019-01-24 株式会社Ｐｒｅｆｅｒｒｅｄ Ｎｅｔｗｏｒｋｓ 時系列データ適合およびセンサ融合のシステム、方法、および装置
WO2018052875A1 (en) 2016-09-15 2018-03-22 Google Llc Image depth prediction neural networks

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
石黒 慎 他，リアルタイム人口を用いたStacked denoising Autoencodersによるタクシー将来需要予測，情報処理学会 論文誌（ジャーナル），情報処理学会，2019年01月15日，第60巻 第1号，第118頁-第128頁

Also Published As

Publication number Publication date
US20240028872A1 (en) 2024-01-25
WO2020158032A1 (ja) 2020-08-06
JPWO2020158032A1 (ja) 2021-09-30

Similar Documents

Publication Publication Date Title
US11797872B2 (en) Quantum bit prediction
US11574164B2 (en) Neural network cooperation
Zhang et al. Remaining useful life estimation for mechanical systems based on similarity of phase space trajectory
US20230316720A1 (en) Anomaly detection apparatus, anomaly detection method, and program
CN108140146A (zh) 用于使用绝热量子计算机进行机器学习的离散变分自动编码器系统和方法
WO2019220620A1 (ja) 異常検出装置、異常検出方法及びプログラム
JP7542821B2 (ja) 学習方法、系列解析方法、学習装置、系列解析装置、及びプログラム
US20220230067A1 (en) Learning device, learning method, and learning program
WO2022009010A1 (en) Model fidelity monitoring and regeneration for manufacturing process decision support
KR20220145007A (ko) 텐서 스트림 데이터에서 손실, 이상치, 패턴을 검출 및 복구하는 데이터 처리 방법
CN117708728A (zh) 传感器测量异常检测
CN120124767A (zh) 一种量子神经网络训练方法、数据处理方法、设备及介质
EP4402608B1 (en) System and method for training of neural network model for control of high dimensional physical systems
JP7163977B2 (ja) 推定装置、学習装置、それらの方法、およびプログラム
KR20190141581A (ko) 데이터 예측을 위한 인공신경망을 학습하는 방법 및 장치
JP7655446B2 (ja) 解析装置及びプログラム
Li et al. Structural damage quantification using long short-term memory (LSTM) auto-encoder and impulse response functions
JP2019095894A (ja) 推定装置、学習装置、学習済みモデル、推定方法、学習方法、及びプログラム
JP7579632B2 (ja) 推定装置、システム及び方法
Li et al. Pygpphs: A python package for Bayesian modeling of port-Hamiltonian systems
JP7672844B2 (ja) 性能予測装置、性能予測方法及び性能予測プログラム
CN110309869A (zh) 面向未知场景的稳定识别方法及装置
JP7575008B2 (ja) 制御方法、制御装置、及びプログラム
JP2019113962A (ja) 解析装置、解析方法及びプログラム
Kocijan System identification with GP models

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20210510

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20220322

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20220519

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20220920

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20221003

R150 Certificate of patent or registration of utility model

Ref document number: 7163977

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150

S533 Written request for registration of change of name

Free format text: JAPANESE INTERMEDIATE CODE: R313533

R350 Written notification of registration of transfer

Free format text: JAPANESE INTERMEDIATE CODE: R350