JP6397380B2 - 時空間変数予測装置及びプログラム - Google Patents

Description

本発明は、時空間変数予測装置及びプログラムに関する。

従来の技術として、時空間変数の時系列データから、クリギングモデルを用いて未観測の時空間変数を推定あるいは予測する手法がある。いま、地理的・時間的広がりを持った空間内のＮ箇所｛ｘ_１，…，ｘ_Ｎ｝でＮ個の観測値｛ｔ_１，…，ｔ_Ｎ｝が得られたとする。クリギングモデルでは、新しい観測地点ｘ_＊における予測値ｔ_＊を確率分布として出力する。

クリギングモデルは、入力変数（位置、時間など）が似ているほど近い値を持つという仮定に基づいたモデルである。「似ている」ということの定義は問題設定に応じて任意に決めることができ、セミバリオグラム、コバリオグラムと呼ばれる関数で定義される。各種のモデル化に使用される代表的なセミバリオグラムの例として、ナゲット効果モデル、球形モデル、指数モデル、ガウスモデルなどが挙げられる。予測の前処理として、全データを用いてセミバリオグラム、コバリオグラムの選択とパラメータ推定を行う。モデルの選択とパラメータの推定には、全データの時空間変数のペアの差分について線形和をとった“経験バリオグラム”を用いる。

上記の定義に従ってデータ点から設計した経験バリオグラムに対しセミバリオグラムモデルの当てはめを行い、加重最小二乗法を用いてモデル選択とパラメータ推定を行う。予測に使用されるクリギングモデルとして単純クリギング、普遍クリギング、ブロッククリギングなどが提案されている。予測に用いるモデルは問題設定に応じて任意に決めることができる。これらのモデルは全て“時空間変数の確率場は本質的定常性あるいは本質的定常性に準ずる定常性を持つ"という仮定に基づいている。

本質的定常性とは、（１）時空間変数の平均値が入力変数（時間帯や地域）によらず一定で、かつ（２）任意の時空間変数のペアの分散がそれらの入力変数の類似度にのみ依存するという性質である。単純クリギングは時空間変数の確率場に本質的定常性を仮定したモデルであり、新しい入力変数｛ｘ_＊｝における時空間変数の予測値ｔ_＊は次式の平均、分散を持つガウス分布に従う。

ここでγはコバリオグラム、Ｃは要素γ（ｘ_ｎ，ｘ_ｍ）＋β^−１δ_ｎｍを持つＮ×Ｎの共分散行列である。ただしδ_ｎｍはクロネッカーのデルタ、βは定数である。単純クリギングは、ガウス過程による回帰の一種とみなすことができる。ガウス過程では、コバリオグラムの代わりにカーネル関数が用いられる。定常性条件（２）は、「似ている」ということの定義が一意に決まるという仮定に対応している。

前述の通り、クリギングモデル（あるいはガウス過程による回帰）は、“「似ている」ということの定義が一意に決まる”という仮定に基づいて時空間変数の値を予測するものである。しかし、実際には「似ている」入力変数から時空間変数が受ける影響は、地域、時間帯など入力変数の特性によって異なるはずである。そこで、クリギングモデルの定常性の過程を緩め、局所非定常性を導入した混合ガウス過程が提案された。混合ガウス過程は、単一のガウス過程を複数個混合したモデルであり、複数のガウス過程の重ね合わせで表現される。混合ガウス過程では、新しい入力変数ｘ_＊における予測値ｔ_＊は次式の平均、分散を持つガウス分布に従う。

ｐ（ｚ_＊＝ｒ｜ｘ_＊）はｔ_＊がｒ個目の要素から生じた確率であり、負担率と呼ばれる。各入力変数の特性（地域や時間帯など）に応じてデータを複数のクラスタに分類し、クラスタごとにデータを最もよく説明するセミバリオグラムの選択、パラメータ推定を行う。各データ点が属するクラスタと各クラスタのパラメータはＥＭアルゴリズムを用いて同時に推定される（非特許文献１，非特許文献２参照）。

S. De Iaco, D.E. Myers, and D. Posa."Space-time analysis using a general product-sum model.", Statistics & Probability Letters, 52(1):p.21−28, 2001. Benedikt Gr¨aler, Lydia E. Gerharz, and Edzer J. Pebesma. "Spatio-temporal analysis and interpolation of PM10 measurements in Europe.", Technical report, ETC/ACM, 2012.

