(57)【特許請求の範囲】
【請求項１】
位置情報及び時間情報を有する入力変数に対する時空間変数の観測値を有する観測データの集合に基づいて、未観測の位置情報及び時間情報に対する時空間変数の値を予測する時空間変数予測装置であって、
前記観測データの集合に基づいて、複数のガウス過程を、空間的な広がり及び時間的な広がりに対応する複数の階層で混合した階層混合ガウス過程でモデル化された、前記入力変数に対する時空間変数の値を予測するためのモデルに含まれる、前記複数のガウス過程の各々についての、前記観測データ同士の類似性を定義する関数であるカーネル関数の各々のハイパーパラメータと、前記観測データの各々に対する、前記複数のガウス過程の各々の寄与度を表すパラメータである負担率とを学習する学習部
を含む時空間変数予測装置。
【請求項２】
前記学習部は、
前記観測データの集合と、前記複数のガウス過程のカーネル関数の各々のハイパーパラメータとに基づいて、前記観測データの各々に対する、複数のガウス過程からなる複数のユニットの各々の寄与度を表すパラメータであるユニット負担率、及び前記観測データの各々に対する、前記複数のガウス過程の各々の寄与度を表すパラメータである負担率を推定する負担率推定部と、
前記観測データの集合と、前記負担率推定部によって推定された、前記観測データの各々に対する、前記複数のユニットの各々のユニット負担率、及び前記観測データの各々に対する、前記複数のガウス過程の各々の負担率とに基づいて、前記複数のガウス過程の各々に対し、前記ガウス過程のカーネル関数の各々のハイパーパラメータを推定するガウス過程パラメータ推定部と、
予め定められた反復終了条件を満たすまで、前記負担率推定部による推定、及び前記ガウス過程パラメータ推定部による推定を繰り返す反復判定部とを含む請求項１記載の時空間変数予測装置。
【請求項３】
入力された未観測の位置情報及び時間情報を有する前記入力変数に基づいて、前記入力変数に対する前記複数のガウス過程の各々の寄与度を表すパラメータである負担率を推定し、
前記学習部によって学習された、前記複数のガウス過程の各々についての前記カーネル関数の各々のハイパーパラメータと、推定された前記入力変数に対する前記複数のガウス過程の各々の負担率とに基づいて、前記入力変数に対する時空間変数の値を予測する時空間変数算出部を更に含む請求項１又は２記載の時空間変数予測装置。
【請求項４】
コンピュータを、請求項１〜請求項３のいずれか１項に記載の時空間変数予測装置を構成する各部として機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、時空間変数予測装置及びプログラムに関する。
【背景技術】
【０００２】
従来の技術として、時空間変数の時系列データから、クリギングモデルを用いて未観測の時空間変数を推定あるいは予測する手法がある。いま、地理的・時間的広がりを持った空間内のＮ箇所｛ｘ１，…，ｘＮ｝でＮ個の観測値｛ｔ１，…，ｔＮ｝が得られたとする。クリギングモデルでは、新しい観測地点ｘ＊における予測値ｔ＊を確率分布として出力する。
【０００３】
クリギングモデルは、入力変数（位置、時間など）が似ているほど近い値を持つという仮定に基づいたモデルである。「似ている」ということの定義は問題設定に応じて任意に決めることができ、セミバリオグラム、コバリオグラムと呼ばれる関数で定義される。各種のモデル化に使用される代表的なセミバリオグラムの例として、ナゲット効果モデル、球形モデル、指数モデル、ガウスモデルなどが挙げられる。予測の前処理として、全データを用いてセミバリオグラム、コバリオグラムの選択とパラメータ推定を行う。モデルの選択とパラメータの推定には、全データの時空間変数のペアの差分について線形和をとった“経験バリオグラム”を用いる。
【０００４】
上記の定義に従ってデータ点から設計した経験バリオグラムに対しセミバリオグラムモデルの当てはめを行い、加重最小二乗法を用いてモデル選択とパラメータ推定を行う。予測に使用されるクリギングモデルとして単純クリギング、普遍クリギング、ブロッククリギングなどが提案されている。予測に用いるモデルは問題設定に応じて任意に決めることができる。これらのモデルは全て“時空間変数の確率場は本質的定常性あるいは本質的定常性に準ずる定常性を持つ"という仮定に基づいている。
【０００５】
本質的定常性とは、（１）時空間変数の平均値が入力変数（時間帯や地域）によらず一定で、かつ（２）任意の時空間変数のペアの分散がそれらの入力変数の類似度にのみ依存するという性質である。単純クリギングは時空間変数の確率場に本質的定常性を仮定したモデルであり、新しい入力変数｛ｘ＊｝における時空間変数の予測値ｔ＊は次式の平均、分散を持つガウス分布に従う。
【０００６】
【数１】
0006397380.tif 000002
【０００７】
ここでγはコバリオグラム、Ｃは要素γ（ｘｎ，ｘｍ）＋β−１δｎｍを持つＮ×Ｎの共分散行列である。ただしδｎｍはクロネッカーのデルタ、βは定数である。単純クリギングは、ガウス過程による回帰の一種とみなすことができる。ガウス過程では、コバリオグラムの代わりにカーネル関数が用いられる。定常性条件（２）は、「似ている」ということの定義が一意に決まるという仮定に対応している。
【０００８】
前述の通り、クリギングモデル（あるいはガウス過程による回帰）は、“「似ている」ということの定義が一意に決まる”という仮定に基づいて時空間変数の値を予測するものである。しかし、実際には「似ている」入力変数から時空間変数が受ける影響は、地域、時間帯など入力変数の特性によって異なるはずである。そこで、クリギングモデルの定常性の過程を緩め、局所非定常性を導入した混合ガウス過程が提案された。混合ガウス過程は、単一のガウス過程を複数個混合したモデルであり、複数のガウス過程の重ね合わせで表現される。混合ガウス過程では、新しい入力変数ｘ＊における予測値ｔ＊は次式の平均、分散を持つガウス分布に従う。
【０００９】
【数２】
0006397380.tif 000003
【００１０】
ｐ（ｚ＊＝ｒ｜ｘ＊）はｔ＊がｒ個目の要素から生じた確率であり、負担率と呼ばれる。各入力変数の特性（地域や時間帯など）に応じてデータを複数のクラスタに分類し、クラスタごとにデータを最もよく説明するセミバリオグラムの選択、パラメータ推定を行う。各データ点が属するクラスタと各クラスタのパラメータはＥＭアルゴリズムを用いて同時に推定される（非特許文献１，非特許文献２参照）。
【先行技術文献】
【非特許文献】
【００１１】
【非特許文献１】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.
【発明の概要】
【発明が解決しようとする課題】
【００１２】
前述の通り、従来技術の混合ガウス過程は、入力変数の次元・個々の入力変数の性質の違いを考慮せず、入力データをクラスタに分類するものである。しかし、実際には時空間変数のクラスタは時間的・空間的な相関を持つはずである。従来技術では、時間的・空間的な相関を持つ現実世界の時空間変数分布を正確に予測することができないという問題が存在した。
【００１３】
本発明は、上記の点に鑑みてなされたものであり、時間的及び空間的相関を持つ時空間変数の値を精度よく予測することができる時空間変数予測装置及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【００１４】
上記目的を達成するために、本発明に係る時空間変数予測装置は、位置情報及び時間情報を有する入力変数に対する時空間変数の観測値を有する観測データの集合に基づいて、未観測の位置情報及び時間情報に対する時空間変数の値を予測する時空間変数予測装置であって、前記観測データの集合に基づいて、複数のガウス過程を、空間的な広がり及び時間的な広がりに対応する複数の階層で混合した階層混合ガウス過程でモデル化された、前記入力変数に対する時空間変数の値を予測するためのモデルに含まれる、前記複数のガウス過程の各々についての、前記観測データ同士の類似性を定義する関数であるカーネル関数の各々のハイパーパラメータと、前記観測データの各々に対する、前記複数のガウス過程の各々の寄与度を表すパラメータである負担率とを学習する学習部を含んで構成されている。
【００１５】
また、本発明に係る時空間変数予測装置において、前記学習部は、前記観測データの集合と、前記複数のガウス過程のカーネル関数の各々のハイパーパラメータとに基づいて、前記観測データの各々に対する、複数のガウス過程からなる複数のユニットの各々の寄与度を表すパラメータであるユニット負担率、及び前記観測データの各々に対する、前記複数のガウス過程の各々の寄与度を表すパラメータである負担率を推定する負担率推定部と、前記観測データの集合と、前記負担率推定部によって推定された、前記観測データの各々に対する、前記複数のユニットの各々のユニット負担率、及び前記観測データの各々に対する、前記複数のガウス過程の各々の負担率とに基づいて、前記複数のガウス過程の各々に対し、前記ガウス過程のカーネル関数の各々のハイパーパラメータを推定するガウス過程パラメータ推定部と、予め定められた反復終了条件を満たすまで、前記負担率推定部による推定、及び前記ガウス過程パラメータ推定部による推定を繰り返す反復判定部とを含むようにすることができる。
【００１６】
また、本発明に係る時空間変数予測装置は、入力された未観測の位置情報及び時間情報を有する前記入力変数に基づいて、前記入力変数に対する前記複数のガウス過程の各々の寄与度を表すパラメータである負担率を推定し、前記学習部によって学習された、前記複数のガウス過程の各々についての前記カーネル関数の各々のハイパーパラメータと、推定された前記入力変数に対する前記複数のガウス過程の各々の負担率とに基づいて、前記入力変数に対する時空間変数の値を予測する時空間変数算出部を更に含むようにすることができる。
【００１７】
また、本発明のプログラムは、コンピュータを、上記の時空間変数予測装置を構成する各部として機能させるためのプログラムである。
【発明の効果】
【００１８】
以上説明したように、本発明の時空間変数予測装置及びプログラムによれば、位置情報及び時間情報を有する入力変数に対する時空間変数の観測値を有する観測データの集合に基づいて、複数のガウス過程を、空間的な広がり及び時間的な広がりに対応する複数の階層で混合した階層混合ガウス過程でモデル化された、入力変数に対する時空間変数の値を予測するためのモデルに含まれる、複数のガウス過程の各々についての、観測データ同士の類似性を定義する関数であるカーネル関数の各々のハイパーパラメータと、観測データの各々に対する、複数のガウス過程の各々の寄与度を表すパラメータである負担率とを学習することにより、時間的及び空間的相関を持つ時空間変数の値を精度よく予測することができる、という効果が得られる。
【図面の簡単な説明】
【００１９】
【図１】本発明の実施の形態における時空間変数予測装置のブロック図である。
【図２】観測データの集合の一例を示す図である。
【図３】入力部と出力部との構成例を示す図である。
【図４】本発明の実施の形態における時空間変数予測装置の学習処理ルーチンを示すフローチャートである。
【図５】本発明の実施の形態における時空間変数予測装置の時空間変数算出処理ルーチンを示すフローチャートである。
【発明を実施するための形態】
【００２０】
以下、図面を参照して本発明の実施の形態を詳細に説明する。
【００２１】
＜概要＞
【００２２】
本発明の実施の形態では、混合ガウス過程を拡張した階層混合ガウス過程を提案する。階層混合ガウス過程は、ガウス過程を階層的に混合したモデルである。階層混合ガウス過程では、新しい入力変数ｘ＊における予測値ｔ＊は次式の平均を持つガウス分布に従う。
【００２３】
【数３】
0006397380.tif 000004
【００２４】
ｐ（ｚ＊’＝ｒ’｜ｚ＊＝ｒ，ｘ＊）はｔ＊がｒ’個目のユニットのｒ番目の要素から生じた確率であり、負担率と呼ばれる。ｚ，ｚ’は入力変数の各クラスタに対応する潜在変数である。各データ点が属するクラスタと各クラスタのパラメータはＥＭアルゴリズムを用いて同時に推定される。
【００２５】
本発明の実施の形態に係る時空間変数予測装置は、時空間変数の時系列データ（人口分布、人流・交通流の速度・向き、金やダイヤモンドなど鉱物資源の埋蔵量、降水量などの気象データ、土地価格など）を対象としたものであり、観測データに応じて柔軟に適用できるものである。以下では、実施の形態として、人口密度分布の時系列データが観測データとして与えられた条件の下で、未観測地点あるいは未来の時空間変数分布を推定・予測するという場合について説明する。
【００２６】
＜本発明の実施の形態に係る時空間変数予測装置の構成＞
【００２７】
次に、本発明の実施の形態に係る時空間変数予測装置の構成について説明する。図１に示すように、本発明の実施の形態に係る時空間変数予測装置１００は、ＣＰＵと、ＲＡＭと、後述する各処理ルーチンを実行するためのプログラムや各種データを記憶したＲＯＭと、を含むコンピュータで構成することが出来る。時空間変数予測装置１００は、位置情報及び時間情報を有する入力変数に対する時空間変数の観測値を有する観測データの集合に基づいて、未観測の位置情報及び時間情報での時空間変数の値を予測する。この時空間変数予測装置１００は、機能的には図１に示すように、操作部１０と、人口密度情報記憶部１２と、入力部１３と、演算部１４と、時空間変数算出部２６と、出力部２８とを備えている。操作部１０及び演算部１４は、人口密度情報記憶部１２と接続されている。
【００２８】
操作部１０は、後述する人口密度情報記憶部１２に格納されているデータに対するユーザからの各種操作を受け付ける。各種操作とは、人口密度情報記憶部１２に格納された情報を登録、修正、削除する操作等である。操作部１０の入力手段は、キーボードやマウス、メニュー画面、タッチパネルによるもの等、何でもよい。操作部１０は、マウス等の入力手段のデバイスドライバや、メニュー画面の制御ソフトウェアで実現され得る。
【００２９】
人口密度情報記憶部１２には、観測データの集合が格納されている。
【００３０】
人口密度情報記憶部１２には、後述する演算部１４が解析する観測データの集合が格納されており、演算部１４からの要求に従って、観測データの集合を読み出し、当該情報を演算部１４に送信する。
【００３１】
人口密度情報記憶部１２に格納される観測データの集合は、入力変数ｘと時空間変数ｔとの組み合わせの集合である。本発明の実施の形態では、入力変数ｘは位置及び時刻を有し、時空間変数ｔは人口密度である。ある位置、ある時刻における人口密度は｛ｘｉ，ｔｉ｝と表すことができる。人口密度情報記憶部１２はＷｅｂページを保持するＷｅｂサーバや、データベースを具備するデータベースサーバ等である。
【００３２】
図２に、観測データの集合の一例を示す。図２に示すように、例えば、位置ＩＤ及び時刻を表す情報が入力変数ｘとして格納され、人口密度を表す情報が時空間変数ｔとして格納される。
【００３３】
入力部１３は、未観測の位置情報及び時間情報を有する入力変数を受け付ける。
【００３４】
本発明の実施の形態では、入力部１３で受け付けられた位置情報及び時間情報に対して、後述する演算部１４によって得られた各パラメータに基づいて、時空間変数である人口密度の予測が行われる。
【００３５】
入力部１３の入力手段は、キーボードやマウス、メニュー画面、タッチパネルによるもの等、何でもよい。入力部１３は、マウス等の入力手段のデバイスドライバや、メニュー画面の制御ソフトウェアで実現され得る。
【００３６】
図３に、本実施の形態における入力部１３の構成例を示す。図３の構成例では、入力部１３と後述する出力部２８とが１つの画面として構成されている場合を示す。
【００３７】
図３に示すように、入力部１３は、予測を行う対象となる地点あるいは領域に関する位置情報と、予測を行う時間情報とを含む入力変数を受け付ける。また、出力部２８は、後述する時空間変数算出部２６によって出力された、入力変数に対する時空間変数の値を予測結果として表示する。
【００３８】
演算部１４は、学習部１６と、負担率パラメータ格納部２２と、ガウス過程パラメータ格納部２４とを備える。演算部１４では、学習データとして、ガウス過程のハイパーパラメータ、各ガウス過程からの負担率を表すパラメータを算出する。
【００３９】
ここで、本発明の実施の形態における原理について説明する。
【００４０】
人口密度情報記憶部１２に、時空間変数データである観測データの集合Ｄ＝｛ｘｉ，ｔｉ｝Ｎｉ＝１が格納されたとする。ここでｘｉは位置及び時刻など複数の変数を含むベクトル、ｔｉは人口密度である。
【００４１】
ここで解くべき問題は、入力部１３によって受け付けた、（１）新しい入力変数（未観測の位置及び時間）ｘ＊における人口密度の値ｔ＊の予測と、（２）未観測地点を含むある領域における人口密度の値の予測である。
【００４２】
本実施の形態では、問題（１）に焦点を当てて説明する。人口密度分布の予測面を求めるには、ｘ＊を動かしながら問題（１）を繰り返し解けばよい。
【００４３】
時空間変数の予測値ｔ＊が従うモデルが、Ｒ個のガウス過程の混合モデルで記述できると仮定する。この仮定に基づけば、新たな入力変数ｘ＊に対応する値ｔ＊の期待値を、以下の式（１）のように書き下すことができる。
【００４４】
【数４】
0006397380.tif 000005
【００４５】
ここでＸは入力変数ｘの集合（Ｘ＝｛ｘｉ｝Ｎｉ＝１）、Θ＝｛θｒｒ’｝Ｒ，Ｒ’ｒ，ｒ’＝１はｒ’番目のユニットのｒ個目のガウス過程のハイパーパラメータ、ｚ，ｚ’は入力変数の特性（地域や時間帯など）によって分類されたクラスに対応する潜在変数である。ｘ＊がｒ番目のガウス過程から生じたと仮定すると、予測値ｔ＊は次式の平均、分散を持つガウス分布に従う。
【００４６】
【数５】
0006397380.tif 000006
【００４７】
ここでｋはカーネル関数、ｋｒｒ’は要素ｋｒｒ’（ｘ＊，Ｘ）を持つベクトル、Ｃｒｒ’は要素ｋｒｒ’（ｘｎ，ｘｍ）＋β−１δｎｍ＋Ψｒｒ’ｎｍを持つＮ×Ｎの共分散行列である。また、ｔは、要素をｔｉとするベクトルである。ただしΨｒｒ’は次式で定義される対角行列である。
【００４８】
【数６】
0006397380.tif 000007
【００４９】
ここでσｎは定数である。ｉ番目のデータ点がｒ番目のガウス過程の寄与を受けないとき、すなわち負担率ｐ（ｚｉ’＝ｒ’，ｚｉ＝ｒ｜ｘｉ）のとき、Ｃｒｒ’は無限大に発散し、（Ｃｒｒ’）−１の（ｉ，ｉ）成分は０となる。すなわちｔｉはｒ’個目のユニットのｒ番目のガウス過程のパラメータ推定に影響しない。
【００５０】
従来手法では、カーネル関数ｋ（ｘｎ，ｘｍ）は予測の前処理段階で一意に決められる。一方、本発明の実施の形態では、カーネル関数を自由に設計し、各々のカーネル関数のパラメータと重みとを自動推定する。本発明の実施の形態では、気象データや環境データなどの時系列地理統計解析で一般的に使われる代表的なコバリオグラムを二つ挙げる。
【００５１】
【数７】
0006397380.tif 000008
【００５２】
ここで入力変数ｘの次元数をｎｘとおくと
【００５３】
【数８】
0006397380.tif 000009
【００５４】
である。また、上記のコバリオグラムに対応するカーネル関数は次式のように書き下せる。
【００５５】
【数９】
0006397380.tif 000010
【００５６】
本発明の実施の形態では、上記のカーネル関数の線形和で定義した、以下の式で示す新たなカーネル関数を用いる。
【００５７】
【数１０】
0006397380.tif 000011
【００５８】
ここで推定すべきガウス過程のハイパーパラメータは
【００５９】
【数１１】
0006397380.tif 000012
【００６０】
である。入力変数の集合Ｘ＝｛ｘｉ｝Ｎｉ＝１、およびガウス過程のハイパーパラメータの集合Θ＝（θ１１，・・・，θＲＲ’）が与えられた条件の下で、観測データの集合Ｄの尤度関数は次式のように書き下せる。
【００６１】
【数１２】
0006397380.tif 000013
【００６２】
ここでｔ＼ｉはｔからｉ番目の要素を除いたものを表す。
【００６３】
本発明の実施の形態では、上記式の尤度関数を最大化するパラメータΘ，ｚ，ｚ’を推定するため、ＥＭアルゴリズムを用いる。ＥＭアルゴリズムでは、次の手順［１］〜［４］でパラメータ推定を行う。なお、本発明の実施の形態では、負担率の事前分布として以下の式（４）〜（５）に示すＳｏｆｔｍａｘ関数を採用し、後述する時空間変数算出部２６による時空間変数の予測で用いられるＳｏｆｔｍａｘ関数のパラメータを算出する。
【００６４】
【数１３】
0006397380.tif 000014
【００６５】
［１］各パラメータの初期値を選択する。
【００６６】
［２］ＥＭアルゴリズムのＥステップにおいて、観測データ集合の各観測データに対し、以下の式に従って、当該観測データの入力変数に対する、ｒ’番目のユニットのｒ番目のガウス過程の負担率ｐ（ｚｉ＝ｒ｜ｚｉ’＝ｒ’，Ｘ，ｔ，θｒｒ’）、当該観測データの入力変数に対する、ｒ番目のガウス過程の負担率ｐ（ｚｉ＝ｒ｜Ｘ，ｔ，θｒｒ’）、当該観測データの入力変数に対する、ｒ’番目のユニットのｒ番目のガウス過程の同時確率ｐ（ｚｉ＝ｒ，ｚｉ’＝ｒ’｜Ｘ，ｔ，θｒｒ’）を計算する。
【００６７】
【数１４】
0006397380.tif 000015
【００６８】
なお、上記式（２）におけるｐ（ｔｉ｜Ｘ，ｔ，θｒｒ’）は、次式に示す平均μｉｒｒ’、分散Σｉｒｒ’を持つガウス分布で表される。
【００６９】
【数１５】
0006397380.tif 000016
【００７０】
［３］ＥＭアルゴリズムのＭステップにおいて、観測データ集合の各観測データに対して計算された、ｒ’番目のユニットのｒ番目のガウス過程の負担率ｐ（ｚｉ＝ｒ｜ｚｉ’＝ｒ’，Ｘ，ｔ，θｒｒ’）、ｒ番目のガウス過程の負担率ｐ（ｚｉ＝ｒ｜Ｘ，ｔ，θｒｒ’）、ｒ’番目のユニットのｒ番目のガウス過程の同時確率ｐ（ｚｉ＝ｒ，ｚｉ’＝ｒ’｜Ｘ，ｔ，θｒｒ’）に基づいて、尤度関数Ｑを最大化するガウス過程のパラメータΘｎｅｗを計算する。ここで尤度関数Ｑを次式で表されるＱ〜で予め近似する。
【００７１】
【数１６】
0006397380.tif 000017
【００７２】
ここで、πｉｒｒ’は、
【００７３】
【数１７】
0006397380.tif 000018
【００７４】
で定義される。
【００７５】
また、目的関数Ｑ〜のガウス過程のハイパーパラメータθｒｒ’に関する微分は次式で書き下せる。
【００７６】
【数１８】
0006397380.tif 000019
【００７７】
ここでθｒｒ’ｊはｒ番目のガウス過程のｊ番目のハイパーパラメータである。降下勾配法、準ニュートン法等を用いれば目的関数Ｑ〜を最大化するガウス過程のパラメータΘを得ることができる。更新式は次式のように書き下せる。
【００７８】
【数１９】
0006397380.tif 000020
【００７９】
また、負担率を表すＳｏｆｔｍａｘ関数のパラメータｖｒ及びｖｒｒ’の更新式は次式のように書き下せる。
【００８０】
【数２０】
0006397380.tif 000021
【００８１】
ここでηは定数である。
【００８２】
［４］ＥＭアルゴリズムの収束条件が満たされているか調べ、満たされていなければ
【００８３】
【数２１】
0006397380.tif 000022
【００８４】
を実行し、上記手順［２］に戻る。収束条件が満たされていれば、処理を終了する。
【００８５】
従って、学習部１６は、人口密度情報記憶部１２に格納された観測データの集合に基づいて、観測データ同士の類似性を定義する関数であるカーネル関数の各々のハイパーパラメータと、観測データの各々に対する、複数のガウス過程の各々の寄与度を表すパラメータである負担率とを学習する。
【００８６】
本発明の実施の形態におけるカーネル関数は、複数のガウス過程を、空間的な広がり及び時間的な広がりに対応する複数の階層で混合した階層混合ガウス過程でモデル化される。また、本発明の実施の形態におけるカーネル関数は、入力変数に対する時空間変数の値を予測するためのモデルに含まれる、複数のガウス過程の各々についての、観測データ同士の類似性を定義する関数である。
【００８７】
学習部１６は、負担率推定部１８と、ガウス過程パラメータ推定部２０と、反復判定部２１とを備える。
【００８８】
負担率推定部１８は、人口密度情報記憶部１２に格納された観測データの集合と、複数のガウス過程のカーネル関数の各々のハイパーパラメータの初期値、又はガウス過程パラメータ推定部２０による前回推定値とに基づいて、上記式（６）〜（８）に従って、観測データの各々に対する、複数のガウス過程からなる複数のユニットの各々の寄与度を表すパラメータであるユニット負担率、及び観測データの各々に対する、複数のガウス過程の各々の寄与度を表すパラメータである負担率を推定する。
【００８９】
ガウス過程パラメータ推定部２０は、人口密度情報記憶部１２に格納された観測データの集合と、負担率推定部１８によって推定された、観測データの各々に対する、複数のユニットの各々のユニット負担率、及び観測データの各々に対する、複数のガウス過程の各々の負担率とに基づいて、上記式（９）〜（１０）、（１３）に従って、複数のガウス過程の各々に対し、ガウス過程のカーネル関数の各々のハイパーパラメータを推定する。
【００９０】
また、ガウス過程パラメータ推定部２０は、人口密度情報記憶部１２に格納された観測データの集合に基づいて、上記式（１１）〜（１２）に従って、Ｓｏｆｔｍａｘ関数のパラメータｖｒ及びｖｒｒ’を推定する。
【００９１】
反復判定部２１は、予め定められた反復終了条件を満たすまで、負担率推定部１８による推定、及びガウス過程パラメータ推定部２０による推定を繰り返す。そして、反復判定部２１は、予め定められた反復終了条件が満たされた場合には、ガウス過程パラメータ推定部２０によって推定されたＳｏｆｔｍａｘ関数のパラメータｖｒ及びｖｒｒ’を負担率パラメータ格納部２２に格納し、ガウス過程パラメータ推定部２０によって推定されたハイパーパラメータをガウス過程パラメータ格納部２４に格納する。
【００９２】
負担率パラメータ格納部２２には、ガウス過程パラメータ推定部２０によって推定された、Ｓｏｆｔｍａｘ関数のパラメータｖｒ及びｖｒｒ’が格納される。
【００９３】
ガウス過程パラメータ格納部２４には、ガウス過程パラメータ推定部２０によって推定されたガウス過程のカーネル関数の各々のハイパーパラメータが格納される。
【００９４】
時空間変数算出部２６は、入力部１３によって受け付けた、未観測の位置情報及び時間情報を有する入力変数と、負担率パラメータ格納部２２に格納されたＳｏｆｔｍａｘ関数のパラメータとに基づいて、入力変数に対する複数のガウス過程の各々の寄与度を表すパラメータである負担率を推定する。
【００９５】
具体的には、時空間変数算出部２６は、入力変数ｘ＊と、負担率パラメータ格納部２２に格納されたＳｏｆｔｍａｘ関数のパラメータｖｒ及びｖｒｒ’とに基づいて、以下の式（１４）〜（１５）に従って、入力変数ｘ＊に対する複数のガウス過程の各々の寄与度を表すパラメータである負担率を推定する。
【００９６】
【数２２】
0006397380.tif 000023
【００９７】
そして、時空間変数算出部２６は、推定した入力変数ｘ＊に対する負担率と、ガウス過程パラメータ格納部２４に格納されたハイパーパラメータに基づいて、上記式（１）〜（３）に従って、入力変数に対する時空間変数の値を予測する。
【００９８】
出力部２８は、時空間変数算出部２６によって予測された、入力された入力変数に対する時空間変数の値を、結果として出力する。
【００９９】
例えば、出力部２８は、上記図３に示すように、予測対象の地点あるいは領域及び時間での時空間変数の予測値である人口密度と、当該人口密度に関連する情報である混雑度とを結果として出力する。
【０１００】
＜本発明の実施の形態に係る時空間変数予測装置の作用＞
【０１０１】
次に、本発明の実施の形態に係る時空間変数予測装置１００の作用について説明する。
【０１０２】
＜学習処理ルーチン＞
まず、時空間変数予測装置１００は、操作部１０より観測データの集合Ｄが入力されると、観測データの集合Ｄを人口密度情報記憶部１２に格納する。そして、時空間変数予測装置１００は、図４に示す学習処理ルーチンを実行する。
【０１０３】
まず、ステップＳ１００では、繰り返し変数ｊと各パラメータとを初期化する。
【０１０４】
次に、ステップＳ１０２において、反復判定部２１は、予め定められた反復終了条件を満たしたか否かを判定する。予め定められた条件として、繰り返し変数ｊが予め定められた値Ｎｉｔｅｒ未満である場合には、ステップＳ１１０へ進む。一方、繰り返し変数ｊがＮｉｔｅｒ以上である場合には、ステップＳ１０４へ進む。
【０１０５】
ステップＳ１０４において、負担率推定部１８は、人口密度情報記憶部１２に格納された観測データの集合と、上記ステップＳ１００で設定された、複数のガウス過程のカーネル関数の各々のハイパーパラメータの初期値、又は後述するステップＳ１０６で前回推定されたガウス過程のカーネル関数の各々のハイパーパラメータとに基づいて、上記式（６）〜（８）に従って、観測データの各々に対するユニット負担率、及び観測データの各々に対する負担率を推定する。
【０１０６】
ステップＳ１０６において、ガウス過程パラメータ推定部２０は、人口密度情報記憶部１２に格納された観測データの集合と、上記ステップＳ１０４で推定された、観測データの各々に対するユニット負担率、及び観測データの各々に対する負担率とに基づいて、上記式（９）〜（１０）、（１３）に従って、複数のガウス過程の各々に対し、ガウス過程のカーネル関数の各々のハイパーパラメータを推定する。また、ガウス過程パラメータ推定部２０は、人口密度情報記憶部１２に格納された観測データの集合に基づいて、上記式（１１）〜（１２）に従って、Ｓｏｆｔｍａｘ関数のパラメータを推定する。
【０１０７】
ステップＳ１０８では、繰り返し変数ｊを１インクリメントしてステップＳ１０２へ戻る。
【０１０８】
ステップＳ１１０において、反復判定部２１は、上記ステップＳ１０６で推定されたＳｏｆｔｍａｘ関数のパラメータを負担率パラメータ格納部２２に格納し、推定されたガウス過程のハイパーパラメータをガウス過程パラメータ格納部２４に格納して、学習処理ルーチンを終了する。
【０１０９】
次に、図５に示す時空間変数算出処理ルーチンについて説明する。
【０１１０】
学習処理ルーチンが実行され、負担率パラメータ格納部２２にＳｏｆｔｍａｘ関数のパラメータｖｒｒ’が格納され、ガウス過程パラメータ格納部２４にガウス過程のハイパーパラメータが格納され、予測を行う対象となる未観測の位置情報及び時間情報を含む入力変数が入力されると、時空間変数予測装置１００は、図５に示す時空間変数算出処理ルーチンを実行する。
【０１１１】
＜時空間変数算出処理ルーチン＞
ステップＳ２００において、入力部１３は、未観測の位置情報及び時間情報を含む入力変数を受け付ける。
【０１１２】
ステップＳ２０２において、時空間変数算出部２６は、負担率パラメータ格納部２２に格納されたＳｏｆｔｍａｘ関数のパラメータと、ガウス過程パラメータ格納部２４に格納されたガウス過程のハイパーパラメータとを読み出す。
【０１１３】
ステップＳ２０４において、時空間変数算出部２６は、上記ステップＳ２００で受け付けた入力変数と、上記ステップＳ２０２で読み込まれたＳｏｆｔｍａｘ関数のパラメータとに基づいて、上記式（１４）〜（１５）に従って、入力変数に対する複数のガウス過程の各々の寄与度を表すパラメータである負担率を推定する。そして、時空間変数算出部２６は、推定された入力変数に対する負担率と、上記ステップＳ２０２で読み込まれたガウス過程のハイパーパラメータとに基づいて、上記式（１）〜（３）に従って、入力変数に対する時空間変数の値である人口密度を予測する。
【０１１４】
出力部２８は、時空間変数算出部２６によって予測された、入力変数に対する時空間変数の値である人口密度を結果として出力して、時空間変数算出処理ルーチンを終了する。
【０１１５】
以上説明したように、本発明の実施の形態に係る時空間変数予測装置によれば、位置情報及び時間情報を有する入力変数に対する時空間変数の観測値を有する観測データの集合に基づいて、複数のガウス過程を、空間的な広がり及び時間的な広がりに対応する複数の階層で混合した階層混合ガウス過程でモデル化された、入力変数に対する時空間変数の値を予測するためのモデルに含まれる、複数のガウス過程の各々についての、観測データ同士の類似性を定義する関数であるカーネル関数の各々のハイパーパラメータと、観測データの各々に対する、複数のガウス過程の各々の寄与度を表すパラメータである負担率とを学習することにより、時間的及び空間的相関を持つ時空間変数の値を精度よく予測することができる。
【０１１６】
また、時間的及び空間的相関を持つ時空間変数の値を正確に予測することができる。
【０１１７】
なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【０１１８】
例えば、上記の実施の形態では、人口密度分布の時系列データが観測データとして与えられた条件の下で、未観測地点あるいは未来の時空間変数分布を推定・予測するという場合について説明したが、他の様々な時空間変数の時系列データ（人口分布、人流・交通流の速度・向き、金やダイヤモンドなど鉱物資源の埋蔵量、降水量などの気象データ、土地価格など）を対象として本発明を適用することができる。
【０１１９】
また、本発明の実施の形態においては、時空間変数予測装置１００によって各パラメータを推定し、推定された各パラメータを用いて、入力変数に対応する時空間変数である人口密度を予測する場合を例に説明したが、これに限定されるものではない。例えば、各パラメータを推定する処理と、推定された各パラメータを用いて時空間変数を予測する処理とを別々の装置として構成してもよい。この場合、各パラメータを推定する装置は、操作部１０と、人口密度情報記憶部１２と、演算部１４とを備え、推定された各パラメータを用いて時空間変数を予測する装置は、入力部１３と、時空間変数算出部２６と、出力部２８とを備える。
【０１２０】
また、本発明の実施の形態では、負担率の事前分布としてＳｏｆｔｍａｘ関数を採用し、Ｓｏｆｔｍａｘ関数のパラメータを推定する場合を例に説明したが、これに限定されるものではなく、他の関数を用いてもよい。
【０１２１】
また、上述の時空間変数予測装置１００は、内部にコンピュータシステムを有しているが、「コンピュータシステム」は、ＷＷＷシステムを利用している場合であれば、ホームページ提供環境（あるいは表示環境）も含むものとする。
【０１２２】
また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能であるし、ネットワークを介して提供することも可能である。
【符号の説明】
【０１２３】
１０   操作部
１２   人口密度情報記憶部
１３   入力部
１４   演算部
１６   学習部
１８   負担率推定部
２０   ガウス過程パラメータ推定部
２１   反復判定部
２２   負担率パラメータ格納部
２４   ガウス過程パラメータ格納部
２６   時空間変数算出部
２８   出力部
１００ 時空間変数予測装置
【図１】
0006397380.tif 000024
【図２】
0006397380.tif 000025
【図３】
0006397380.tif 000026
【図４】
0006397380.tif 000027
【図５】
0006397380.tif 000028