(57)【要約】【課題】ネットワークを構成する複数の要素のうち未観測の要素のコストを適切かつ効率よく推定することができる、推定装置、推定方法、及びプログラムを提供する。【解決手段】推定装置１００は、空間及び時間に関するモードに対応する複数の列ベクトルを用いて表されたモデルを用いる。推定装置１００は、データ受付部１０及び要素コスト算出部１２を備える。データ受付部１０には、観測データ及び列ベクトルごとに設定された事前分布のデータが入力される。要素コスト算出部１２は、尤度関数、観測データ及び事前分布のデータに基づいて、列ベクトルごとに条件付事後分布を生成する。また、要素コスト算出部１２は、複数回のサンプリングを実行し、サンプリングごとに、各列ベクトルの条件付事後分布から列ベクトルを抽出する。更に、要素コスト算出部１２は、抽出した列ベクトルと上述のモデルとを用いて、複数の要素のコストの推定値を算出する。【選択図】図１
2018112884.tif 000002
【特許請求の範囲】
【請求項１】
ネットワークを構成する複数の要素のコストを、予め設定されたモデルを用いて推定する推定装置であって、
前記複数の要素のコストは、時間に関するモード及び空間に関するモードを少なくとも含む複数のモードを有するテンソルで示され、
前記モデルは、前記複数のモードにそれぞれ対応する複数の列ベクトルを用いて表され、
前記ネットワークを構成する複数の要素のうちの所定数の要素のコストの観測データの入力、及び前記列ベクトルごとに予め設定された事前分布に関するデータの入力を受け付ける、データ受付部と、
予め設定された尤度関数、並びに前記データ受付部が受け付けた前記観測データ及び前記事前分布に関する前記データに基づいて前記列ベクトルごとに条件付事後分布を生成し、複数回のサンプリングを実行することによって、サンプリングごとに、各列ベクトルの前記条件付事後分布から列ベクトルを抽出し、抽出した列ベクトルと前記モデルとを用いて、前記複数の要素のコストの推定値を算出する、要素コスト算出部と、
を備えることを特徴とする、推定装置。
【請求項２】
前記要素コスト算出部は、要素ごとに前記複数回のサンプリングに基づいて算出された複数の前記推定値の平均値を算出し、算出した平均値を各要素のコストとして出力する、
請求項１に記載の推定装置。
【請求項３】
前記要素コスト算出部は、前記サンプリングとして、ギブスサンプリングを実行する、
請求項１または２に記載の推定装置。
【請求項４】
前記空間に関するモードに対応する前記列ベクトルの前記事前分布の分散は、前記ネットワークのトポロジーに基づいて設定される、
請求項１から３のいずれかに記載の推定装置。
【請求項５】
前記テンソルは、前記時間に関するモード及び前記空間に関するモードを含む３つのモードを有する３階のテンソルであり、
前記モデルは、係数と前記３つのモードにそれぞれ対応する列ベクトルとを用いて表され、
前記データ受付部は、更に、前記係数に対して予め設定された事前分布に関するデータの入力を受け付け、
前記要素コスト算出部は、更に、前記尤度関数、並びに前記データ受付部が受け付けた前記観測データ及び前記係数に対して設定された前記事前分布に関する前記データに基づいて前記係数の条件付事後分布を生成し、複数回のサンプリングを実行することによって、サンプリングごとに、前記係数の前記条件付事後分布から係数を抽出し、抽出した係数を用いて、前記複数の要素のコストの前記推定値を算出する、
請求項１から４のいずれかに記載の推定装置。
【請求項６】
ネットワークを構成する複数の要素のコストを、予め設定されたモデルを用いて推定する方法であって、
前記複数の要素のコストは、時間に関するモード及び空間に関するモードを少なくとも含む複数のモードを有するテンソルで示され、
前記モデルは、前記複数のモードにそれぞれ対応する複数の列ベクトルを用いて表され、
（ａ）前記ネットワークを構成する複数の要素のうちの所定数の要素のコストの観測データの入力、及び前記列ベクトルごとに予め設定された事前分布に関するデータの入力を受け付ける、ステップと、
（ｂ）予め設定された尤度関数、並びに前記（ａ）のステップで受け付けた前記観測データ及び前記事前分布に関する前記データに基づいて前記列ベクトルごとに条件付事後分布を生成する、ステップと、
（ｃ）複数回のサンプリングを実行することによって、サンプリングごとに、前記（ｂ）のステップで生成した各列ベクトルの前記条件付事後分布から列ベクトルを抽出する、ステップと、
（ｄ）前記（ｃ）のステップで抽出した列ベクトルと前記モデルとを用いて、前記複数の要素のコストの推定値を算出する、ステップと、
を備えることを特徴とする、推定方法。
【請求項７】
前記（ｄ）のステップでは、要素ごとに前記複数回のサンプリングに基づいて算出された複数の前記推定値の平均値を算出し、算出した平均値を各要素のコストとして出力する、
請求項６に記載の推定方法。
【請求項８】
前記（ｃ）のステップでは、前記サンプリングとして、ギブスサンプリングを実行する、
請求項６または７に記載の推定方法。
【請求項９】
前記空間に関するモードに対応する前記列ベクトルの前記事前分布の分散は、前記ネットワークのトポロジーに基づいて設定される、
請求項６から８のいずれかに記載の推定方法。
【請求項１０】
前記テンソルは、前記時間に関するモード及び前記空間に関するモードを含む３つのモードを有する３階のテンソルであり、
前記モデルは、係数と前記３つのモードにそれぞれ対応する列ベクトルとを用いて表され、
前記（ａ）のステップでは、更に、前記係数に対して予め設定された事前分布に関するデータの入力を受け付け、
前記（ｂ）のステップでは、更に、前記尤度関数、並びに前記（ａ）のステップで受け付けた前記観測データ及び前記係数に対して設定された前記事前分布に関する前記データに基づいて前記係数の条件付事後分布を生成し、
前記（ｃ）のステップでは、更に、複数回のサンプリングを実行することによって、サンプリングごとに、前記（ｂ）のステップで生成した前記係数の前記条件付事後分布から係数を抽出し、
前記（ｄ）のステップでは、前記（ｃ）のステップで抽出した係数を用いて、前記複数の要素のコストの前記推定値を算出する、
請求項６から９のいずれかに記載の推定方法。
【請求項１１】
コンピュータに、ネットワークを構成する複数の要素のコストを、予め設定されたモデルを用いて推定させるプログラムであって、
前記複数の要素のコストは、時間に関するモード及び空間に関するモードを少なくとも含む複数のモードを有するテンソルで示され、
前記モデルは、前記複数のモードにそれぞれ対応する複数の列ベクトルを用いて表され、
前記コンピュータに、
（ａ）前記ネットワークを構成する複数の要素のうちの所定数の要素のコストの観測データの入力、及び前記列ベクトルごとに予め設定された事前分布に関するデータの入力を受け付ける、ステップと、
（ｂ）予め設定された尤度関数、並びに前記（ａ）のステップで受け付けた前記観測データ及び前記事前分布に関する前記データに基づいて前記列ベクトルごとに条件付事後分布を生成する、ステップと、
（ｃ）複数回のサンプリングを実行することによって、サンプリングごとに、前記（ｂ）のステップで生成した各列ベクトルの前記条件付事後分布から列ベクトルを抽出する、ステップと、
（ｄ）前記（ｃ）のステップで抽出した列ベクトルと前記モデルとを用いて、前記複数の要素のコストの推定値を算出する、ステップと、
を実行させることを特徴とする、プログラム。
【請求項１２】
前記（ｄ）のステップでは、要素ごとに前記複数回のサンプリングに基づいて算出された複数の前記推定値の平均値を算出し、算出した平均値を各要素のコストとして出力する、
請求項１１に記載のプログラム。
【請求項１３】
前記（ｃ）のステップでは、前記サンプリングとして、ギブスサンプリングを実行する、
請求項１１または１２に記載のプログラム。
【請求項１４】
前記空間に関するモードに対応する前記列ベクトルの前記事前分布の分散は、前記ネットワークのトポロジーに基づいて設定される、
請求項１１から１３のいずれかに記載のプログラム。
【請求項１５】
前記テンソルは、前記時間に関するモード及び前記空間に関するモードを含む３つのモードを有する３階のテンソルであり、
前記モデルは、係数と前記３つのモードにそれぞれ対応する列ベクトルとを用いて表され、
前記（ａ）のステップでは、更に、前記係数に対して予め設定された事前分布に関するデータの入力を受け付け、
前記（ｂ）のステップでは、更に、前記尤度関数、並びに前記（ａ）のステップで受け付けた前記観測データ及び前記係数に対して設定された前記事前分布に関する前記データに基づいて前記係数の条件付事後分布を生成し、
前記（ｃ）のステップでは、更に、複数回のサンプリングを実行することによって、サンプリングごとに、前記（ｂ）のステップで生成した前記係数の前記条件付事後分布から係数を抽出し、
前記（ｄ）のステップでは、前記（ｃ）のステップで抽出した係数を用いて、前記複数の要素のコストの前記推定値を算出する、
請求項１１から１４のいずれかに記載のプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、ネットワークを構成する複数の要素（ノード及びリンク等）のコストを推定することができる、推定装置、推定方法、及びプログラムに関する。
【背景技術】
【０００２】
交通網及び配水網等のネットワークの最適化においては、ネットワークを構成する要素（ノード及びリンク）のコスト（重み）を明らかにする必要がある。しかしながら、大規模ネットワークでは、ネットワークを構成する要素の数が多くなるので、コストを観測するためのセンサを全ての要素に設置するのは困難である。
【０００３】
このため、大規模ネットワークの最適化においては、センサによって観測することができない要素のコストを推定する必要がある。このような、未観測の要素のコストを推定するための技術が、例えば、特許文献１に開示されている。
【０００４】
特許文献１に開示された交通情報推定システムでは、センサによって交通情報が観測された道路リンクにおけるリンク旅行速度に基づいて、交通情報が未観測の道路リンクにおけるリンク旅行速度を推定する。具体的には、例えば、観測済みの道路リンクと相関度が高い未観測の道路リンクのリンク旅行速度が、前記観測済みの道路リンクのリンク旅行速度に基づいて推測される。上記相関度は、交通状態が類似する道路リンク同士程高くなるように計算される。
【０００５】
また、特許文献１においては、今回のタイムスパンでは交通情報が観測されていない道路リンクについても、過去に観測された交通情報に基づいて、リンク旅行速度を推定することができる。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】特開２０１０−１９１６１４号公報
【非特許文献】
【０００７】
【非特許文献１】Rai, P., Wang, Y., Guo, S., Chen, G., Dunson, D. and Carin, L. (2014). Scalable Bayesian Low-Rank Decomposition of Incomplete Multiway Tensors. ICML.
【発明の概要】
【発明が解決しようとする課題】
【０００８】
ところで、特許文献１に記載されたシステムでは、タイムスパンごとに道路リンクのリンク旅行速度を推定する必要がある。このため、各道路リンクのリンク旅行速度を複数のタイムスパンに亘って推定する場合には、推定のための処理回数が多くなり、リンク旅行速度を効率よく推定することができない。
【０００９】
また、上記のように、特許文献１に記載されたシステムでは、交通情報が観測されていない道路リンクについては、過去に観測された交通情報に基づいてリンク旅行速度を推定する必要がある。このため、過去に交通情報が一度も観測されていない道路リンクについては、リンク旅行速度を適切に推定することができない。
【００１０】
本発明の目的の一例は、ネットワークを構成する複数の要素のうち未観測の要素のコストを適切かつ効率よく推定することができる、推定装置、推定方法、及びプログラムを提供することにある。
【課題を解決するための手段】
【００１１】
上記目的を達成するため、本発明の一側面における推定装置は、ネットワークを構成する複数の要素のコストを、予め設定されたモデルを用いて推定する推定装置であって、
前記複数の要素のコストは、時間に関するモード及び空間に関するモードを少なくとも含む複数のモードを有するテンソルで示され、
前記モデルは、前記複数のモードにそれぞれ対応する複数の列ベクトルを用いて表され、
前記ネットワークを構成する複数の要素のうちの所定数の要素のコストの観測データの入力、及び前記列ベクトルごとに予め設定された事前分布に関するデータの入力を受け付ける、データ受付部と、
予め設定された尤度関数、並びに前記データ受付部が受け付けた前記観測データ及び前記事前分布に関する前記データに基づいて前記列ベクトルごとに条件付事後分布を生成し、複数回のサンプリングを実行することによって、サンプリングごとに、各列ベクトルの前記条件付事後分布から列ベクトルを抽出し、抽出した列ベクトルと前記モデルとを用いて、前記複数の要素のコストの推定値を算出する、要素コスト算出部と、
を備えることを特徴とする。
【００１２】
また、上記目的を達成するため、本発明の一側面における推定方法は、ネットワークを構成する複数の要素のコストを、予め設定されたモデルを用いて推定する方法であって、
前記複数の要素のコストは、時間に関するモード及び空間に関するモードを少なくとも含む複数のモードを有するテンソルで示され、
前記モデルは、前記複数のモードにそれぞれ対応する複数の列ベクトルを用いて表され、
（ａ）前記ネットワークを構成する複数の要素のうちの所定数の要素のコストの観測データの入力、及び前記列ベクトルごとに予め設定された事前分布に関するデータの入力を受け付ける、ステップと、
（ｂ）予め設定された尤度関数、並びに前記（ａ）のステップで受け付けた前記観測データ及び前記事前分布に関する前記データに基づいて前記列ベクトルごとに条件付事後分布を生成する、ステップと、
（ｃ）複数回のサンプリングを実行することによって、サンプリングごとに、前記（ｂ）のステップで生成した各列ベクトルの前記条件付事後分布から列ベクトルを抽出する、ステップと、
（ｄ）前記（ｃ）のステップで抽出した列ベクトルと前記モデルとを用いて、前記複数の要素のコストの推定値を算出する、ステップと、
を備えることを特徴とする。
【００１３】
更に、上記目的を達成するため、本発明の一側面におけるプログラムは、コンピュータに、ネットワークを構成する複数の要素のコストを、予め設定されたモデルを用いて推定させるプログラムであって、
前記複数の要素のコストは、時間に関するモード及び空間に関するモードを少なくとも含む複数のモードを有するテンソルで示され、
前記モデルは、前記複数のモードにそれぞれ対応する複数の列ベクトルを用いて表され、
前記コンピュータに、
（ａ）前記ネットワークを構成する複数の要素のうちの所定数の要素のコストの観測データの入力、及び前記列ベクトルごとに予め設定された事前分布に関するデータの入力を受け付ける、ステップと、
（ｂ）予め設定された尤度関数、並びに前記（ａ）のステップで受け付けた前記観測データ及び前記事前分布に関する前記データに基づいて前記列ベクトルごとに条件付事後分布を生成する、ステップと、
（ｃ）複数回のサンプリングを実行することによって、サンプリングごとに、前記（ｂ）のステップで生成した各列ベクトルの前記条件付事後分布から列ベクトルを抽出する、ステップと、
（ｄ）前記（ｃ）のステップで抽出した列ベクトルと前記モデルとを用いて、前記複数の要素のコストの推定値を算出する、ステップと、
を実行させることを特徴とする。
【発明の効果】
【００１４】
以上のように本発明によれば、ネットワークを構成する複数の要素のうち未観測の要素のコストを適切かつ効率よく推定することができる。
【図面の簡単な説明】
【００１５】
【図１】図１は、本発明の実施の形態における推定装置の概略構成を示すブロック図である。
【図２】図２は、本発明の実施の形態における推定装置の具体的構成を示すブロック図である。
【図３】図３は、ネットワークの一例を示す図である。
【図４】図４は、本発明の実施の形態における推定装置の動作を示すフロー図である。
【図５】図５は、本発明の実施の形態における推定装置を実現するコンピュータの一例を示すブロック図である。
【発明を実施するための形態】
【００１６】
（実施の形態）
以下、本発明の実施の形態における、推定装置、推定方法、及びプログラムについて、図１〜５を参照しながら説明する。
【００１７】
［装置構成］
最初に、図１を用いて、本発明の実施の形態における推定装置の概略構成について説明する。図１は、本発明の実施の形態における推定装置の概略構成を示すブロック図である。
【００１８】
図１に示すように、本実施の形態における推定装置１００は、データ受付部１０及び要素コスト算出部１２を備えている。推定装置１００は、ネットワークを構成する複数の要素のコストを推定する装置である。推定装置１００は、例えば、交通網、配水網、及び航空網等の種々のネットワークにおいて、複数の要素のコストを推定するために利用される。
【００１９】
本実施形態では、推定装置１００は、ネットワークを構成する複数の要素のコストを、予め設定されたモデルを用いて推定する。複数の要素のコストは、空間に関するモード及び時間に関するモードを少なくとも含む複数のモードを有するテンソルで示される。予め設定された上記モデルは、上記複数のモードにそれぞれ対応する複数の列ベクトルを用いて表される。モデルについての詳細は後述する。
【００２０】
データ受付部１０は、ネットワークを構成する複数の要素のうちの所定数の要素のコストの観測データの入力、及び上述の列ベクトルごとに予め設定された事前分布に関するデータの入力を受け付ける。観測データは、例えば、上記所定数の要素ごとに設けられたセンサ（図示せず）によって観測されたデータである。観測データは、各モードに対応する情報に関連付けられている。本実施形態では、観測データは、観測位置（ネットワークにおける要素の位置）及び観測時間に関連付けられている。データ受付部１０は、受け付けた観測データ及び列ベクトルごとの事前分布に関するデータを、要素コスト算出部１２へ送信する。
【００２１】
要素コスト算出部１２は、予め設定された尤度関数、並びにデータ受付部１０から入力された観測データ及び列ベクトルごとの事前分布に関するデータに基づいて、列ベクトルごとに条件付事後分布を生成する。また、要素コスト算出部１２は、複数回のサンプリングを実行することによって、サンプリングごとに、各列ベクトルの条件付事後分布から列ベクトルを抽出する。更に、要素コスト算出部１２は、抽出した列ベクトルと予め設定された上述のモデルとを用いて、複数の要素のコストの推定値を算出する。
【００２２】
以上のように、本実施形態に係る推定装置１００においては、空間に関するモードに対応する列ベクトル及び時間に関するモードに対応する列ベクトルを用いて表されたモデルが用いられる。また、列ベクトルごとに事前分布に関するデータが与えられ、その事前分布に基づいて条件付事後分布が生成され、複数回のサンプリングによって、サンプリングごとに、各列ベクトルの条件付事後分布から列ベクトルが抽出される。そして、抽出された列ベクトルと上記モデルとを用いて、複数の要素のコストの推定値が算出される。したがって、本実施形態に係る推定装置１００においては、複数の異なる観測位置にまたがって、各要素のコストを推定することができる。これにより、例えば、過去に一度もデータが観測されていない要素についても、コストを適切に推定することができる。また、本実施形態に係る推定装置１００においては、複数の異なる観測時間にまたがって、各要素のコストを推定することができる。これにより、観測時間ごとに各要素のコストを推定しなくてよいので、各要素のコストを効率よく推定することができる。
【００２３】
続いて、図２を用いて、本発明の実施の形態における推定装置の構成及び機能について更に具体的に説明する。図２には、本発明の実施の形態における推定装置の具体的構成を示すブロック図が示されている。
【００２４】
図２に示すように、本実施形態では、推定装置１００は、上述のデータ受付部１０及び要素コスト算出部１２を備えている。本実施形態では、データ受付部１０は、観測データ受付部１４と、事前分布受付部１６とを備えている。要素コスト算出部１２は、推定部１８と、推定値算出部２０と、コスト表示部２２とを備えている。推定部１８は、因子行列推定部２４と、時空間相関推定部２６とを備えている。以下においては、まず、ネットワーク、モデル及び尤度関数について説明し、データ受付部１０及び要素コスト算出部１２の詳細については後述する。
【００２５】
図３は、推定装置１００によるコスト推定が行なわれるネットワークの一例を示す図である。図３に示すネットワーク３０は、複数のノードｎ１〜ｎ７と、複数のリンクｆ１〜ｆ６とを備える。例えば、ネットワーク３０が配水網を表している場合には、ノードｎ１〜ｎ７がそれぞれ、浄水場又は配水場等の施設に対応し、リンクｆ１〜ｆ６がそれぞれ、施設間において水が流れる管路に対応する。また、例えば、ネットワーク３０が交通網を表している場合には、ノードｎ１〜ｎ７がそれぞれ、地点（交差点等）に対応し、リンクｆ１〜ｆ６がそれぞれ、２つの地点を結ぶ交通路に対応する。
【００２６】
なお、図３においてリンクｆ１〜ｆ６に付された矢印は、ネットワーク３０が配水網であると仮定した場合には水の流通方向を示し、ネットワーク３０が交通網であると仮定した場合には車の進行方向を示す。
【００２７】
本実施形態では、複数のノードをネットワークを構成する複数の要素として、推定装置１００によって各要素のコストを推定してもよく、複数のリンク（エッジ）をネットワークを構成する複数の要素として、推定装置１００によって各要素のコストを推定してもよい。
【００２８】
例えば、ネットワーク３０が配水網である場合には、推定装置１００によって、各管路（リンクｆ１〜ｆ６）における水の流量を、各要素のコストとして推定してもよい。また、例えば、ネットワーク３０が交通網である場合、推定装置１００によって、車が各交通路（リンクｆ１〜ｆ６）を通るのに要する時間（２つの地点間の移動時間）を、各要素のコストとして推定してもよい。更に、例えば、ネットワーク３０が交通網である場合、推定装置１００によって、各地点（ノードｎ１〜ｎ７）を通行する車の量（交通量）、又は各交通路（リンクｆ１〜ｆ６）を通行する車の量（交通量）を、各要素のコストとして推定してもよい。
【００２９】
上述したように、ネットワークを構成する複数の要素のコストは、空間に関するモード及び時間に関するモードを少なくとも含む複数のモードを有するテンソルで示される。本実施形態では、ネットワークを構成する複数の要素のコストが、例えば、観測位置に関する第１モード（空間に関するモード）、観測時間に関する第２モード（時間に関するモード）、及び観測日に関する第３モードを含む３つのモードを有する３階のテンソルで示される。なお、本実施形態では、第３モードは、観測日に関するモードに限定されず、観測時間及び観測位置以外の他の変数に関するモードであればよい。
【００３０】
ネットワークを構成する複数の要素のコストを示す３階のテンソルＸは、例えば、下記式（ｉ）によって表される。なお、下記式（ｉ）及び後述の（ｘiｖ）において、白抜き文字の「Ｒ」は実数を示し、「ｎ１」は、空間における観測数（観測位置の数）を示し、「ｎ２」は、観測時点数（本実施形態では、１日当たりの観測回数）を示し、「ｎ３」は、観測を実施した日数を示している。すなわち、下記式（ｉ）においては、「ｎ１」が第１モードに対応し、「ｎ２」が第２モードに対応し、「ｎ３」が第３モードに対応する。なお、下記式（ｉ）及び後述の（ｘiｖ）中の白抜き文字の「Ｒ」は、後述の他の式中の「Ｒ」とは異なる。
2018112884.tif 000003
【００３１】
上述したように、本実施形態では、ネットワークを構成する複数の要素のコストを、予め設定されたモデルを用いて推定する。モデルは、例えば、上記テンソルを分解することによって得られる。具体的には、例えば、モデルは、上述の３階のテンソルＸがＣＰ分解可能であると仮定して、下記式（ii）によって示される。なお、下記式（ii）において、「λｒ」は係数を示し、「ｕｒ（１）」は上述の第１モードに対応する列ベクトルを示し、「ｕｒ（２）」は上述の第２モードに対応する列ベクトルを示し、「ｕｒ（３）」は上述の第３モードに対応する列ベクトルを示している。また、下記式（ii）において、「Ｒ」は、テンソルのランクを示す。本実施形態では、後述するように、係数λｒに非特許文献１に開示された「multiplicative gamma process」を設定することによって、ランクＲを自動的に決定することができる。
2018112884.tif 000004
【００３２】
なお、上記テンソルで示される複数の要素のうち、データ（コスト）が観測されている要素（欠損していない要素。以下、観測要素とも記載する。）のインデックス“ｉ”を下記式（iii）で示した場合、任意の観測要素の観測データ“ｙｉ”は、下記式（iv）で示すことができる。なお、下記式（iii）において、「Ｉ」は、観測要素のインデックスの集合を示し、「Ｉ」の要素数は、観測要素の総数に等しい。また、下記式（iv）において、「ｘｉ」は任意の観測要素の実際のコストを示し、「εｉ」は測定誤差を示す。
2018112884.tif 000005
2018112884.tif 000006
【００３３】
上記の関係から、尤度関数は、下記式（ｖ）で示すことができる。
2018112884.tif 000007
【００３４】
以下、データ受付部１０及び要素コスト算出部１２について具体的に説明する。なお、以下においては、上記式（ii）〜（ｖ）に基づいて、複数の要素のコストを推定する場合について説明する。
【００３５】
図２を参照して、データ受付部１０の観測データ受付部１４には、複数の観測要素の観測データが入力される。本実施形態では、観測データは、観測位置（ネットワークにおける要素の位置）、観測時間及び観測日に関連付けられている。観測データ受付部１４は、受け付けた観測データを、要素コスト算出部１２の推定部１８へ送信する。なお、上述したように、本実施形態では、データ受付部１０に入力される観測データには、測定誤差が含まれている。
【００３６】
データ受付部１０の事前分布受付部１６には、上記式（ii）〜（ｖ）の各パラメータに対する事前分布のデータが入力される。事前分布受付部１６は、受け付けた事前分布のデータを、要素コスト算出部１２の推定部１８へ送信する。
【００３７】
本実施形態では、各パラメータに対する事前分布は、ネットワークトポロジー及び推定対象となるネットワークの種別等に基づいて適宜設定される。各パラメータの事前分布は、下記式（vi）〜（xiii）によって表される。なお、下記式において、Ｇａ（）は、ガンマ分布を示し、Ｆ（ｋ）は、特徴量を示している。また、以下においては、下記式（xiii）で示される行列Ｕｋを、因子行列と称する。
2018112884.tif 000008
2018112884.tif 000009
2018112884.tif 000010
2018112884.tif 000011
2018112884.tif 000012
2018112884.tif 000013
2018112884.tif 000014
2018112884.tif 000015
但し、上記（viii）式において、
2018112884.tif 000016
とする。
【００３８】
なお、本実施形態では、例えば、因子行列の第１モードに対応する列ベクトルｕｒ（１）の事前分布が、ネットワークトポロジーに基づいて設定されてもよい。以下、第１モードに対応する列ベクトルｕｒ（１）の事前分布をネットワークトポロジーに基づいて設定する場合の一例を簡単に説明する。
【００３９】
図３のネットワーク３０において、リンクｆ１〜ｆ６のフローをそれぞれフローｕ１〜ｕ６とすると、フロー保存側より、下記式（xv），（xvi）が成立する。
ｕ１＝ｕ３＋ｕ４      ・・・（xv）
ｕ２＋ｕ３＝ｕ５＋ｕ６      ・・・（xvi）
【００４０】
上記式（xv），（xvi）の関係から、下記式（xvii）で示す関係が得られる。なお、下記式（xvii）の「σ２」は、非常に小さい値（例えば、１０−５）である。
2018112884.tif 000017
【００４１】
ここで、上記式（xvii）を下記式（xviii）で表す。
Ａｕ＝ε      ・・・（xviii）
【００４２】
上記式（xviii）より、転置行列ＡＴを用いた下記式（xix）が得られる。
ＡＴＡｕ＝ＡＴε      ・・・（xix）
【００４３】
上記式（xix）において、非常に小さい正のλ（例えば、１０−５）より、ｕを、（ＡＴＡ＋λＩ）−１ＡＴεと近似する。このとき、ｕの分散Ｖaｒ（ｕ）は、下記式（xx）で得られる。
Ｖaｒ（ｕ）＝（ＡＴＡ＋λＩ）−１ＡＴＡ（ＡＴＡ＋λＩ）−１・σ２      ・・・（xx）
【００４４】
上記式（xx）で得られる分散Ｖaｒ（ｕ）は、（ＡＴＡ＋λＩ）−１・σ２と近似することができる。以上の考え方に基づいて、第１モードに対応する列ベクトルｕｒ（１）の事前分布は、下記式（xxi）によって設定することができる。
2018112884.tif 000018
【００４５】
以上のように、本実施形態では、ネットワークトポロジーに基づいて決定される行列Ａを用いて、第１モードに対応する列ベクトルｕｒ（１）の事前分布を設定することができる。すなわち、ネットワークトポロジーを考慮して事前分布を設定することができる。なお、上記の事前分布の設定方法は一例であり、上述したように、各パラメータに対する事前分布は、ネットワークトポロジー及び推定対象となるネットワークの種別等に基づいて適宜設定すればよい。
【００４６】
推定部１８の因子行列推定部２４は、予め設定された上述の尤度関数、並びに観測データ受付部１４及び事前分布受付部１６から送信されてきたデータに基づいて、複数の条件付事後分布を生成する。具体的には、因子行列推定部２４は、下記式（xxii）〜（xxxii）で示されるように、上記式（ii）の係数λｒ、係数λｒの事前分布のパラメータδｒ（上記式（vi）参照）、因子行列の列ベクトルｕｒ（ｋ）、列ベクトルｕｒ（ｋ）の平均部分に含まれるβｒ（ｋ）（上記式（viii）参照）、及び上記式（iv）の測定誤差εiの分散を示すパラメータτεについてそれぞれ、条件付事前分布を生成する。なお、本実施形態では、因子行列推定部２４は、後述する時空間相関推定部２６において抽出されたパラメータτｋ，θｋを用いて、条件付事前分布を生成する。
2018112884.tif 000019
2018112884.tif 000020
2018112884.tif 000021
2018112884.tif 000022
2018112884.tif 000023
2018112884.tif 000024
2018112884.tif 000025
2018112884.tif 000026
2018112884.tif 000027
2018112884.tif 000028
2018112884.tif 000029
但し、上記式において、
2018112884.tif 000030
2018112884.tif 000031
2018112884.tif 000032
2018112884.tif 000033
2018112884.tif 000034
2018112884.tif 000035
2018112884.tif 000036
2018112884.tif 000037
である。
【００４７】
更に、因子行列推定部２４は、上記複数の条件付事後分布それぞれに対して複数回のサンプリングを実行することによって、サンプリングごとに、係数λｒの事前分布のパラメータδｒ、測定誤差εiの分散部分のパラメータτε、係数λｒ、因子行列の各列ベクトルｕｒ（ｋ）の平均部分に含まれるβｒ（ｋ）、及び各列ベクトルｕｒ（ｋ）をそれぞれ抽出する。本実施形態では、因子行列推定部２４によって、ギブスサンプリングが実行される。なお、サンプリング回数は、予め設定される。また、上記式（xxxii）において、ｘiは、今までにサンプリングした係数λｒ及び各列ベクトルｕｒ（ｋ）（１≦ｒ≦Ｒ，１≦ｋ≦３）を用いて、上記式(ii)に基づいてテンソルＸを再構成することによって得られる。
【００４８】
推定部１８の時空間相関推定部２６は、予め設定された上述の尤度関数、並びに観測データ受付部１４及び事前分布受付部１６から送信されてきたデータに基づいて、複数の条件付事後分布を生成する。具体的には、時空間相関推定部２６は、下記式（xli）及び（xlii）で示されるように、各モードに対応する列ベクトルｕｒ（ｋ）のパラメータτｋ（上記式（xxvii）等参照）、及び相関行列Σθｋのパラメータθｋ（上記式（xxvii）等参照）についてそれぞれ、条件付事前分布を生成する。
2018112884.tif 000038
2018112884.tif 000039
【００４９】
時空間相関推定部２６は、上記複数の条件付事後分布それぞれに対して複数回のサンプリングを実行することによって、サンプリングごとに、パラメータτｋ，θｋを抽出する。本実施形態では、時空間相関推定部２６によって、ギブスサンプリングが実行される。なお、サンプリング回数は、予め設定される。時空間相関推定部２６によって抽出されたパラメータτｋ，θｋは、上述のように、因子行列推定部２４によって用いられる。なお、パラメータτｋ，θｋは、ネットワークトポロジー及び推定対象となるネットワークの種別によっては、ギブスサンプリングするパラメータとして設定する必要はない。例えば、上記式（xxi）では、θ１は設定されていない。
【００５０】
推定値算出部２０は、因子行列推定部２４によってサンプリングされた係数λｒ、及び各列ベクトルｕｒ（ｋ）を用いて、サンプリングごとに、上記式（ii）に基づいて、テンソルＸを再構成する。すなわち、因子行列推定部２４によるサンプリングごとに、推定値算出部２０によって、複数の要素のコストの推定値が算出される。本実施形態では、推定値算出部２０は、要素ごとに複数回のサンプリングに基づいて算出された複数の推定値の平均値を算出する。推定値算出部２０は、算出した上記平均値を、各要素のコストとしてコスト表示部２２へ出力する。コスト表示部２２は、例えば、液晶表示パネル又は有機ＥＬパネル等の表示パネルを含み、推定値算出部２０から出力されたコストを表示する。
【００５１】
このようにして、本実施形態では、推定値算出部２０によって、ネットワークの複数の要素のうち、実際にデータ観測が行なわれていない要素のデータ（コスト）が補完される。すなわち、欠損データが補完される。また、本実施形態では、測定誤差を考慮した観測データに基づいて、各要素のコストが算出される。これにより、各要素の観測データのデノイジングを行なうことができる。
【００５２】
なお、ネットワークが配水網である場合、推定値算出部２０において、管路抵抗を各要素のコストとして算出してもよい。この場合、推定値算出部２０は、算出した上記平均値（各管路の流量）に基づいて、各管路の管路抵抗を算出する。管路抵抗は、例えば、管路に水を供給するためのポンプの１日の消費電力（観測値）と、管路を流れる水の１日の流量とに基づいて求めることができる。
【００５３】
具体的には、例えば、「１日における各ポンプの消費電力（観測値）と理論消費電力との差ｙ」、「推定値算出部２０によって算出された各管路の流量Ｘ（ただし、定数項を導入するため、全ての要素が１である列ベクトルを含む）」、及び「各管路の管路抵抗β」の関係は、目的変数ｙ、説明変数Ｘ、及び回帰係数βとして、下記式（xliii）に示す線形回帰モデルで示すことができる。
2018112884.tif 000040
【００５４】
本実施形態では、例えば、推定値算出部２０において、上記の回帰係数βをベイズ推定する。回帰係数βの事前分布が下記式（xliv）で表されるとすると、事後分布は、下記式（xlv）で表される。このように、推定値算出部２０において回帰係数βのベイズ推定を行なうことによって、各管路の管路抵抗を、各要素のコストとして算出することができる。
2018112884.tif 000041
2018112884.tif 000042
【００５５】
なお、上述の実施形態では、ネットワークを構成する複数の要素が３階のテンソルで示される場合について説明した。しかしながら、ネットワークを構成する複数の要素は、２階のテンソルで示されてもよく、４階以上のテンソルで示されてもよい。この場合、テンソルに応じたモデルが用いられる。
【００５６】
［装置動作］
次に、本発明の実施の形態における推定装置の動作について図４を用いて説明する。図４は、本発明の実施の形態における推定装置の動作を示すフロー図である。以下の説明においては、適宜図１〜図３を参酌する。また、本実施の形態では、推定装置１００を動作させることによって、推定方法が実施される。よって、本実施の形態における推定方法の説明は、以下の推定装置１００の動作説明に代える。
【００５７】
図４に示すように、捜索装置１００においては、まず、データ受付部１０が、複数の観測要素の観測データ及び事前分布のデータを取得する（ステップＡ１）。本実施形態では、上述したように、要素コスト算出部１２によって観測データが取得され、観測データ受付部１４によって、上記式（ii）〜（ｖ）の各パラメータに対する事前分布のデータが取得される。
【００５８】
次に、推定部１８が、複数の条件付事前分布を生成する（ステップＡ２）。本実施形態では、上述したように、因子行列推定部２４によって、上記式（ii）の係数λｒ、係数λｒの事前分布のパラメータδｒ（上記式（vi）参照）、因子行列の列ベクトルｕｒ（ｋ）、列ベクトルｕｒ（ｋ）の平均部分に含まれるβｒ（ｋ）（上記式（viii）参照）、及び上記式（iv）の測定誤差εiの分散を示すパラメータτεについてそれぞれ、条件付事前分布が生成される。また、上述したように、時空間相関推定部２６によって、各モードに対応する列ベクトルｕｒ（ｋ）のパラメータτｋ（上記式（xxvii）参照）、及び相関行列Σθｋのパラメータθｋ（上記式（xxvii）参照）についてそれぞれ、条件付事前分布が生成される。
【００５９】
次に、推定部１８が、ステップＡ２において生成された条件付事前分布から、各パラメータを抽出する（ステップＡ３）。本実施形態では、上述したように、因子行列推定部２４によって複数回のサンプリングが実行され、サンプリングごとに、係数λｒの事前分布のパラメータδｒ、測定誤差εiの分散部分のパラメータτε、係数λｒ、因子行列の各列ベクトルｕｒ（ｋ）の平均部分に含まれるβｒ（ｋ）、及び各列ベクトルｕｒ（ｋ）がそれぞれ抽出される。また、上述したように、時空間相関推定部２６によって複数回のサンプリングが実行され、サンプリングごとに、パラメータτｋ，θｋが抽出される。
【００６０】
最後に、推定値算出部２０が、因子行列推定部２４によってサンプリングされた係数λｒ、及び各列ベクトルｕｒ（ｋ）を用いて、複数の要素のコストの推定値を算出する（ステップＡ４）。本実施形態では、上述したように、推定値算出部２０は、因子行列推定部２４によるサンプリングごとに上記式（ii）に基づいてテンソルＸを再構成し、要素ごとに算出された複数の推定値の平均値を算出する。さらに、推定値算出部２０は、算出した上記平均値を、各要素のコストとしてコスト表示部２２へ出力する。
【００６１】
［プログラム］
本発明の実施の形態におけるプログラムは、コンピュータに、図４に示すステップＡ１〜Ａ４を実行させるプログラムであればよい。このプログラムをコンピュータにインストールし、実行することによって、本実施の形態における推定装置と推定方法とを実現することができる。この場合、推定装置となるコンピュータのＣＰＵ（Central Processing Unit）は、観測データ受付部１４、事前分布受付部１６、推定値算出部２０、コスト表示部２２、因子行列推定部２４及び時空間相関推定部２６として機能し、処理を行なう。
【００６２】
また、本実施の形態におけるプログラムは、複数のコンピュータによって構築されたコンピュータシステムによって実行されてもよい。この場合は、例えば、各コンピュータが、それぞれ、観測データ受付部１４、事前分布受付部１６、推定値算出部２０、コスト表示部２２、因子行列推定部２４及び時空間相関推定部２６のいずれかとして機能してもよい。
【００６３】
［物理構成］
ここで、本実施の形態におけるプログラムを実行することによって、推定装置を実現するコンピュータについて図を用いて説明する。図５は、本発明の実施の形態における推定装置を実現するコンピュータの一例を示すブロック図である。
【００６４】
図５に示すように、コンピュータ１１０は、ＣＰＵ１１１と、メインメモリ１１２と、記憶装置１１３と、入力インターフェイス１１４と、表示コントローラ１１５と、データリーダ／ライタ１１６と、通信インターフェイス１１７とを備える。これらの各部は、バス１２１を介して、互いにデータ通信可能に接続される。
【００６５】
ＣＰＵ１１１は、記憶装置１１３に格納された、本実施の形態におけるプログラム（コード）をメインメモリ１１２に展開し、これらを所定順序で実行することにより、各種の演算を実施する。メインメモリ１１２は、典型的には、ＤＲＡＭ（Dynamic Random Access Memory）等の揮発性の記憶装置である。また、本実施の形態におけるプログラムは、コンピュータ読み取り可能な記録媒体１２０に格納された状態で提供される。なお、本実施の形態におけるプログラムは、通信インターフェイス１１７を介して接続されたインターネット上で流通するものであってもよい。
【００６６】
また、記憶装置１１３の具体例としては、ハードディスクドライブの他、フラッシュメモリ等の半導体記憶装置が挙げられる。入力インターフェイス１１４は、ＣＰＵ１１１と、キーボードおよびマウスといった入力機器１１８との間のデータ伝送を仲介する。表示コントローラ１１５は、ディスプレイ装置１１９と接続され、ディスプレイ装置１１９での表示を制御する。
【００６７】
データリーダ／ライタ１１６は、ＣＰＵ１１１と記録媒体１２０との間のデータ伝送を仲介し、記録媒体１２０からのプログラムの読み出し、およびコンピュータ１１０における処理結果の記録媒体１２０への書き込みを実行する。通信インターフェイス１１７は、ＣＰＵ１１１と、他のコンピュータとの間のデータ伝送を仲介する。
【００６８】
また、記録媒体１２０の具体例としては、ＣＦ（Compact Flash（登録商標））およびＳＤ（Secure Digital）等の汎用的な半導体記憶デバイス、フレキシブルディスク（Flexible Disk）等の磁気記憶媒体、またはＣＤ−ＲＯＭ（Compact Disk Read Only Memory）などの光学記憶媒体が挙げられる。
【００６９】
なお、本実施の形態における推定装置は、プログラムがインストールされたコンピュータではなく、各部に対応したハードウェアを用いることによっても実現可能である。更に、推定装置は、一部がプログラムで実現され、残りの部分がハードウェアで実現されていてもよい。
【産業上の利用可能性】
【００７０】
以上のように、本発明によれば、ネットワークを構成する複数の要素のうち未観測の要素のコストを適切かつ効率よく推定することができる。
【符号の説明】
【００７１】
１０      データ受付部
１２      要素コスト算出部
１４      観測データ受付部
１６      事前分布受付部
１８      推定部
２０      推定値算出部
２２      コスト表示部
２４      因子行列推定部
２６      時空間相関推定部
３０      ネットワーク
１００    推定装置
１１０    コンピュータ
１１１    ＣＰＵ
１１２    メインメモリ
１１３    記憶装置
１１４    入力インターフェイス
１１５    表示コントローラ
１１６    データリーダ／ライタ
１１７    通信インターフェイス
１１８    入力機器
１１９    ディスプレイ装置
１２０    記録媒体
１２１    バス
【図１】
2018112884.tif 000043
【図２】
2018112884.tif 000044
【図３】
2018112884.tif 000045
【図４】
2018112884.tif 000046
【図５】
2018112884.tif 000047