前述の通り、従来技術の混合ガウス過程は、入力変数の次元・個々の入力変数の性質の違いを考慮せず、入力データをクラスタに分類するものである。しかし、実際には時空間変数のクラスタは時間的・空間的な相関を持つはずである。従来技術では、時間的・空間的な相関を持つ現実世界の時空間変数分布を正確に予測することができないという問題が存在した。

本発明は、上記の点に鑑みてなされたものであり、時間的及び空間的相関を持つ時空間変数の値を精度よく予測することができる時空間変数予測装置及びプログラムを提供することを目的とする。

上記目的を達成するために、本発明に係る時空間変数予測装置は、位置情報及び時間情報を有する入力変数に対する時空間変数の観測値を有する観測データの集合に基づいて、未観測の位置情報及び時間情報に対する時空間変数の値を予測する時空間変数予測装置であって、前記観測データの集合に基づいて、複数のガウス過程を、空間的な広がり及び時間的な広がりに対応する複数の階層で混合した階層混合ガウス過程でモデル化された、前記入力変数に対する時空間変数の値を予測するためのモデルに含まれる、前記複数のガウス過程の各々についての、前記観測データ同士の類似性を定義する関数であるカーネル関数の各々のハイパーパラメータと、前記観測データの各々に対する、前記複数のガウス過程の各々の寄与度を表すパラメータである負担率とを学習する学習部を含んで構成されている。

また、本発明に係る時空間変数予測装置において、前記学習部は、前記観測データの集合と、前記複数のガウス過程のカーネル関数の各々のハイパーパラメータとに基づいて、前記観測データの各々に対する、複数のガウス過程からなる複数のユニットの各々の寄与度を表すパラメータであるユニット負担率、及び前記観測データの各々に対する、前記複数のガウス過程の各々の寄与度を表すパラメータである負担率を推定する負担率推定部と、前記観測データの集合と、前記負担率推定部によって推定された、前記観測データの各々に対する、前記複数のユニットの各々のユニット負担率、及び前記観測データの各々に対する、前記複数のガウス過程の各々の負担率とに基づいて、前記複数のガウス過程の各々に対し、前記ガウス過程のカーネル関数の各々のハイパーパラメータを推定するガウス過程パラメータ推定部と、予め定められた反復終了条件を満たすまで、前記負担率推定部による推定、及び前記ガウス過程パラメータ推定部による推定を繰り返す反復判定部とを含むようにすることができる。

また、本発明に係る時空間変数予測装置は、入力された未観測の位置情報及び時間情報を有する前記入力変数に基づいて、前記入力変数に対する前記複数のガウス過程の各々の寄与度を表すパラメータである負担率を推定し、前記学習部によって学習された、前記複数のガウス過程の各々についての前記カーネル関数の各々のハイパーパラメータと、推定された前記入力変数に対する前記複数のガウス過程の各々の負担率とに基づいて、前記入力変数に対する時空間変数の値を予測する時空間変数算出部を更に含むようにすることができる。

また、本発明のプログラムは、コンピュータを、上記の時空間変数予測装置を構成する各部として機能させるためのプログラムである。

以上説明したように、本発明の時空間変数予測装置及びプログラムによれば、位置情報及び時間情報を有する入力変数に対する時空間変数の観測値を有する観測データの集合に基づいて、複数のガウス過程を、空間的な広がり及び時間的な広がりに対応する複数の階層で混合した階層混合ガウス過程でモデル化された、入力変数に対する時空間変数の値を予測するためのモデルに含まれる、複数のガウス過程の各々についての、観測データ同士の類似性を定義する関数であるカーネル関数の各々のハイパーパラメータと、観測データの各々に対する、複数のガウス過程の各々の寄与度を表すパラメータである負担率とを学習することにより、時間的及び空間的相関を持つ時空間変数の値を精度よく予測することができる、という効果が得られる。

本発明の実施の形態における時空間変数予測装置のブロック図である。 観測データの集合の一例を示す図である。 入力部と出力部との構成例を示す図である。 本発明の実施の形態における時空間変数予測装置の学習処理ルーチンを示すフローチャートである。 本発明の実施の形態における時空間変数予測装置の時空間変数算出処理ルーチンを示すフローチャートである。

以下、図面を参照して本発明の実施の形態を詳細に説明する。

＜概要＞

