(57)【要約】【課題】非定常的な時空間データに対しても任意地点の近未来の時空間変数値を高精度に推定可能な予測装置、予測方法及びプログラムを提供することである。【解決手段】本発明の一態様は、逐次観測される時空間データに基づき推定される回帰係数により時空間変数値を予測する予測装置であって、逐次更新される各観測期間において観測された時空間データに対して空間内挿推定を実行することによって、前記観測期間における所定の地点での時空間変数値を推定する空間内挿推定部と、前記観測期間において観測された時空間データと推定された時空間変数値とに対して空間回帰分析を実行することによって、前記観測期間における任意地点での時空間変数値を決定するための回帰関数を推定する空間回帰分析部と、前記推定された回帰関数の回帰係数から前記観測期間の将来時点の回帰係数を推定し、前記推定した将来時点の回帰係数による前記回帰関数に従って将来時点における任意地点での時空間変数値を予測する時空間変数予測部と、を有する予測装置に関する。【選択図】図１
2018022237.tif 000002
【特許請求の範囲】
【請求項１】
逐次観測される時空間データに基づき推定される回帰係数により時空間変数値を予測する予測装置であって、
逐次更新される各観測期間において観測された時空間データに対して空間内挿推定を実行することによって、前記観測期間における所定の地点での時空間変数値を推定する空間内挿推定部と、
前記観測期間において観測された時空間データと推定された時空間変数値とに対して空間回帰分析を実行することによって、前記観測期間における任意地点での時空間変数値を決定するための回帰関数を推定する空間回帰分析部と、
前記推定された回帰関数の回帰係数から前記観測期間の将来時点の回帰係数を推定し、前記推定した将来時点の回帰係数による前記回帰関数に従って将来時点における任意地点での時空間変数値を予測する時空間変数予測部と、
を有する予測装置。
【請求項２】
前記空間内挿推定は、クリギング又はガウス過程回帰に基づく、請求項１記載の予測装置。
【請求項３】
前記空間回帰分析は、前記回帰関数としてカーネル関数を用いる、請求項１又は２記載の予測装置。
【請求項４】
前記空間回帰分析は、ｋ平均法を用いたクラスタリング手法により前記回帰関数の中心位置パラメータに関する目的変数を最小化する、請求項１乃至３何れか一項記載の予測装置。
【請求項５】
時空間変数値を決定するための回帰関数の過去の回帰係数から将来の回帰係数を予測する予測装置であって、
前記将来の回帰係数を表す第１の行列が前記過去の回帰係数を表す第２の行列と未知パラメータを表す第３の行列との積を用いて表現可能であるとき、所定のランクを有する２つの低ランク行列によって前記第３の行列を分解する分解部と、
前記第２の行列と前記第３の行列との積と前記第１の行列との差分を最小化する前記２つの低ランク行列を決定する分解行列決定部と、
前記決定された低ランク行列から前記将来の回帰係数を予測する回帰係数予測部と、
を有する予測装置。
【請求項６】
逐次観測される時空間データに基づき推定される回帰係数により時空間変数値を予測する予測装置により実行される予測方法であって、
逐次更新される各観測期間において観測された時空間データに対して空間内挿推定を実行することによって、前記観測期間における所定の地点での時空間変数値を推定するステップと、
前記観測期間において観測された時空間データと推定された時空間変数値とに対して空間回帰分析を実行することによって、前記観測期間における任意地点での時空間変数値を決定するための回帰関数を推定するステップと、
前記推定された回帰関数の回帰係数から前記観測期間の将来時点の回帰係数を推定し、前記推定した将来時点の回帰係数による前記回帰関数に従って将来時点における任意地点での時空間変数値を予測するステップと、
を有する予測方法。
【請求項７】
前記空間内挿推定は、クリギング又はガウス過程回帰に基づく、請求項６記載の予測方法。
【請求項８】
前記空間回帰分析は、前記回帰関数としてカーネル関数を用いる、請求項６又は７記載の予測方法。
【請求項９】
前記空間回帰分析は、ｋ平均法を用いたクラスタリング手法により前記回帰関数の中心位置パラメータに関する目的変数を最小化する、請求項６乃至８何れか一項記載の予測方法。
【請求項１０】
時空間変数値を決定するための回帰関数の過去の回帰係数から将来の回帰係数を予測する予測装置により実行される予測方法であって、
前記将来の回帰係数を表す第１の行列が前記過去の回帰係数を表す第２の行列と未知パラメータを表す第３の行列との積を用いて表現可能であるとき、所定のランクを有する２つの低ランク行列によって前記第３の行列を分解するステップと、
前記第２の行列と前記第３の行列との積と前記第１の行列との差分を最小化する前記２つの低ランク行列を決定するステップと、
前記決定された低ランク行列から前記将来の回帰係数を予測するステップと、
を有する予測方法。
【請求項１１】
請求項１乃至５何れか一項記載の予測装置の各部としてプロセッサを機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、空間的な広がりを有する統計量（人の流れ、交通の流れ、雨量、大気の化学物質量など）がある時間間隔で観測される時空間データにおいて、現在の時刻から過去のある時刻までに観測された時空間データを用いて、現在の時刻から近未来の時刻までの任意地点における時空間変数値を予測する予測装置に関する。
【背景技術】
【０００２】
現在からある一定期間遡った時間と空間にわたる統計量（人の流れ、交通の流れ、雨量、大気の化学物質量など）に基づき現在からある一定期間先の将来の値を予測する問題は、時空間統計解析のための要素技術として重要である。このような時空間予測問題の既存技術として、クリギングやガウス過程回帰に基づく手法がこれまで多数提案されている（非特許文献１）。
【先行技術文献】
【非特許文献】
【０００３】
【非特許文献１】Noel Cressie and Christopher K. Wikle: Statistics for spatio-temporal data, Wiley, 2011
【発明の概要】
【発明が解決しようとする課題】
【０００４】
しかしながら、非特許文献１に記載される従来手法は、空間内挿法としては有効であるが、時空間データの定常性を仮定しているため、将来予測には適用限界がある。なお、ここでの「定常性」とは、統計量の平均値及び分散が一定であって、且つ空間的相関及び時間的相関が特定の空間及び時刻に依らず一定であるという性質である。従って、イベント会場のようにある時刻において大勢の人々がある場所から別の場所に一斉に移動するような非定常的な時空間データでの近未来予測に上述した従来手法が適用されると、人数の急激な変化にほとんど追随できないという問題があった。
【０００５】
上述した問題点に鑑み、本発明の課題は、非定常的な時空間データに対しても任意地点の近未来の時空間変数値を高精度に推定可能な予測装置、予測方法及びプログラムを提供することである。
【課題を解決するための手段】
【０００６】
上記課題を解決するため、本発明の一態様は、逐次観測される時空間データに基づき推定される回帰係数により時空間変数値を予測する予測装置であって、逐次更新される各観測期間において観測された時空間データに対して空間内挿推定を実行することによって、前記観測期間における所定の地点での時空間変数値を推定する空間内挿推定部と、前記観測期間において観測された時空間データと推定された時空間変数値とに対して空間回帰分析を実行することによって、前記観測期間における任意地点での時空間変数値を決定するための回帰関数を推定する空間回帰分析部と、前記推定された回帰関数の回帰係数から前記観測期間の将来時点の回帰係数を推定し、前記推定した将来時点の回帰係数による前記回帰関数に従って将来時点における任意地点での時空間変数値を予測する時空間変数予測部と、を有する予測装置に関する。
【発明の効果】
【０００７】
本発明によると、非定常的な時空間データに対しても任意地点の近未来の時空間変数値を高精度に推定することができる。
【図面の簡単な説明】
【０００８】
【図１】図１は、本発明の一実施例による予測装置の機能構成を示すブロック図である。
【図２】図２は、本発明の一実施例による時空間変数オンライン予測処理の流れを示す説明図である。
【図３】図３は、本発明の一実施例によるカーネルパラメータの推定を示す説明図である。
【図４】図４は、本発明の一実施例による過去の回帰係数群からの回帰係数の推定を示す説明図である。
【図５】図５は、本発明の一実施例による時空間変数予測に用いる学習データを示す説明図である。
【図６】図６は、本発明の一実施例による時空間変数オンライン予測処理を示すフロー図である。
【図７】図７は、本発明と従来技術との予測精度比較を示す図である。
【図８】図８は、本発明の一実施例による予測装置のハードウェア構成を示すブロック図である。
【発明を実施するための形態】
【０００９】
以下、図面に基づいて本発明の実施の形態を説明する。
【００１０】
以下の実施例では、逐次観測される時空間データに基づき推定される回帰係数により時空間変数値をオンライン予測する予測装置が開示される。
【００１１】
まず、図１〜５を参照して、本発明の一実施例による予測装置の構成を説明する。図１は、本発明の一実施例による予測装置の機能構成を示すブロック図である。
【００１２】
図１に示されるように、予測装置１００は、空間内挿推定部１１０、空間回帰分析部１２０及び時空間変数予測部１３０を有し、時空間変数予測部１３０は、分解部１３１、分解行列決定部１３２及び回帰係数予測部１３３を有する。以下の実施例では、時刻が離散時間で表現されているとし、現時点をｔ０とし、現時点ｔ０から過去のτ期間まで（ｔ０,ｔ０−１，...，ｔ０−τ＋１）の時空間データから将来時点の時空間変数値を予測する予測装置１００を説明する。
【００１３】
空間内挿推定部１１０は、逐次更新される各観測期間において観測された時空間データに対して空間内挿推定を実行することによって、当該観測期間における所定の地点での時空間変数値を推定する。空間的に広がりを有する時空間データの場合、観測対象領域の全ての場所で観測することは実用上困難であり、実際には、図２（ａ）に示されるように、少数の観測地点で観測することになる。そのため、非観測地点における時空間変数値を何らかの方法で推定する必要がある。例えば、空間内挿推定部１１０は、クリギング又はガウス過程回帰に基づき、観測地点で観測された時空間データから非観測地点における時空間変数値を推定してもよい。クリギング及びガウス過程回帰は、当該技術分野において公知の手法であり（非特許文献１）、それらの詳細についてはここでは割愛する。これら公知の手法の何れかを用いることによって、空間内挿推定部１１０は、図２（ｂ）に示されるように、観測地点での時空間データから非観測地点での時空間変数値を推定することが可能である。
【００１４】
具体的には、空間内挿推定部１１０は、ｔ＝ｔ０，ｔ０−１，...，ｔ０−τ＋１の各時点で独立に、所定の非観測地点に対する空間内挿推定を実行する。当該空間内挿推定は、ｔ＝ｔ０，ｔ０−１，...，ｔ０−τ＋１の各時点で独立に実行されるため、並列処理が可能となる。ただし、非観測地点の位置は当該観測期間にわたって共通とされる。観測点数及び非観測点数の総数をｎとすると、当該空間内挿推定によって、ｘ１，...，ｘｎの各地点（通常は、２次元座標により表現される）での時空間変数値が、
【００１５】
【数１】
2018022237.tif 000003
として得られる。
【００１６】
空間回帰分析部１２０は、観測期間において観測された時空間データと推定された時空間変数値とに対して空間回帰分析を実行することによって、当該観測期間における任意地点での時空間変数値を決定するための回帰関数を推定する。すなわち、空間回帰分析部１２０は、空間内挿推定部１１０により内挿補完された時空間データに基づく関数回帰によって空間全体にわたって時空間変数を位置の関数として回帰する。具体的には、図２（ｃ）に示されるように、地点ｘでの時空間変数値ｆ（ｘ）を、中心位置パラメータμ及びカーネル幅パラメータβを有するｍ個のカーネル関数Ｋ（ｘ；θｊ）（ただし、ｊ＝１，...，ｍ）の線形和
【００１７】
【数２】
2018022237.tif 000004
として回帰する。例えば、カーネル関数として、Radial Basis Function(RBF)カーネル
【００１８】
【数３】
2018022237.tif 000005
が用いられてもよい。ここで、θｊ＝（μｊ，βｊ）は未知のパラメータであり、μｊは中心位置パラメータであり、βｊはカーネル幅パラメータである。すなわち、２次元平面の場合、カーネル回帰とは、２次元平面上の値をｍ個のカーネル関数Ｋ（ｘ；θ１），...，Ｋ（ｘ；θｍ）の線形和で近似する線形回帰分析であり、μｊは２次元ベクトルとなる。図２（ｃ）において、各三角錐が１つのカーネル関数に相当し、三角錐の位置及び幅がそれぞれ中心位置及びカーネル幅に相当する。
【００１９】
まず、中心位置パラメータμは、例えば、以下の手順で推定できる。現在時点ｔ０から過去のτ期間までに観測された時空間変数（時空間観測データ）から将来の時空間変数値を推定する場合、図３に示されるように、τ＜τ"なる時間幅の観測データを学習データとする。μｊは後述される重み付きｋ平均法により求められる。
【００２０】
具体的には、空間回帰分析部１２０はまず、
【００２１】
【数４】
2018022237.tif 000006
に従って、時間パラメータｉのそれぞれについて、時空間変数値の時間方向の分散を算出する。ここで、
【００２２】
【数５】
2018022237.tif 000007
は、地点ｘｉでの期間τ"にわたる時空間変数値の平均値
【００２３】
【数６】
2018022237.tif 000008
である。
【００２４】
次に、空間回帰分析部１２０は、目的関数
【００２５】
【数７】
2018022237.tif 000009
を最小化するｍ個の中心位置パラメータμ１，...，μｍを求める。ただし、
【００２６】
【数８】
2018022237.tif 000010
であり、γは正の定数である。式（４）を最小化する中心位置パラメータμ１，...，μｍを求めることは、直感的には、時空間変数値の時間変動の大きな位置に可能な限り多数のカーネルを配置することと等価である。
【００２７】
なお、式（４）の最小化問題は、公知の重み付きｋ平均法と呼ばれるクラスタリング手法と同様の解法で解くことができる。当該手法の詳細については、例えば、J.B. MacQueen: Some methods for classification and analysis of multivariate observations, Proceedings of 5th Berkeley Symposium on Mathematical Statistics and Probability, Berkeley University Press, vol.1, pp.281-297,1967を参照されたい。
【００２８】
具体的には、空間回帰分析部１２０は、中心位置パラメータμ１，...，μｍを適切に初期化し、ボロノイ分割及び重み付き重心計算
【００２９】
【数９】
2018022237.tif 000011
を収束するまで反復してもよい。ここで、Ｃｋはボロノイ分割した際に第ｋ分割クラスタに属するｘｉの集合を表す。σ＝１のとき、これは通常のｋ平均法となるのは明らかであろう。ここで、重みσを導入するのは、カーネル関数の代表位置を時空間変数値の分散の大きな位置に可能な限り多くのカーネル関数を配置し、回帰精度を向上させるためである。なお、ｍ個のカーネル関数の中心位置パラメータμ１，...，μｍは、時点ｔ０から時点ｔ０−τ"＋１までの期間にわたって同一である。
【００３０】
次に、カーネル幅パラメータβは、以下の手順で推定できる。空間回帰分析部１２０は、取得した中心位置パラメータμ１，...，μｍの位置関係から各中心位置でのカーネル幅を決定する。具体的には、空間回帰分析部１２０は、μｊに対して、当該μｊを除く中心位置パラメータμｊにｋ番目に近い中心位置パラメータとのユークリッド距離をβｊとする。すなわち、μｊ'をμｊにｋ番目に近い中心位置パラメータとすると、空間回帰分析部１２０は、
【００３１】
【数１０】
2018022237.tif 000012
を求めることになる。各位置についてβｊを直接最適化するのは高コストであるが、上述した手法では、ｍ個のカーネルのカーネル幅が、１つのハイパーパラメータｋのみに依存する形で位置毎に適応的に求められる。なお、最適なハイパーパラメータｋの求め方は、以降において説明される。
【００３２】
このようにして、空間回帰分析部１２０は、時点ｔ０−τ"＋１からｔ０までにわたる共通のカーネル関数を取得でき、すなわち、ｍ個のカーネルの中心位置パラメータμｊ及びカーネル幅パラメータβｊは、時点ｔ０−τ"＋１からｔ０までの期間にわたって同一である。
【００３３】
次に、カーネル回帰係数ｗは、以下の手順で推定できる。空間回帰分析部１２０は、各時点について独立にカーネル回帰係数ｗを推定する。時点ｔでのカーネル回帰係数をｗｔとしたとき、空間回帰分析部１２０は、時点ｔについて空間内挿推定部１１０により取得された時空間変数値
【００３４】
【数１１】
2018022237.tif 000013
を用いて、
【００３５】
【数１２】
2018022237.tif 000014
により表される正則化項付き自乗誤差の最小化によって、時点ｔでの回帰係数を求める。ただし、
【００３６】
【数１３】
2018022237.tif 000015
であり、Ｔはベクトルの転置を表し、λは正則化パラメータ（正の定数）であり、Ｋはｎ×ｍ行列を表し、その第（ｉ，ｊ）要素はカーネル関数値Ｋ（ｘｉ；θｊ）である。式（７）を最小化するｗｔをｗｔ＊とすると、
【００３７】
【数１４】
2018022237.tif 000016
となる。ここで、Ｉｍはｍ次元単位行列を表す。このようにして、ｔ０からｔ０−τ"＋１までの各時点ｔについて最適なカーネル関数回帰を実現できる。ｔ０からｔ０−τ"＋１までの期間においてカーネル関数は共通であるが、カーネル回帰係数については時点ｔ毎に最適な値が推定される。
【００３８】
上記説明から理解されるように、カーネル回帰を実施するためには、上述したハイパーパラメータｋを適切に設定する必要がある。通常、ハイパーパラメータｋを設定するため、交差検定（Cross-Validation: CV）法、特に、Leave-one-out CV(LOOCV)法が多用される。これらの詳細については、例えば、K. P. Murphy: Machine learning, MIT press, 2012を参照されたい。しかしながら、LOOCV法では、（ｋ，ｍ，λ）の各候補値の組に対して上述したカーネル回帰推定の計算が必要となるため、計算時間がかかり、本発明によるオンライン処理には適していない。
【００３９】
この問題に対して、カーネル回帰の場合、LOOCV法と比較してはるかに効率的な一般化交差検定（Generalized CV: GCV）法が適用可能であることが理論的に示すことができるため、空間回帰分析部１２０は、GCV法によりハイパーパラメータｋを推定する。なお、GCV法の詳細については、例えば、P. Craven and G. Wahba: Smoothing noisy data with spline functions, Numerical Mathematics, Vo.31, pp.377-403, 1979を参照されたい。
【００４０】
具体的には、空間回帰分析部１２０は、（ｋ，ｍ，λ）の候補値の組に対して、
【００４１】
【数１５】
2018022237.tif 000017
を最小化する（ｋ，ｍ，λ）を選択する。ここで、
【００４２】
【数１６】
2018022237.tif 000018
は、ある（ｋ，ｍ，λ）の候補値の組について推定した推定値である。また、Ｈｉｉは、カーネル回帰における平滑化行列
【００４３】
【数１７】
2018022237.tif 000019
の第（ｉ，ｉ）要素を表す。GCV法では、
【００４４】
【数１８】
2018022237.tif 000020
の推定は、各(k, m, λ)に対して１回実行すればよく、LOOCV法と比較して計算量が1/nτとなる。
【００４５】
時空間変数予測部１３０は、推定された回帰関数の回帰係数から観測期間の将来時点の回帰パラメータを推定し、推定した将来時点の回帰係数に基づく回帰関数に従って将来時点における任意地点での時空間変数値を予測する。具体的には、時空間変数予測部１３０は、図４に示されるように、現時点及び現時点から過去のτ期間までの既知の回帰係数
【００４６】
【数１９】
2018022237.tif 000021
を用いて、現時点からτ'期間先の将来時点の回帰係数
【００４７】
【数２０】
2018022237.tif 000022
を予測する。この予測問題に対して、時空間変数予測部１３０は、第ｊ番目のカーネル関数の時点ｔ０からτ'期間先の将来時点の回帰係数を、現時点ｔ０から過去のτ期間までの全カーネル関数の回帰係数の線形和、すなわち、
【００４８】
【数２１】
2018022237.tif 000023
により近似する。ここで、εｊは誤差項であり、
【００４９】
【数２２】
2018022237.tif 000024
は未知の線形重みパラメータである。
【００５０】
さらに、式（１１）は時点ｔ０に依らず任意時点でも成立すると仮定し、τｍ次元の説明変数ベクトル
【００５１】
【数２３】
2018022237.tif 000025
と、それの各要素の回帰係数を時間ｌだけシフトした
【００５２】
【数２４】
2018022237.tif 000026
とが類似している場合、前者に対する目的変数の値
【００５３】
【数２５】
2018022237.tif 000027
と後者に対する目的変数の値
【００５４】
【数２６】
2018022237.tif 000028
とが類似しているという仮定に基づき、時空間変数予測部１３０は、図５に示されるように、現時点から過去のτ"期間までの既知の回帰係数学習データを用いて、未知の線形重みパラメータを学習する。当該未知の線形重みパラメータはｊ＝１，...，ｍ毎に独立であるため、時空間変数予測部１３０は、各パラメータｊ毎に独立に線形重み未知パラメータの学習を並列に実行可能である。
【００５５】
ここで、式（１１）は、τｍ次元の説明変数（ｘ１，...，ｘτｍ）とスカラーｙｊとの関係式
【００５６】
【数２７】
2018022237.tif 000029
と数学的に等価である。表記の簡単化のため、以下では式（１２）を用いて未知パラメータθｊ，ｐの推定法を説明する。
【００５７】
ｎ個の学習データの組を
【００５８】
【数２８】
2018022237.tif 000030
とする。明らかに、
【００５９】
【数２９】
2018022237.tif 000031
は上述した回帰係数に相当し、ｙｊ（ｉ）は予測対象の将来時点の回帰係数に相当する。これらが式（１２）を満たすことから、
【００６０】
【数３０】
2018022237.tif 000032
が成り立つ。ここで、Ｙは第（ｉ，ｊ）要素がｙｊ（ｉ）であるｎ×ｍ行列を表し、Ｘは第（ｉ，ｋ）要素がｘｋ（ｉ）であるｎ×τｍ行列を表し、Θは第（ｋ，ｊ）要素がθｊ，ｋであるτｍ×ｍ行列を表し、Ｅは全ての列ベクトルが同じであって、（ε１・・・εｍ）であるｎ×ｍ行列を表す。
【００６１】
未知パラメータθは、最小自乗法により
【００６２】
【数３１】
2018022237.tif 000033
として解析的に求められるが、学習データ数ｎがｎ＜＜τｍであるとき、ＸＴＸの逆行列が特異となるという問題が生じる。
【００６３】
この問題を回避するため、分解部１３１は、ｍ×ｒ行列Ａ及びτｍ×ｒ行列Ｂを用いて、
【００６４】
【数３２】
2018022237.tif 000034
と行列Θを低ランク（ランクｒ）行列に分解する。すなわち、分解部１３１は、将来の回帰係数を表す行列Ｙが過去の回帰係数を表す行列Ｘと未知パラメータを表す行列Θとの積を用いて表現可能であるとき、所定のランクを有する２つの低ランク行列Ａ，Ｂによって行列Θを分解する。
【００６５】
ここで、
【００６６】
【数３３】
2018022237.tif 000035
のΘに関する最小化問題は、
【００６７】
【数３４】
2018022237.tif 000036
の制約の下、
【００６８】
【数３５】
2018022237.tif 000037
を最小化するＡ，Ｂを求める問題となる。分解行列決定部１３２は、行列の積ＸΘと行列Ｙとの差分を最小化する低ランク行列Ａ，Ｂを決定する。当該最適化問題について、未知パラメータ決定部１３２は、
【００６９】
【数３６】
2018022237.tif 000038
とし、さらに、τｍ×τｍ行列
【００７０】
【数３７】
2018022237.tif 000039
の固有ベクトル行列を
【００７１】
【数３８】
2018022237.tif 000040
とし（ｖｊは第ｊ番目の固有値に対する固有ベクトル（ｍ次元列ベクトル）である）、Ａ，Ｂの最適解を、
【００７２】
【数３９】
2018022237.tif 000041
として求めることができる。
【００７３】
回帰係数予測部１３３は、決定された低ランク行列から将来の回帰係数を予測する。具体的には、回帰係数予測部１３３は、求められた低ランク行列Ａ，Ｂの最適解から式（１１）の未知パラメータを求めることができ、式（１１）を用いて将来時点での回帰係数ｗを取得できる。
【００７４】
このようにして、時空間変数予測部１３０は、取得した将来時点での回帰係数ｗを式（１）に適用し、適用後の式（１）を用いて任意地点での将来時点における時空間変数値を予測することができる。
【００７５】
なお、分解部１３１、分解行列決定部１３２及び回帰係数予測部１３３は、必ずしも予測装置１００に内蔵される必要はなく、時空間変数値を決定するための回帰関数の過去の回帰係数から将来の回帰係数を予測する独立した予測装置として構成されてもよい。
【００７６】
上述した空間内挿推定部１１０、空間回帰分析部１２０及び時空間変数予測部１３０による処理は、現時点から過去の一定の期間までの時間ウィンドウに対して実行された。時間の経過に従って当該時間ウィンドウもまた移動し、例えば、時点ｔ０からｔ０−τ＋１までの時間ウィンドウは、次に時点ｔ０＋１からｔ０−τ＋２までの時間ウィンドウに移動するなどである。これら各時間ウィンドウに対して空間内挿推定部１１０、空間回帰分析部１２０及び時空間変数予測部１３０による上記処理が繰り返され、直近の時間ウィンドウに対して新たな将来時点での回帰係数ｗが逐次求められる。このようにして、直近の学習データに基づき算出された回帰係数ｗに基づき、直近の学習データを反映した時空間変数値をオンライン予測することができる。
【００７７】
次に、図６を参照して、本発明の一実施例によるオンライン予測処理を説明する。図６は、本発明の一実施例による時空間変数オンライン予測処理を示すフロー図である。
【００７８】
図６に示されるように、ステップＳ１０１において、予測装置１００は、観測対象エリアの観測地点から観測された時空間データを取得する。例えば、予測装置１００は、通信ネットワークを介し各観測地点に設置された計測手段から観測データを受信する。
【００７９】
ステップＳ１０２において、予測装置１００は、受信した観測データを時空間データとして記憶する。具体的には、予測装置１００は、受信した観測データを送信元の計測手段の位置情報及び観測時間と関連付けて記憶する。
【００８０】
ステップＳ１０３において、予測装置１００は、記憶した時空間データに対して空間内挿推定を実行する。具体的には、予測装置１００は、ある観測期間内の各時点において観測された時空間データに対して空間内挿推定を実行する。例えば、当該空間内挿推定は、クリギング又はガウス過程回帰などの何れか適切な公知技術に基づくものであってもよい。
【００８１】
ステップＳ１０４において、予測装置１００は、所定の非観測地点について内挿補完された時空間データを取得する。典型的には、当該所定の地点は当該観測期間の各時点にわたって固定的に設定された共通の位置である。
【００８２】
ステップＳ１０５において、予測装置１００は、観測地点及び非観測地点の時空間データに対して空間回帰分析を実行し、当該観測期間における任意地点での時空間変数値を決定するための回帰関数を推定する。例えば、当該空間回帰分析は、回帰関数としてカーネル関数を用いてもよいし、あるいは、ｋ平均法を用いたクラスタリング手法により回帰関数（カーネル関数）の中心位置パラメータに関する目的変数を最小化することであってもよい。
【００８３】
ステップＳ１０６において、予測関数１００は、推定された回帰係数を保持する。上述した実施例では、回帰係数はカーネル関数のカーネル回帰係数ｗである。
【００８４】
ステップＳ１０７において、予測関数１００は、推定された過去の回帰係数に対して時空間回帰分析を実行する。具体的には、予測関数１００は、式（１３）に関して上述したように、将来の回帰係数を表す行列Ｙが過去の回帰係数を表す行列Ｘと未知パラメータを表す行列Θとの積を用いて表現可能であるとき、所定のランクを有する２つの分解行列Ａ，Ｂによって行列Θを分解する。その後、予測関数１００は、
【００８５】
【数４０】
2018022237.tif 000042
を最小化する分解行列Ａ，Ｂを決定し、上述した最適化問題の式（１８）及び（１９）により表される最適解から、式（１１）を用いて将来の回帰係数を予測する。
【００８６】
ステップＳ１０８において、予測装置１００は、予測した将来の回帰係数を保持すると共に、これを適用した式（１）に従って将来時点における任意地点の時空間変数値を予測する。
【００８７】
ステップＳ１０９において、予測装置１００は、次の時間ウィンドウに対してステップＳ１０１〜Ｓ１０８の処理を実行するか判断し、当該処理を継続する場合、すなわち、次の観測期間に基づき回帰係数を更新する場合、ステップＳ１０１に戻る。他方、当該処理を停止する場合、予測装置１００は当該処理を終了する。これにより、予測装置１００は、直近の時間ウィンドウに対して新たな将来時点での回帰係数ｗを逐次求めることができ、直近の学習データに基づき算出された回帰係数ｗに基づき時空間変数値をオンライン予測することができる。
【００８８】
図７は、本発明と従来技術との予測精度比較を示す図である。図示されるデータは、あるイベント会場での来場者のカウントデータを用いた実験結果を示す。当該データは、会場に設置された２４カ所のＷｉＦｉアクセスポイントにおいて、１分毎に観測された合計で４００分からなるデータである。上述した実施例において、τ＝１０、τ"＝５０とし、τ'＝５，１０，２０，３０に対する２４カ所での予測精度を評価した。実験評価では、誤差の評価尺度して、平均相対誤差
【００８９】
【数４１】
2018022237.tif 000043
を用いた。ここで、
【００９０】
【数４２】
2018022237.tif 000044
は現在時点ｔに対する将来時点ｔ＋τ'での時空間変数の予測値を表し、ｙｔ＋τ'は真値を表す。上記のＭＲＰＥを２４カ所について算出し、当該２４カ所での平均値及び標準偏差によって評価した。なお、既存手法としてガウス過程と比較した。図６の実験結果から理解されるように、近未来予測（τ'＝５）では、両者に顕著な差はないが、τ'＝１０以降では、本発明がガウス過程と比較して予測精度及び安定性（予測精度のばらつき）の観点で顕著な優位性を有することが確認された。
【００９１】
次に、図８を参照して、本発明の一実施例による予測装置のハードウェア構成を説明する。図８は、本発明の一実施例による予測装置のハードウェア構成を示すブロック図である。
【００９２】
図８に示されるように、予測装置１００は、典型的には、サーバにより実現されてもよく、例えば、バスを介し相互接続されるドライブ装置１０１、補助記憶装置１０２、メモリ装置１０３、プロセッサ１０４、インタフェース装置１０５及び通信装置１０６から構成される。予測装置１００における上述した各種機能及び処理を実現するプログラムを含む各種コンピュータプログラムは、ＣＤ−ＲＯＭ（Ｃｏｍｐａｃｔ  Ｄｉｓｋ−Ｒｅａｄ  Ｏｎｌｙ  Ｍｅｍｏｒｙ）、ＤＶＤ（Ｄｉｇｉｔａｌ  Ｖｅｒｓａｔｉｌｅ  Ｄｉｓｋ）、フラッシュメモリなどの記録媒体１０７によって提供されてもよい。プログラムを記憶した記録媒体１０７がドライブ装置１０１にセットされると、プログラムが記録媒体１０７からドライブ装置１０１を介して補助記憶装置１０２にインストールされる。但し、プログラムのインストールは必ずしも記録媒体１０７により行う必要はなく、ネットワークなどを介し何れかの外部装置からダウンロードするようにしてもよい。補助記憶装置１０２は、インストールされたプログラムを格納すると共に、必要なファイルやデータなどを格納する。メモリ装置１０３は、プログラムの起動指示があった場合に、補助記憶装置１０２からプログラムやデータを読み出して格納する。プロセッサ１０４は、メモリ装置１０３に格納されたプログラムやプログラムを実行するのに必要なパラメータなどの各種データに従って、予測装置１００の各種機能及び処理を実行する。インタフェース装置１０５は、ネットワーク又は外部装置に接続するための通信インタフェースとして用いられる。通信装置１０６は、インターネットなどのネットワークと通信するための各種通信処理を実行する。しかしながら、上述したハードウェア構成は単なる一例であり、予測装置１００は、上述したハードウェア構成に限定されるものでなく、他の何れか適切なハードウェア構成により実現されてもよい。
【００９３】
なお、上述した予測装置１００の各部及びステップＳ１０１〜Ｓ１０９は、コンピュータのメモリ装置１０３に記憶されたプログラムをプロセッサ１０４が実行することによって実現されてもよい。
【００９４】
以上、本発明の実施例について詳述したが、本発明は上述した特定の実施形態に限定されるものではなく、特許請求の範囲に記載された本発明の要旨の範囲内において、種々の変形・変更が可能である。
【符号の説明】
【００９５】
１００  予測装置
１１０  空間内挿推定部
１２０  空間回帰分析部
１３０  時空間変数予測部
１３１  分解部
１３２  分解行列決定部
１３３  回帰係数予測部
【図１】
2018022237.tif 000045
【図２】
2018022237.tif 000046
【図３】
2018022237.tif 000047
【図４】
2018022237.tif 000048
【図５】
2018022237.tif 000049
【図６】
2018022237.tif 000050
【図７】
2018022237.tif 000051
【図８】
2018022237.tif 000052