(57)【特許請求の範囲】
【請求項１】
少なくともエージェントの初期位置の情報に基づいて、対象エリアにおけるエージェントの分布の情報を含むシミュレーション結果を出力するシミュレーション手段と、
前記シミュレーション結果におけるエージェントの分布を観測結果に基づくエージェントの分布に近づけるように再配置することであって、前記再配置による移動距離の合計値を最小化するという制約下で再配置することによってデータ同化を行う再配置手段と、
前記再配置後のエージェントの分布に基づいて推定された前記対象エリアにおけるエージェントの分布の情報を出力する出力手段と
を備える推定装置。
【請求項２】
前記出力手段は、パーティクルフィルタによりフィルタリングされた前記再配置後のエージェントの分布に基づいて推定された、前記対象エリアにおけるエージェントの分布の情報を出力する、請求項１に記載の推定装置。
【請求項３】
前記移動距離は、前記シミュレーション結果におけるエージェントの分布と、前記再配置後のエージェントの分布との間のＥＭＤ（Earth Mover's Distance）に基づく距離である、請求項１に記載の推定装置。
【請求項４】
予め設定された確率分布に基づいて、エージェントの行動の種別を示す行動パラメータをエージェントの行動の推定値として出力する行動推定手段をさらに備え、
前記出力手段は、前記行動推定手段により出力された行動パラメータをそれぞれのエージェントに対応付けて、前記推定されたエージェントの分布の情報を出力する、請求項１から３のいずれか一項に記載の推定装置。
【請求項５】
制御部を備える情報処理装置により実施される方法であって、
前記制御部が、少なくともエージェントの初期位置の情報に基づいて、対象エリアにおけるエージェントの分布の情報を含むシミュレーション結果を出力すること、
前記制御部が、前記シミュレーション結果におけるエージェントの分布を観測結果に基づくエージェントの分布に近づけるように再配置することであって、前記再配置による移動距離の合計値を最小化するという制約下で再配置することによってデータ同化を行うこと、
前記制御部が、前記再配置後のエージェントの分布に基づいて推定された前記対象エリアにおけるエージェントの分布の情報を出力することと
を備える推定方法。
【請求項６】
コンピュータを、
少なくともエージェントの初期位置の情報に基づいて、対象エリアにおけるエージェントの分布の情報を含むシミュレーション結果を出力するシミュレーション手段、
前記シミュレーション結果におけるエージェントの分布を観測結果に基づくエージェントの分布に近づけるように再配置することであって、前記再配置による移動距離の合計値を最小化するという制約下で再配置することによってデータ同化を行う再配置手段、
前記再配置後のエージェントの分布に基づいて推定された前記対象エリアにおけるエージェントの分布の情報を出力する出力手段、
として機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、推定技術に関する。
【背景技術】
【０００２】
従来、様々なシミュレーション技術が知られている。例えば、特許文献１には、低次元シミュレーションについての技術が開示されている。特許文献１には、低次元シミュレーションを高精度で行うための熱流体シミュレーション装置について開示されている。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開２０１４−２６４４０号公報
【発明の概要】
【発明が解決しようとする課題】
【０００４】
しかしながら、上記の従来技術では、シミュレーションの精度には限界があり、より高い精度で行うことが可能な技術が望まれている。
【０００５】
本発明は、上記に鑑みてなされたものである。本発明の目的は、より高い精度でシミュレーションを行うことが可能な技術を提供することにある。
【課題を解決するための手段】
【０００６】
本発明における推定装置は、少なくともエージェントの初期位置の情報に基づいて、対象エリアにおけるエージェントの分布の情報を含むシミュレーション結果を出力するシミュレーション手段と、前記シミュレーション結果におけるエージェントの分布を観測結果に基づくエージェントの分布に近づけるように再配置することであって、前記再配置による移動距離の合計値を最小化するという制約下で再配置することによってデータ同化を行う再配置手段と、前記再配置後のエージェントの分布に基づいて推定された前記対象エリアにおけるエージェントの分布の情報を出力する出力手段とを備える。
【０００７】
本発明における推定方法は、制御部を備える情報処理装置により実施される方法であって、前記制御部が、少なくともエージェントの初期位置の情報に基づいて、対象エリアにおけるエージェントの分布の情報を含むシミュレーション結果を出力すること、前記制御部が、前記シミュレーション結果におけるエージェントの分布を観測結果に基づくエージェントの分布に近づけるように再配置することであって、前記再配置による移動距離の合計値を最小化するという制約下で再配置することによってデータ同化を行うこと、前記制御部が、前記再配置後のエージェントの分布に基づいて推定された前記対象エリアにおけるエージェントの分布の情報を出力することとを備える。
【０００８】
本発明におけるプログラムは、コンピュータを、少なくともエージェントの初期位置の情報に基づいて、対象エリアにおけるエージェントの分布の情報を含むシミュレーション結果を出力するシミュレーション手段、前記シミュレーション結果におけるエージェントの分布を観測結果に基づくエージェントの分布に近づけるように再配置することであって、前記再配置による移動距離の合計値を最小化するという制約下で再配置することによってデータ同化を行う再配置手段、前記再配置後のエージェントの分布に基づいて推定された前記対象エリアにおけるエージェントの分布の情報を出力する出力手段、として機能させる。
【０００９】
また、本発明のプログラムは、ＣＤ−ＲＯＭ等の光学ディスク、磁気ディスク、半導体メモリなどの各種の記録媒体を通じて、又は通信ネットワークなどを介してダウンロードすることにより、コンピュータにインストール又はロードすることができる。
【発明の効果】
【００１０】
本発明によれば、より高い精度でシミュレーションを行うことが可能な技術を提供することができる。
【図面の簡単な説明】
【００１１】
【図１】一実施形態における推定装置の構成を示す概念図である。
【図２】一実施形態における推定装置の処理構成の一部を示す概念図である。
【図３】一実施形態における推定装置による処理の例を示すフローチャートである。
【図４】一実施形態における推定装置による処理に用いられる係数の設定の例を示す図である。
【発明を実施するための形態】
【００１２】
以下、本発明の実施形態について、図面を参照しつつ詳細に説明する。ただし、発明の範囲をこれらに限定するものではない。
【００１３】
また、本明細書等において、「部」とは、単に物理的構成を意味するものではなく、その構成が有する機能をソフトウェアによって実現する場合も含む。また、１つの構成が有する機能が２つ以上の物理的構成により実現されても、２つ以上の構成の機能が１つの物理的構成により実現されてもよい。
【００１４】
［前提］
本実施形態では、推定装置による人、自動車、電車などを含む移動体（移動可能なオブジェクト）の分布の推定方法の例について説明する。まず、本実施形態における移動体の分布の推定方法の前提について説明する。なお、本実施形態において、「エージェント」という用語を使用しているが、説明の便宜上、「エージェント」は、コンピュータ上でシミュレートされる対象としてのエージェントや、実世界における人、自動車、電車などの移動体を含む概念であることとして使用されている。
【００１５】
近年、ＧＰＳ（Global Positioning System）などによりユーザの位置データをリアルタイムで取得可能である。具体的な分析方法として、例えば、１辺が２５０ｍの矩形領域に滞在するエージェントの数のデータを１時間ごとに記録して分析する方法が知られており、このように記録されたデータを集中分布（aggregated-distribution）データと呼ぶ。
【００１６】
リアルタイム観察によって取得されるデータが集中分布データである場合、膨大なエージェントの数及び広範囲なエリアについて時系列にシミュレーションを行うことが必要となるエージェントの分布の推定は、高次元な多変数時系列予測問題として定式化することができる。
【００１７】
一般に、状態空間モデルのパラメータを推定するための非線形的な方法は小次元の問題に対して特に有効である。例えば、高次元の問題であるエージェントの分布の推定するために、非線形という観点で、シミュレーションに基づく推定法であるパーティクルフィルタが用いられることが多くあるが、単純にパーティクルフィルタのみを用いて推定を行った場合、高い推定精度が得られない場合がある。
【００１８】
上記のような課題を解決するために、本願の発明者らは、パーティクルフィルタにおいて、シミュレーション結果のエージェントの位置と、再配置後のエージェントの位置との間の距離を最小化するという制約下で、尤度の高い位置にエージェントを再配置する処理を行い、当該処理結果を用いてエージェントの分布の推定を行うことにより、高い精度の推定結果が得られることを見出した。以下に説明する実施形態は、このような考えに基づくものである。
【００１９】
次に、本実施形態において用いられる各種の条件、変数及び数式等について定義する。
ユニバーサルエリアＵ１は、エージェントが移動可能な地理的空間全体である。対象エリアＵは、ユニバーサルエリアＵ１のサブセットであり、エージェントの行動が推定される地理的空間である。ユニットエリアｕiは、対象エリアＵのパーティションＰ１の要素である。つまり、次式が成り立つ。
0006632054.tif 000002
【００２０】
0006632054.tif 000003
【００２１】
ユニットエリアは任意の方法で区分されたエリアであるが、ユニットエリアの基本的な例は、地理的グリッドシステムにおけるセルである。以下に地理的グリッドシステム及びセルについて説明する。
【００２２】
地理的グリッドシステムは、おおよそ直交するグリット線により地理的空間をおおよそ矩形形状の複数のセルに区切る。対象セルＣは、対象エリアＵとの間で空ではない共通集合を有する。すなわち、次式が成り立つ。
0006632054.tif 000004
【００２３】
0006632054.tif 000005
【００２４】
本実施形態における状態空間モデルにおいて、時間ｔにおける状態変数ｘtは、エージェントの位置と、後述する行動モデルのパラメータとからなる。状態変数ｘtは、状態変数の時間発展（時間の経過に伴って変化した状態変数）を提供するシステムモデルｆtにより求めることができ、次式で示される。
ｘt＝ｆt（ｘt-1，ｖt）
ここで、ｖtはシステムノイズである。システムノイズｖtを例えば、ガウスノイズとすることができる。
【００２５】
また、時間ｔにおける観測データｙtは、観測データを提供するプロセスを表す観測モデルｈtにより求めることができ、次式で示される。
ｙt＝ｈt（ｘt，ｗt）
ここで、ｗtは観測ノイズである。観測ノイズｗtは例えば、ＧＰＳ等による測位に基づく観測誤差を考慮して設定される。
【００２６】
マルチエージェントシミュレーションによりエージェントの分布を推定するために、状態変数ｘtは、要素としてエージェントｌの位置ａt含む。ａtは、次式により示される。
0006632054.tif 000006
【００２７】
また、行動モデルのパラメータのリストθtは、次式により示される。
0006632054.tif 000007
【００２８】
0006632054.tif 000008
【００２９】
上記を考慮して、状態空間モデルは、次のように記述することができる（ここではシステムノイズを無視している。）。
0006632054.tif 000009
【００３０】
0006632054.tif 000010
【００３１】
式（１）を用いたエージェントの分布推定問題は、モデリング段階と、推定段階とで次のように定式化される。
【００３２】
0006632054.tif 000011
【００３３】
［構成］
本実施形態における推定装置の構成について説明する。
図１を参照して、本実施形態に係る推定装置のハードウェアの概略構成の例を説明する。同図に示すように、推定装置１０は、制御部１１、通信部１４、記憶部１５及び入出力部１６を主に備える。制御部１１は、ＣＰＵ（Ｃｅｎｔｒａｌ  Ｐｒｏｃｅｓｓｉｎｇ  Ｕｎｉｔ）１２及びメモリ１３を主に備えて構成される。推定装置１０は、専用若しくは汎用のサーバ・コンピュータ又はパーソナルコンピュータなどの情報処理装置を用いて実現することができる。また、推定装置１０は、例えば、ＣＰＵ１２がメモリ１３等に格納された所定のプログラムを実行することにより、各種の機能実現手段として機能する。なお、推定装置１０は、単一の情報処理装置より構成されるものであっても、ネットワーク上に分散した複数の情報処理装置より構成されるものであってもよい。
【００３４】
制御部１１では、ＣＰＵ１２は、記憶部１５等に記憶されたプログラムをメモリ１３に展開して実行することにより、推定装置１０が備える各種構成の動作を制御し、また、各種処理の実行を制御する。制御部１１において実行される処理の詳細は後述する。
【００３５】
通信部１４は、ネットワークを介して外部の各種情報処理装置との間で通信をするための通信インタフェースである。通信部１４は、例えば、推定装置１０における処理に必要な情報を外部装置から受信し、また、推定装置１０における処理結果の情報を外部装置に送信することができる。
【００３６】
記憶部１５は、ハードディスク等の記憶装置によって構成される。記憶部１５は、制御部１１における処理の実行に必要な各種プログラムや、制御部１１による処理結果の情報など、各種の情報を記憶する。
【００３７】
入出力部１６は、各種の情報の入出力を行うための処理部である。入出力部１６には、例えば、操作部及び表示部などのユーザインタフェースとしての構成や、外付けの記憶装置を接続するためのインタフェースなどが含まれる。
【００３８】
［機能］
本実施形態における推定装置１０の機能について説明する。推定装置１０は、対象エリアにおけるエージェントの分布を推定する機能を有する。
【００３９】
１．処理の概要
まず、推定装置１０において、パーティクルフィルタによるエージェントの分布の推定のために実行される処理のアルゴリズムを擬似的なプログラム言語で以下に示す。なお、このアルゴリズムは本実施形態における推定処理の概念を示すための例に過ぎず、当該推定処理の概念に矛盾しない範囲でこのアルゴリズムにおける各ステップの処理順序や処理方法を変更することができる。
0006632054.tif 000012
【００４０】
0006632054.tif 000013
【００４１】
アルゴリズムの行５から始まるｗｈｉｌｅ制御文内では、推定する期間内における各時間における推定処理が行われる。行７から始まるｆｏｒ制御文内では、各パーティクルに対応するシミュレータに存在するエージェントついて、推定処理を行う。ｆｏｒ制御文内における推定処理において、まず、後述する行動モデルに基づいて、時間ｔにおける当該シミュレータに存在する全てのエージェントの行動を決定し、当該行動及び時間ｔ−１におけるエージェントの位置に基づいて時間ｔにおける複数のエージェントの位置をシミュレーションする。なお、各パーティクルにおいて、行動モデルに基づいて決定される行動は変更されうる。さらに、時間ｔにおける観測データを取得し、時間ｔ−１における行動パラメータに基づいて時間ｔにおける行動パラメータを決定する。その後、後述するＥＭＤ部により、シミュレーションの結果得られた複数のエージェントの位置を再配置するときの移動先の座標を計算し、当該座標に基づいて、再配置後のエージェントの位置を決定する。当該パーティクルの重みが観測データに基づいて決定され、その後、次のパーティクルの処理のために、パーティクルがリサンプリングされる。ｆｏｒ制御文内における各処理の詳細は後述する。
【００４２】
パーティクルの状態は、モデリング段階において、行動モデル、マルチエージェントシミュレーション、観測モデルが決定されているという条件下で推定される。パーティクルの状態ベクトル（エージェントの位置及び行動パラメータを含む。）の時間発展は、マルチエージェントシミュレーションと、行動パラメータに関する差分方程式により与えられる。各エージェントの目的地は、行動モデルにより決定される。観測時に、パーティクルの状態ベクトルは、観測データを用いて更新され（更新方法の詳細は後述する。）、観測データを用いてパーティクルの重み付けがなされ、これらの結果に基づいて、推定結果が算出される。なお、上記のアルゴリズムに示されたｅｑ．（３）及びｅｑ．（４）は、後述する式（３）及び式（４）にそれぞれ対応する。
【００４３】
２．処理の詳細
図２は、上記のアルゴリズムにおけるｆｏｒ制御文内の処理を実行するために主に必要な処理構成の例を示している。これらの処理構成は、例えば、制御部１１において、ＣＰＵ１２が、記憶部１５等に記憶されたプログラムをメモリ１３に展開して実行することにより実現される。図３を参照しながら、推定処理の詳細を説明する。
【００４４】
２．１．ＢＭ
ＢＭは、予め選択された行動モデルと行動パラメータに基づいて、パーティクルｉについて、時間ｔにおけるエージェントの行動を決定する。例えばＢＭは、エージェントの行動として、エージェントの時間ｔにおける目的地を決定する。ＢＭは、行動を決めるために、任意の行動モデルを適用することができる。また、ＢＭは、シミュレータにエージェントの初期位置及び目的地の情報を与える。
【００４５】
図３を参照して、ＢＭが、予め選択された行動モデルに基づいてエージェントの目的地を決定する方法の例を説明する。この例における行動モデルでは、行動パラメータに基づいて、家、駅、店及び現在地のいずれかが目的地としてある確率で選択される。図３においては、行動パラメータに含まれる要素が、ｐ１からｐ４で示されている。ｐ１からｐ４にはそれぞれ、０から１の間の数値が予め設定されている。また、図３に示す処理は、制御部１１が記憶部１５に記憶されたプログラムをメモリ１３に展開して実行することによって制御される。
【００４６】
まず、ステップＳ３１でエージェントの現在地が駅であるか否かが判断され、駅である場合、処理はステップＳ３２へ進み、駅でない場合、処理はステップＳ３３へ進む。
【００４７】
ステップＳ３２において、０から１の間の乱数を出力する乱数生成器ｒ（１）の出力がｐ４未満であるか否かが判断される。ｐ４未満である場合、現在地が目的地であると判断される、他の場合、目的地は家であると判断される。
【００４８】
ステップＳ３３において、乱数生成器ｒ（１）は乱数を出力し、当該出力がｐ１未満であるか否かが判断される。ｐ１未満である場合、目的地は家であると判断され、他の場合、処理はステップＳ３４へ進む。
【００４９】
ステップＳ３４において、乱数生成器ｒ（１）は乱数を出力し、当該出力がｐ２未満であるか否かが判断される。ｐ２未満である場合、目的地は駅であると判断され、他の場合、処理はステップＳ３５へ進む。
【００５０】
ステップＳ３５において、乱数生成器ｒ（１）は乱数を出力し、当該出力がｐ３未満であるか否かが判断される。ｐ３未満である場合、目的地は店であると判断され、他の場合、現在地が目的地であると判断される。
【００５１】
なお、目的地を判断するための閾値として上記では行動パラメータをそのまま用いたが、エリア別、時間帯別、曜日別、通常時／災害時別、及び年齢別など、の各種条件と行動パラメータの関数として閾値を設定することもできる。また、上記の条件を組み合わせて、確率を設定することができる。例えば、家を目的地とするかの判断時に用いられる値として、年齢と現在地から自宅までの距離とに応じて設定される属性別帰宅確率係数ｅを定義することができる。このとき、例えば、図３のステップＳ３３において、ｐ１に属性別帰宅確率係数ｅを乗じた値ｅｐ１が乱数生成器ｒ（１）の出力未満であるかを判断し、ｅｐ１未満である場合、目的地は家であると判断され、他の場合、処理はステップＳ３４に進むようにすることができる。
【００５２】
図４には、年齢と現在地から自宅までの距離とに応じて設定される属性別帰宅確率係数ｅの設定方法の例がグラフで示されている。この例では、「２０才未満」、「２０才代及び４０才代」、及び「３０才代及び５０才以上」という年齢の範囲別に、現在地から自宅までの距離に応じて設定される属性別帰宅確率係数ｅが示されている。
【００５３】
２．２．ＳＩＭ
図２におけるＳＩＭは、エージェントの初期位置及び目的地等の情報を入力として、複数のエージェントの動きをシミュレーションすることにより、時間ｔにおける複数のエージェントの位置ａtから、時間ｔ＋１におけるエージェントの位置ａt+1を特定し、設定する。推定装置１０上でＳＩＭによりエージェントの動きをシミュレーションすることにより、実世界における人、車、電車などを含む移動体としてのエージェントの動き及び分布をシミュレーションすることができる。なお、ＳＩＭには、任意のシミュレータを適用することができる。
【００５４】
本実施形態では、広域道路網交通流シミュレーションシステム「SOUND」（http://www.i-transportlab.jp/products/sound/）にルート検索機能を追加したシミュレータを用いる。エージェントが動くルートは、混雑を加味した通過時間を重みとして用いてダイクストラで決定される。ルートの通過時間は、混雑による待ち時間と走行時間の和とする。シミュレーションの実行には、これらのパラメータに加えて、道路ネットワーク、エージェントの初期位置（例えば、シミュレーション開始時の位置。または、時間ｔ＋１のときの位置を推定する場合における時間ｔのときの位置）、出発地（トリップの開始位置）及び目的地等の情報のセットが与えられている必要がる。あるエージェントの初期位置、出発地及び目的地等の情報は、例えば、人の流れデータからサンプリングされたデータに基づいて作成される。目的地の作成方法については、前述のＢＭの説明において説明したとおりである。
【００５５】
ここで、人の流れデータとは、何らかのデータソースをもとに粒度を揃えた時空間位置で、人々の流動をデータ化したものである。データソースとして、例えば、パーソントリップ調査の結果（パーソントリップデータ）が用いられる。パーソントリップ調査とは、対象の地域における人の動きを調べ、交通機関の実態を把握する調査である。
【００５６】
２．３．ＥＭＤ部
図２におけるＥＭＤ（以下、「ＥＭＤ部」という。）は、ＳＩＭにより特定されたエージェントの位置ａt+1を、時間ｔ＋１におけるエージェントの位置の観測データｙt+1に近づけるように再配置し、再配置後の位置ａt+1を出力する。すなわち、ＥＭＤ部は、ＳＩＭによるシミュレーション結果におけるエージェントの分布を観測結果に基づくエージェントの分布に近づけるように再配置することによってデータ同化を行う再配置手段として機能する。また、当該再配置は、再配置による移動距離の合計値を最小化するという制約下で行われる。以下に、このような再配置の処理の意義と、処理の詳細を説明する。
【００５７】
通常のパーティクルフィルタの処理においては、シミュレーションの結果、尤度の低いパーティクルは破棄されてしまうため、尤度の低いパーティクルが多い場合、推定結果の精度が下がる可能性がある。特に、エージェントの分布の推定のような高次元の問題の場合、全てのパーティクルの状態の尤度が低くなる可能性が高い。従って、本実施形態においては、ＥＭＤ部が、観測データを用いて尤度が高くなるようにパーティクルに対応するエージェントの位置を再配置することにより、破棄されてしまうパーティクルの数を減らすことができ、さらに破棄されないパーティクルの尤度が高くなり、その結果、推定結果の精度を高くすることができる。なお、観測データに近づける再配置の方法は無限に存在する。従って、本実施形態では、後述するように、シミュレーション結果を観測値にどれだけ近づけるかを決めるパラメータＫを用いて、再配置後の各ユニットエリアのエージェントの数を決める。さらに、その決められたエージェント数になるような無数の再配置の方法から、移動距離が最小となる再配置の方法を選び、当該選ばれた方法によりエージェントの再配置を行う。
【００５８】
すなわち、本実施形態におけるＥＭＤ部による再配置の方法の基本的なコンセプトとして、シミュレーション結果を観測値にどれだけ近づけるかを予め設定し、さらに、再配置によるエージェントの移動量の合計距離ができるだけ小さくなるように、エージェントの分布を観測データに近づける。
【００５９】
0006632054.tif 000014
【００６０】
0006632054.tif 000015
【００６１】
0006632054.tif 000016
【００６２】
0006632054.tif 000017
ここで、パラメータＫは０から１の間の任意の実数である。
【００６３】
0006632054.tif 000018
【００６４】
式（３）において、再配置におけるエージェントの移動距離が最小になるという制約下で、各パーティクルにおける複数のエージェントの位置が観測データに近づけられる。本実施形態では、式（３）をの時間発展として採用する。に関する分布は、ガウス分
0006632054.tif 000019
【００６５】
２．４．ＥＭＤ（Earth Mover's Distance）による状態ベクトルの計算
上記のように、ＥＭＤ部は、パーティクルの時間発展を求めるために式（３）を用いるが、式（３）を用いるためには、式（２）最適化問題を解く方法を決める必要がある。ここで、以下に説明するように、本願の発明者らは、ＥＭＤ（Earth Mover's Distance）を求める最適化問題と、式（２）の最適化問題が同等であることを見出した。従って、式（
0006632054.tif 000020
定し、ａからａ’に移動するエージェントを荷物量（ＥＭＤの文脈におけるフロー）として設定し、ユニットエリア間の距離を分布間の距離として設定することにより、ＥＭＤの最適化問題を解くことができる。このようにして解かれたＥＭＤを与えるフローについて
0006632054.tif 000021
般に、ＥＭＤの最適化問題を解くための解法は、多くのものが知られているが、本実施形態において、任意の解法を採用することができる。
【００６６】
ＥＭＤを求める最適化問題と、式（２）の最適化問題が同等であることについて説明する。ＥＭＤは、２つの分布の差異を測る尺度と解釈することができる。ＥＭＤの文脈では、分布をシグネチャと呼ばれる形式で表現し、次元数がＮのシグネチャは、
0006632054.tif 000022
と定義される。ここで、ｍjはｊ番目の要素の位置、ｗjはその位置の重みである。シグネチャＰとＱの間のＥＭＤは、次に示すような制約された最適化問題の解を用いて定義される。
0006632054.tif 000023
【００６７】
フローは、位置ｊ１から位置ｊ２へ移動する荷物量の合計を表す。コストは、位置ｊ１から位置ｊ２の距離を表す。ここで、式（５）のＥＭＤを求める式の分母が定数であることに留意し、ＥＭＤを与えるフローを求める式に書き換えると次のとおりとなる。
0006632054.tif 000024
【００６８】
0006632054.tif 000025
【００６９】
0006632054.tif 000026
なお、観測データｙtは、様々な方法で取得することができ、例えば、携帯端末の測位データを用いて取得することができる。このとき、観測データｙtは、測位データ、測位誤差、携帯電話の普及率、測位データ利用許可率などの値を用いて算出される。
【００７０】
0006632054.tif 000027
ここで、Ａは正規化係数である。ＥＭＤ部により出力された再配置後の複数のエージェントの位置（すなわち、エージェントの分布）は、ＷＥＩにより算出された重みに基づいて、後の処理でフィルタリングされ、当該フィルタリングの結果に基づいてエージェントの分布の推定結果の情報が出力される。
【００７１】
0006632054.tif 000028
、新たなパーティクルにおいてエージェントのシミュレーション等の処理が行われる。全てのパーティクルの処理が終わった後、時刻ｔが観測値が得られる時刻である場合はリサンプリングの処理を行い、各ユニットエリア毎の各時間におけるエージェントの分布の推定値が算出される。入出力部１６は、当該算出されたエージェントの分布の情報を推定結果として出力する出力手段として機能する。また、入出力部１６は、推定結果として出力される分布に含まれるエージェントのそれぞれに対して、ｂにより設定された行動パラメータを対応付けて出力することができる。
【００７２】
以上のように、本実施形態によれば、パーティクルフィルタにおいて、シミュレーション結果におけるエージェントの位置と、再配置後のエージェントの位置との間の距離を最小化するという制約下で、エージェントの位置を観測データに近づけることにより、尤度の高い位置にエージェントを再配置する処理を行う。その後、当該処理結果を用いてエージェントの分布の推定を行う。その結果、高い精度でエージェントの分布の推定を行うことができる。
【００７３】
［変形例］
本発明は、上記した実施形態に限定されるものではなく、本発明の要旨を逸脱しない範囲内において、他の様々な形で実施することができる。上記実施形態はあらゆる点で単なる例示にすぎず、限定的に解釈されるものではない。
【００７４】
例えば、上記の実施形態において、式（２）最適化問題を解くために、ＥＭＤを用いたが、これに限定されない。以下に説明するＥＭＤ−Ｌ1又はＥＭＤ−Ｌ1＋が用いられてもよい。
【００７５】
Ａ．ＥＭＤ−Ｌ1
ＥＭＤ−Ｌ1は、例えば「H. Ling and K. Okada. An efficient earth mover's distance algorithm for robust histogram comparison. Pattern Analysis and Machine Intelligence, IEEE Transactions on, 29(5):840-853, 2007.」により提案されているより効果的なＥＭＤである。ＥＭＤ−Ｌ1における未知変数の数は、ヒストグラムにおいてＮ個のビンがあるときに、通常のＥＭＤにおけるＯ（Ｎ2）からＯ（Ｎ）に減る。従って、計算速度の向上が期待できる。
【００７６】
Ｂ．ＥＭＤ−Ｌ1＋
ＥＭＤ−Ｌ1＋は、上記の文献において提案されているＴｒｅｅ−ＥＭＤを本願の発明者らが拡張したものである。当該拡張は、ＥＭＤ−Ｌ1における不必要なセル（ユニットエリア）を無視し、出力を単純化するという２点である。
【００７７】
１点目の拡張に関し、不必要なセルとは、海、山、森林などのエージェントが移動しない（もしくは、移動していたとしても少数である）エリアにおけるセルである。ＥＭＤ−Ｌ1＋では、このようなセルを処理の対象エリアから除く。その結果、処理速度が向上する。
【００７８】
２点目の拡張では、生成されるツリーの数を２つに限定することである。Ｔｒｅｅ−ＥＭＤの変数は、隣り合う２つのビンの間のフローの量を示すので、フローの表現は冗長である。このように冗長性を削減することにより、処理速度を向上させることができる。
【符号の説明】
【００７９】
１０  推定装置、１１  制御部、１２  ＣＰＵ、１３  メモリ、１４  通信部、１５  記憶部、１６  入出力部
【図１】
0006632054.tif 000029
【図２】
0006632054.tif 000030
【図３】
0006632054.tif 000031
【図４】
0006632054.tif 000032