本発明の実施の形態では、混合ガウス過程を拡張した階層混合ガウス過程を提案する。階層混合ガウス過程は、ガウス過程を階層的に混合したモデルである。階層混合ガウス過程では、新しい入力変数ｘ_＊における予測値ｔ_＊は次式の平均を持つガウス分布に従う。

ｐ（ｚ_＊’＝ｒ’｜ｚ_＊＝ｒ，ｘ_＊）はｔ_＊がｒ’個目のユニットのｒ番目の要素から生じた確率であり、負担率と呼ばれる。ｚ，ｚ’は入力変数の各クラスタに対応する潜在変数である。各データ点が属するクラスタと各クラスタのパラメータはＥＭアルゴリズムを用いて同時に推定される。

本発明の実施の形態に係る時空間変数予測装置は、時空間変数の時系列データ（人口分布、人流・交通流の速度・向き、金やダイヤモンドなど鉱物資源の埋蔵量、降水量などの気象データ、土地価格など）を対象としたものであり、観測データに応じて柔軟に適用できるものである。以下では、実施の形態として、人口密度分布の時系列データが観測データとして与えられた条件の下で、未観測地点あるいは未来の時空間変数分布を推定・予測するという場合について説明する。

＜本発明の実施の形態に係る時空間変数予測装置の構成＞

次に、本発明の実施の形態に係る時空間変数予測装置の構成について説明する。図１に示すように、本発明の実施の形態に係る時空間変数予測装置１００は、ＣＰＵと、ＲＡＭと、後述する各処理ルーチンを実行するためのプログラムや各種データを記憶したＲＯＭと、を含むコンピュータで構成することが出来る。時空間変数予測装置１００は、位置情報及び時間情報を有する入力変数に対する時空間変数の観測値を有する観測データの集合に基づいて、未観測の位置情報及び時間情報での時空間変数の値を予測する。この時空間変数予測装置１００は、機能的には図１に示すように、操作部１０と、人口密度情報記憶部１２と、入力部１３と、演算部１４と、時空間変数算出部２６と、出力部２８とを備えている。操作部１０及び演算部１４は、人口密度情報記憶部１２と接続されている。

操作部１０は、後述する人口密度情報記憶部１２に格納されているデータに対するユーザからの各種操作を受け付ける。各種操作とは、人口密度情報記憶部１２に格納された情報を登録、修正、削除する操作等である。操作部１０の入力手段は、キーボードやマウス、メニュー画面、タッチパネルによるもの等、何でもよい。操作部１０は、マウス等の入力手段のデバイスドライバや、メニュー画面の制御ソフトウェアで実現され得る。

人口密度情報記憶部１２には、観測データの集合が格納されている。

人口密度情報記憶部１２には、後述する演算部１４が解析する観測データの集合が格納されており、演算部１４からの要求に従って、観測データの集合を読み出し、当該情報を演算部１４に送信する。

人口密度情報記憶部１２に格納される観測データの集合は、入力変数ｘと時空間変数ｔとの組み合わせの集合である。本発明の実施の形態では、入力変数ｘは位置及び時刻を有し、時空間変数ｔは人口密度である。ある位置、ある時刻における人口密度は｛ｘ_ｉ，ｔ_ｉ｝と表すことができる。人口密度情報記憶部１２はＷｅｂページを保持するＷｅｂサーバや、データベースを具備するデータベースサーバ等である。

図２に、観測データの集合の一例を示す。図２に示すように、例えば、位置ＩＤ及び時刻を表す情報が入力変数ｘとして格納され、人口密度を表す情報が時空間変数ｔとして格納される。

入力部１３は、未観測の位置情報及び時間情報を有する入力変数を受け付ける。

本発明の実施の形態では、入力部１３で受け付けられた位置情報及び時間情報に対して、後述する演算部１４によって得られた各パラメータに基づいて、時空間変数である人口密度の予測が行われる。

入力部１３の入力手段は、キーボードやマウス、メニュー画面、タッチパネルによるもの等、何でもよい。入力部１３は、マウス等の入力手段のデバイスドライバや、メニュー画面の制御ソフトウェアで実現され得る。

図３に、本実施の形態における入力部１３の構成例を示す。図３の構成例では、入力部１３と後述する出力部２８とが１つの画面として構成されている場合を示す。

図３に示すように、入力部１３は、予測を行う対象となる地点あるいは領域に関する位置情報と、予測を行う時間情報とを含む入力変数を受け付ける。また、出力部２８は、後述する時空間変数算出部２６によって出力された、入力変数に対する時空間変数の値を予測結果として表示する。

演算部１４は、学習部１６と、負担率パラメータ格納部２２と、ガウス過程パラメータ格納部２４とを備える。演算部１４では、学習データとして、ガウス過程のハイパーパラメータ、各ガウス過程からの負担率を表すパラメータを算出する。

ここで、本発明の実施の形態における原理について説明する。

人口密度情報記憶部１２に、時空間変数データである観測データの集合Ｄ＝｛ｘ_ｉ，ｔ_ｉ｝^Ｎ _ｉ＝１が格納されたとする。ここでｘ_ｉは位置及び時刻など複数の変数を含むベクトル、ｔ_ｉは人口密度である。

ここで解くべき問題は、入力部１３によって受け付けた、（１）新しい入力変数（未観測の位置及び時間）ｘ_＊における人口密度の値ｔ_＊の予測と、（２）未観測地点を含むある領域における人口密度の値の予測である。

本実施の形態では、問題（１）に焦点を当てて説明する。人口密度分布の予測面を求めるには、ｘ_＊を動かしながら問題（１）を繰り返し解けばよい。

時空間変数の予測値ｔ_＊が従うモデルが、Ｒ個のガウス過程の混合モデルで記述できると仮定する。この仮定に基づけば、新たな入力変数ｘ_＊に対応する値ｔ_＊の期待値を、以下の式（１）のように書き下すことができる。

ここでＸは入力変数ｘの集合（Ｘ＝｛ｘ_ｉ｝^Ｎ _ｉ＝１）、Θ＝｛θ_ｒｒ’｝^Ｒ，Ｒ’ _{ｒ，ｒ’＝１}はｒ’番目のユニットのｒ個目のガウス過程のハイパーパラメータ、ｚ，ｚ’は入力変数の特性（地域や時間帯など）によって分類されたクラスに対応する潜在変数である。ｘ_＊がｒ番目のガウス過程から生じたと仮定すると、予測値ｔ_＊は次式の平均、分散を持つガウス分布に従う。

ここでｋはカーネル関数、ｋ^ｒｒ’は要素ｋ^ｒｒ’（ｘ_＊，Ｘ）を持つベクトル、Ｃ^ｒｒ’は要素ｋ^ｒｒ’（ｘ_ｎ，ｘ_ｍ）＋β^−１δ_ｎｍ＋Ψ^ｒｒ’ _ｎｍを持つＮ×Ｎの共分散行列である。また、ｔは、要素をｔ_ｉとするベクトルである。ただしΨ^ｒｒ’は次式で定義される対角行列である。

ここでσ_ｎは定数である。ｉ番目のデータ点がｒ番目のガウス過程の寄与を受けないとき、すなわち負担率ｐ（ｚ_ｉ’＝ｒ’，ｚ_ｉ＝ｒ｜ｘ_ｉ）のとき、Ｃ^ｒｒ’は無限大に発散し、（Ｃ^ｒｒ’）^−１の（ｉ，ｉ）成分は０となる。すなわちｔ_ｉはｒ’個目のユニットのｒ番目のガウス過程のパラメータ推定に影響しない。

従来手法では、カーネル関数ｋ（ｘ_ｎ，ｘ_ｍ）は予測の前処理段階で一意に決められる。一方、本発明の実施の形態では、カーネル関数を自由に設計し、各々のカーネル関数のパラメータと重みとを自動推定する。本発明の実施の形態では、気象データや環境データなどの時系列地理統計解析で一般的に使われる代表的なコバリオグラムを二つ挙げる。

ここで入力変数ｘの次元数をｎ_ｘとおくと

である。また、上記のコバリオグラムに対応するカーネル関数は次式のように書き下せる。

本発明の実施の形態では、上記のカーネル関数の線形和で定義した、以下の式で示す新たなカーネル関数を用いる。

ここで推定すべきガウス過程のハイパーパラメータは

である。入力変数の集合Ｘ＝｛ｘ_ｉ｝^Ｎ _ｉ＝１、およびガウス過程のハイパーパラメータの集合Θ＝（θ^１１，・・・，θ^ＲＲ’）が与えられた条件の下で、観測データの集合Ｄの尤度関数は次式のように書き下せる。

ここでｔ_＼ｉはｔからｉ番目の要素を除いたものを表す。

本発明の実施の形態では、上記式の尤度関数を最大化するパラメータΘ，ｚ，ｚ’を推定するため、ＥＭアルゴリズムを用いる。ＥＭアルゴリズムでは、次の手順［１］〜［４］でパラメータ推定を行う。なお、本発明の実施の形態では、負担率の事前分布として以下の式（４）〜（５）に示すＳｏｆｔｍａｘ関数を採用し、後述する時空間変数算出部２６による時空間変数の予測で用いられるＳｏｆｔｍａｘ関数のパラメータを算出する。

［１］各パラメータの初期値を選択する。

［２］ＥＭアルゴリズムのＥステップにおいて、観測データ集合の各観測データに対し、以下の式に従って、当該観測データの入力変数に対する、ｒ’番目のユニットのｒ番目のガウス過程の負担率ｐ（ｚ_ｉ＝ｒ｜ｚ_ｉ’＝ｒ’，Ｘ，ｔ，θ^ｒｒ’）、当該観測データの入力変数に対する、ｒ番目のガウス過程の負担率ｐ（ｚ_ｉ＝ｒ｜Ｘ，ｔ，θ^ｒｒ’）、当該観測データの入力変数に対する、ｒ’番目のユニットのｒ番目のガウス過程の同時確率ｐ（ｚ_ｉ＝ｒ，ｚ_ｉ’＝ｒ’｜Ｘ，ｔ，θ^ｒｒ’）を計算する。

なお、上記式（２）におけるｐ（ｔ_ｉ｜Ｘ，ｔ，θ^ｒｒ’）は、次式に示す平均μ_ｉ ^ｒｒ’、分散Σ_ｉ ^ｒｒ’を持つガウス分布で表される。

［３］ＥＭアルゴリズムのＭステップにおいて、観測データ集合の各観測データに対して計算された、ｒ’番目のユニットのｒ番目のガウス過程の負担率ｐ（ｚ_ｉ＝ｒ｜ｚ_ｉ’＝ｒ’，Ｘ，ｔ，θ^ｒｒ’）、ｒ番目のガウス過程の負担率ｐ（ｚ_ｉ＝ｒ｜Ｘ，ｔ，θ^ｒｒ’）、ｒ’番目のユニットのｒ番目のガウス過程の同時確率ｐ（ｚ_ｉ＝ｒ，ｚ_ｉ’＝ｒ’｜Ｘ，ｔ，θ^ｒｒ’）に基づいて、尤度関数Ｑを最大化するガウス過程のパラメータΘ^ｎｅｗを計算する。ここで尤度関数Ｑを次式で表されるＱ^〜で予め近似する。

ここで、π_ｉ ^ｒｒ’は、

で定義される。

また、目的関数Ｑ^〜のガウス過程のハイパーパラメータθ^ｒｒ’に関する微分は次式で書き下せる。

ここでθ^ｒｒ’ _ｊはｒ番目のガウス過程のｊ番目のハイパーパラメータである。降下勾配法、準ニュートン法等を用いれば目的関数Ｑ^〜を最大化するガウス過程のパラメータΘを得ることができる。更新式は次式のように書き下せる。

また、負担率を表すＳｏｆｔｍａｘ関数のパラメータｖ^ｒ及びｖ^ｒｒ’の更新式は次式のように書き下せる。

ここでηは定数である。

［４］ＥＭアルゴリズムの収束条件が満たされているか調べ、満たされていなければ

を実行し、上記手順［２］に戻る。収束条件が満たされていれば、処理を終了する。

従って、学習部１６は、人口密度情報記憶部１２に格納された観測データの集合に基づいて、観測データ同士の類似性を定義する関数であるカーネル関数の各々のハイパーパラメータと、観測データの各々に対する、複数のガウス過程の各々の寄与度を表すパラメータである負担率とを学習する。

本発明の実施の形態におけるカーネル関数は、複数のガウス過程を、空間的な広がり及び時間的な広がりに対応する複数の階層で混合した階層混合ガウス過程でモデル化される。また、本発明の実施の形態におけるカーネル関数は、入力変数に対する時空間変数の値を予測するためのモデルに含まれる、複数のガウス過程の各々についての、観測データ同士の類似性を定義する関数である。

学習部１６は、負担率推定部１８と、ガウス過程パラメータ推定部２０と、反復判定部２１とを備える。

負担率推定部１８は、人口密度情報記憶部１２に格納された観測データの集合と、複数のガウス過程のカーネル関数の各々のハイパーパラメータの初期値、又はガウス過程パラメータ推定部２０による前回推定値とに基づいて、上記式（６）〜（８）に従って、観測データの各々に対する、複数のガウス過程からなる複数のユニットの各々の寄与度を表すパラメータであるユニット負担率、及び観測データの各々に対する、複数のガウス過程の各々の寄与度を表すパラメータである負担率を推定する。

ガウス過程パラメータ推定部２０は、人口密度情報記憶部１２に格納された観測データの集合と、負担率推定部１８によって推定された、観測データの各々に対する、複数のユニットの各々のユニット負担率、及び観測データの各々に対する、複数のガウス過程の各々の負担率とに基づいて、上記式（９）〜（１０）、（１３）に従って、複数のガウス過程の各々に対し、ガウス過程のカーネル関数の各々のハイパーパラメータを推定する。

また、ガウス過程パラメータ推定部２０は、人口密度情報記憶部１２に格納された観測データの集合に基づいて、上記式（１１）〜（１２）に従って、Ｓｏｆｔｍａｘ関数のパラメータｖ^ｒ及びｖ^ｒｒ’を推定する。

反復判定部２１は、予め定められた反復終了条件を満たすまで、負担率推定部１８による推定、及びガウス過程パラメータ推定部２０による推定を繰り返す。そして、反復判定部２１は、予め定められた反復終了条件が満たされた場合には、ガウス過程パラメータ推定部２０によって推定されたＳｏｆｔｍａｘ関数のパラメータｖ^ｒ及びｖ^ｒｒ’を負担率パラメータ格納部２２に格納し、ガウス過程パラメータ推定部２０によって推定されたハイパーパラメータをガウス過程パラメータ格納部２４に格納する。

負担率パラメータ格納部２２には、ガウス過程パラメータ推定部２０によって推定された、Ｓｏｆｔｍａｘ関数のパラメータｖ^ｒ及びｖ^ｒｒ’が格納される。

ガウス過程パラメータ格納部２４には、ガウス過程パラメータ推定部２０によって推定されたガウス過程のカーネル関数の各々のハイパーパラメータが格納される。

時空間変数算出部２６は、入力部１３によって受け付けた、未観測の位置情報及び時間情報を有する入力変数と、負担率パラメータ格納部２２に格納されたＳｏｆｔｍａｘ関数のパラメータとに基づいて、入力変数に対する複数のガウス過程の各々の寄与度を表すパラメータである負担率を推定する。

具体的には、時空間変数算出部２６は、入力変数ｘ_＊と、負担率パラメータ格納部２２に格納されたＳｏｆｔｍａｘ関数のパラメータｖ^ｒ及びｖ^ｒｒ’とに基づいて、以下の式（１４）〜（１５）に従って、入力変数ｘ_＊に対する複数のガウス過程の各々の寄与度を表すパラメータである負担率を推定する。

そして、時空間変数算出部２６は、推定した入力変数ｘ_＊に対する負担率と、ガウス過程パラメータ格納部２４に格納されたハイパーパラメータに基づいて、上記式（１）〜（３）に従って、入力変数に対する時空間変数の値を予測する。

出力部２８は、時空間変数算出部２６によって予測された、入力された入力変数に対する時空間変数の値を、結果として出力する。

例えば、出力部２８は、上記図３に示すように、予測対象の地点あるいは領域及び時間での時空間変数の予測値である人口密度と、当該人口密度に関連する情報である混雑度とを結果として出力する。

＜本発明の実施の形態に係る時空間変数予測装置の作用＞

次に、本発明の実施の形態に係る時空間変数予測装置１００の作用について説明する。

＜学習処理ルーチン＞
まず、時空間変数予測装置１００は、操作部１０より観測データの集合Ｄが入力されると、観測データの集合Ｄを人口密度情報記憶部１２に格納する。そして、時空間変数予測装置１００は、図４に示す学習処理ルーチンを実行する。

まず、ステップＳ１００では、繰り返し変数ｊと各パラメータとを初期化する。

次に、ステップＳ１０２において、反復判定部２１は、予め定められた反復終了条件を満たしたか否かを判定する。予め定められた条件として、繰り返し変数ｊが予め定められた値Ｎ_ｉｔｅｒ未満である場合には、ステップＳ１１０へ進む。一方、繰り返し変数ｊがＮ_ｉｔｅｒ以上である場合には、ステップＳ１０４へ進む。

ステップＳ１０４において、負担率推定部１８は、人口密度情報記憶部１２に格納された観測データの集合と、上記ステップＳ１００で設定された、複数のガウス過程のカーネル関数の各々のハイパーパラメータの初期値、又は後述するステップＳ１０６で前回推定されたガウス過程のカーネル関数の各々のハイパーパラメータとに基づいて、上記式（６）〜（８）に従って、観測データの各々に対するユニット負担率、及び観測データの各々に対する負担率を推定する。

ステップＳ１０６において、ガウス過程パラメータ推定部２０は、人口密度情報記憶部１２に格納された観測データの集合と、上記ステップＳ１０４で推定された、観測データの各々に対するユニット負担率、及び観測データの各々に対する負担率とに基づいて、上記式（９）〜（１０）、（１３）に従って、複数のガウス過程の各々に対し、ガウス過程のカーネル関数の各々のハイパーパラメータを推定する。また、ガウス過程パラメータ推定部２０は、人口密度情報記憶部１２に格納された観測データの集合に基づいて、上記式（１１）〜（１２）に従って、Ｓｏｆｔｍａｘ関数のパラメータを推定する。

ステップＳ１０８では、繰り返し変数ｊを１インクリメントしてステップＳ１０２へ戻る。

ステップＳ１１０において、反復判定部２１は、上記ステップＳ１０６で推定されたＳｏｆｔｍａｘ関数のパラメータを負担率パラメータ格納部２２に格納し、推定されたガウス過程のハイパーパラメータをガウス過程パラメータ格納部２４に格納して、学習処理ルーチンを終了する。

次に、図５に示す時空間変数算出処理ルーチンについて説明する。

学習処理ルーチンが実行され、負担率パラメータ格納部２２にＳｏｆｔｍａｘ関数のパラメータｖ^ｒｒ’が格納され、ガウス過程パラメータ格納部２４にガウス過程のハイパーパラメータが格納され、予測を行う対象となる未観測の位置情報及び時間情報を含む入力変数が入力されると、時空間変数予測装置１００は、図５に示す時空間変数算出処理ルーチンを実行する。

＜時空間変数算出処理ルーチン＞
ステップＳ２００において、入力部１３は、未観測の位置情報及び時間情報を含む入力変数を受け付ける。

ステップＳ２０２において、時空間変数算出部２６は、負担率パラメータ格納部２２に格納されたＳｏｆｔｍａｘ関数のパラメータと、ガウス過程パラメータ格納部２４に格納されたガウス過程のハイパーパラメータとを読み出す。

ステップＳ２０４において、時空間変数算出部２６は、上記ステップＳ２００で受け付けた入力変数と、上記ステップＳ２０２で読み込まれたＳｏｆｔｍａｘ関数のパラメータとに基づいて、上記式（１４）〜（１５）に従って、入力変数に対する複数のガウス過程の各々の寄与度を表すパラメータである負担率を推定する。そして、時空間変数算出部２６は、推定された入力変数に対する負担率と、上記ステップＳ２０２で読み込まれたガウス過程のハイパーパラメータとに基づいて、上記式（１）〜（３）に従って、入力変数に対する時空間変数の値である人口密度を予測する。

出力部２８は、時空間変数算出部２６によって予測された、入力変数に対する時空間変数の値である人口密度を結果として出力して、時空間変数算出処理ルーチンを終了する。

以上説明したように、本発明の実施の形態に係る時空間変数予測装置によれば、位置情報及び時間情報を有する入力変数に対する時空間変数の観測値を有する観測データの集合に基づいて、複数のガウス過程を、空間的な広がり及び時間的な広がりに対応する複数の階層で混合した階層混合ガウス過程でモデル化された、入力変数に対する時空間変数の値を予測するためのモデルに含まれる、複数のガウス過程の各々についての、観測データ同士の類似性を定義する関数であるカーネル関数の各々のハイパーパラメータと、観測データの各々に対する、複数のガウス過程の各々の寄与度を表すパラメータである負担率とを学習することにより、時間的及び空間的相関を持つ時空間変数の値を精度よく予測することができる。

また、時間的及び空間的相関を持つ時空間変数の値を正確に予測することができる。

なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。

例えば、上記の実施の形態では、人口密度分布の時系列データが観測データとして与えられた条件の下で、未観測地点あるいは未来の時空間変数分布を推定・予測するという場合について説明したが、他の様々な時空間変数の時系列データ（人口分布、人流・交通流の速度・向き、金やダイヤモンドなど鉱物資源の埋蔵量、降水量などの気象データ、土地価格など）を対象として本発明を適用することができる。

また、本発明の実施の形態においては、時空間変数予測装置１００によって各パラメータを推定し、推定された各パラメータを用いて、入力変数に対応する時空間変数である人口密度を予測する場合を例に説明したが、これに限定されるものではない。例えば、各パラメータを推定する処理と、推定された各パラメータを用いて時空間変数を予測する処理とを別々の装置として構成してもよい。この場合、各パラメータを推定する装置は、操作部１０と、人口密度情報記憶部１２と、演算部１４とを備え、推定された各パラメータを用いて時空間変数を予測する装置は、入力部１３と、時空間変数算出部２６と、出力部２８とを備える。

また、本発明の実施の形態では、負担率の事前分布としてＳｏｆｔｍａｘ関数を採用し、Ｓｏｆｔｍａｘ関数のパラメータを推定する場合を例に説明したが、これに限定されるものではなく、他の関数を用いてもよい。

また、上述の時空間変数予測装置１００は、内部にコンピュータシステムを有しているが、「コンピュータシステム」は、ＷＷＷシステムを利用している場合であれば、ホームページ提供環境（あるいは表示環境）も含むものとする。

また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能であるし、ネットワークを介して提供することも可能である。

１０ 操作部
１２ 人口密度情報記憶部
１３ 入力部
１４ 演算部
１６ 学習部
１８ 負担率推定部
２０ ガウス過程パラメータ推定部
２１ 反復判定部
２２ 負担率パラメータ格納部
２４ ガウス過程パラメータ格納部
２６ 時空間変数算出部
２８ 出力部
１００ 時空間変数予測装置

Claims (4)

1. 位置情報及び時間情報を有する入力変数に対する時空間変数の観測値を有する観測データの集合に基づいて、未観測の位置情報及び時間情報に対する時空間変数の値を予測する時空間変数予測装置であって、
   前記観測データの集合に基づいて、複数のガウス過程を、空間的な広がり及び時間的な広がりに対応する複数の階層で混合した階層混合ガウス過程でモデル化された、前記入力変数に対する時空間変数の値を予測するためのモデルに含まれる、前記複数のガウス過程の各々についての、前記観測データ同士の類似性を定義する関数であるカーネル関数の各々のハイパーパラメータと、前記観測データの各々に対する、前記複数のガウス過程の各々の寄与度を表すパラメータである負担率とを学習する学習部
   を含む時空間変数予測装置。
2. 前記学習部は、
   前記観測データの集合と、前記複数のガウス過程のカーネル関数の各々のハイパーパラメータとに基づいて、前記観測データの各々に対する、複数のガウス過程からなる複数のユニットの各々の寄与度を表すパラメータであるユニット負担率、及び前記観測データの各々に対する、前記複数のガウス過程の各々の寄与度を表すパラメータである負担率を推定する負担率推定部と、
   前記観測データの集合と、前記負担率推定部によって推定された、前記観測データの各々に対する、前記複数のユニットの各々のユニット負担率、及び前記観測データの各々に対する、前記複数のガウス過程の各々の負担率とに基づいて、前記複数のガウス過程の各々に対し、前記ガウス過程のカーネル関数の各々のハイパーパラメータを推定するガウス過程パラメータ推定部と、
   予め定められた反復終了条件を満たすまで、前記負担率推定部による推定、及び前記ガウス過程パラメータ推定部による推定を繰り返す反復判定部とを含む請求項１記載の時空間変数予測装置。
3. 入力された未観測の位置情報及び時間情報を有する前記入力変数に基づいて、前記入力変数に対する前記複数のガウス過程の各々の寄与度を表すパラメータである負担率を推定し、
   前記学習部によって学習された、前記複数のガウス過程の各々についての前記カーネル関数の各々のハイパーパラメータと、推定された前記入力変数に対する前記複数のガウス過程の各々の負担率とに基づいて、前記入力変数に対する時空間変数の値を予測する時空間変数算出部を更に含む請求項１又は２記載の時空間変数予測装置。
4. コンピュータを、請求項１〜請求項３のいずれか１項に記載の時空間変数予測装置を構成する各部として機能させるためのプログラム。

Images