WO2019069865A1

WO2019069865A1 - パラメータ推定システム、パラメータ推定方法およびパラメータ推定プログラム記録媒体

Info

Publication number: WO2019069865A1
Application number: PCT/JP2018/036731
Authority: WO
Inventors: 隆文梶原; 山崎　啓介
Original assignee: 日本電気株式会社; 国立研究開発法人産業技術総合研究所
Priority date: 2017-10-03
Filing date: 2018-10-01
Publication date: 2019-04-11
Also published as: US20200250586A1; JP6853968B2; US11636399B2; JPWO2019069865A1

Abstract

第１のサンプル取得部は事前分布からパラメータのサンプルを取得する。関数実行部はサンプルを関数に供給してある分布からデータを取得する。類似度算出部はカーネル関数を用いてデータと正解データとの間の類似度を算出する。カーネル平均構築部は類似度とサンプルとカーネル関数とから事後分布のカーネル平均を構築する。第２のサンプル取得部はカーネル平均とカーネル関数とからパラメータの新たなサンプルを取得する。サンプル評価部は、新たなサンプルから選択した１つのサンプルを関数に供給することによって得られた新たなデータと正解データとの間の差分が、所定の閾値より小さいか否かを判断する。サンプル評価部は、差分が所定の閾値より小さいと判断した場合に、選択したサンプルを、パラメータとして推定する。本発明により、関数の高次元のパラメータ推定が可能となり、計算時間を減少することが可能となる。

Description

パラメータ推定システム、パラメータ推定方法およびパラメータ推定プログラム記録媒体

　本発明は、関数のパラメータを推定するパラメータ推定システム、パラメータ推定方法およびパラメータ推定プログラム記録媒体に関する。

　関数のパラメータを推定したいという要望がある。例えば、パラメータが与えられたある地域における人口分布（人流パターン）であって、関数が人流シミュレータであったとする。その場合、その人流パターンにもっとも近いデータを生み出す、人流シミュレータのパラメータを求めたいという要望がある。ここで、関数が人流シミュレータである場合、「パラメータ」とは、その人流パターンの「集団の人数、どこから来たか（出発地点）、何時に来たか（出発時刻）、どこから帰っていくのか（最終到達地点）、どこを経由するのか（経由地）、経由値の滞在時間」から成る。尚、出発地点と最終到達地点とは同じであってもよい。

　そのような関数のパラメータを求める手法として、種々の手法が知られている。

　例えば、近似ベイズ計算法（ＡＢＣ：Approximate Bayesian Computation）が知られている。ＡＢＣとは、モンテカルロ法の一種で、ランダムに選んだパタメータの値から、パフォーマンスの高いパラメータを順次求めていく方法である。例えば、最もシンプルなタイプのＡＢＣでは、関数から疑似データを生成し、その擬似データが現実のデータと近ければ、関数に使用されたパラメータを正しいものとして認める。

　ＡＢＣでは、事前知識に基づいてサンプリングを決めてシミュレーションを行う。そして、ＡＢＣでは、結果（擬似データ）が現実のデータにうまく適合するものだけをフィルタリングしている。

　しかしながら、ＡＢＣでは、パラメータが高次元であるとなかなか適合しないという問題がある。そのため、ＡＢＣでは、多数のサンプルが必要となる。

　Kernel ABCは、全サンプルを活用して、適合度で重み付けしたカーネル平均（kernel mean）として事後分布を表現する（非特許文献１参照）。多くの場合、事後分布からサンプリングを行いたいという動機が存在する。しかしながら、事後分布のkernel meanからのサンプリングは容易ではない。

　Kernel Herdingは、事後分布のkernel meanからのサンプリングを可能にする。特に、Kernel Herdingによるサンプリングは分布の特徴を少ない数でうまくカバーでき、スーパサンプルと呼ばれる。

　Automatic Relevance Determination (ARD)は、与えられたパラメータ集合から重要なものを残し、冗長性を取り除く有効なツールである。

　一方、パラメータとして人流パターンを推定する方法が種々提案されている。

　例えば、特許文献１は、都市の効率的な運用を支援する人流分析システムを開示している。特許文献１に開示された人流分析システムは、地域指定部と、人流算出部と、表示制御部とを備える。地域指定部は、地図上の地域を指定する。人流算出部は、時刻および位置の情報を含んだ複数の人流情報に基づき、地域指定部の指定地域に関わる人流を算出する。表示制御部は、表示部に表示された地図上に、人流算出部が算出した主たる人流の動線を表示させる。人流分析システムは、シミュレーション部を備える。シミュレーション部は、人流をシミュレーションによって算出するものであり、例えば、人流データベースに基づき、新たにエージェントを生成した場合の人流を算出する。ここでエージェントとは、環境の変化に基づく人流や交通流の変化を発生させるため、実際の環境の代わりに設定される仮想的な存在であり、例えば居住者の変化や施設の開設や閉鎖などをいう。

　また特許文献２は、時間的及び空間的相関を持つ時空間変数の値を精度良く予測することができる、時空間変数予測装置を開示している。特許文献２では、時空間変数の時系列データとして、人口分布、人流・交通流の速度・向き、金やダイヤモンドなど鉱物資源の埋蔵量、降水量などの気象データ、土地価格などを想定している。特許文献２では、実施の形態として、人口密度分布の時系列データが観測データとして与えられた条件の下で、未観測地点あるいは未来の時空間変数分布を推定・予測している。

　特許文献２に開示された時空間変数予測装置は学習部を備える。学習部は、人口密度情報記憶部に格納された観測データの集合に基づいて、観察データ同士の類似度を定義する関数であるカーネル関数の各々のハイパーパラメータと、観測データの各々に対する、複数のガウス過程の各々の寄与度を表すパラメータである負担率とを学習する。学習部は、負担率推定部と、ガウス過程パラメータ推定部と、反復判定部とを含む。負担率推定部は、観測データの集合と、複数のガウス過程のカーネル関数の各々のハイパーパラメータとに基づいて、複数のガウス過程からなる複数のユニットの各々の寄与度を表すパラメータであるユニット負担率と、複数のガウス過程の各々の寄与度を表すパラメータである負担率を推定する。ガウスパラメータ推定部は、観測データの集合と、複数のユニットの各々のユニット負担率と、複数のガウス過程の各々の負担率とに基づいて、複数のガウス過程の各々に対し、ガウス過程のカーネル関数の各々のハイパーパラメータを推定する。反復判定部は、予め定められた反復終了条件を満たすまで、負担率推定部による推定、およびガウス過程パラメータ推定部による推定を繰り返す。

　さらに特許文献３は、大規模な人流分布推定処理をより効率的に実行する人流調査支援システムを開示している。この特許文献３に開示された人流調査支援システムでは、分散処理プラットフォームを利用して人流の変化を表現する人流モデルにもとづく人流分布の仮説を大量生成している。そして、これらの各仮説と計測データとの適合性を評価することにより、計測データと矛盾しない人流分布を算出している。この際、人流モデルにもとづく人流分布の時間発展処理を類似する人流モデルごとに計算することで、計算を並列化しかつ重複計算を排除している。この人流モデル毎の処理並列化によって計測データとの時空間索引によって実測値との評価を高速化している。

　特許文献３に開示された人流調査支援システムは、複数のコンピュータに計算を分散する分散処理プラットフォームと、人流分布の仮説データから人流分布を推定した人流分布データを配信する機能を持つ人流分布配信システムと、人流のリアルタイムの計測データを計測する人流計測器とを備える。特許文献３において、分散処理プラットフォームは、人流全数調査データのうち現在時刻と一致するもの全体を所定数分複製し、各仮説の各人物の位置をランダムに少しずらし、人流仮説データとしている。

　また、特許文献４は、実店舗における商品の購入を支援する商品情報提供システムを開示している。特許文献４では、カメラもしくはＧＰＳ（Global Positioning System）情報から動線を分析・推定することで、来訪者の興味を把握し、推薦商品を決定している。

　非特許文献２は、モバイル空間統計を活用して、２５０ｍメッシュ単位でエリア間での回遊行動を分析する、技術的思想を開示している。

　非特許文献３は、アンケート情報を基にして、階層的な回遊モデルで立ち寄り場所を確率的に推定する、技術的思想を開示している。

国際公開第２０１６／０６７３６９号特開２０１７－０３３１９８号公報特開２０１４－１１２２９２号公報特開２０１５－０２２３５４号公報

Shigeki Nakagome, Shuhei Mano and Kenji Fukumizu (2013) Kernel Approximate Bayesian Computation in Population Genetic Inferences, Statistical Applications in Genetics and Molecular Biology, Volume 12, Issue 6, Pages 667-678 都心の歩行者回遊分析（その２）～都心滞在者の属性について～、名古屋市都市センタ　研究報告書、No. 117 2015. まちなか回遊行動の詳細分析と政策シミュレーションのための予測モデル、土木学会論文集、Vol. 71, No.5 323-335, 2015

　先行技術（ＡＢＣ）には、次に述べるような課題がある。

　第１の課題は、推定すべき関数のパラメータ数が多いと、パラメータの精度が落ちることである。

　また、ＡＢＣでは、パラメータ候補について、非効率に探索している。したがって、第２の課題は、あまり意味のないパラメータ候補についてもシミュレーションを行うので、無駄な計算時間がかかってしまうことである。

　さらに、第３の課題は、実際には必要のないパラメータを間違った知識のもとに探索したとき、それを無条件に受け入れてしまうことである。

　前述した先行技術（特許文献１～４、非特許文献１～３）には、それぞれ、次に述べるような問題がある。

　特許文献１は、単に、人流データベースに基づき、新たにエージェントを生成した場合の人流を算出する技術的思想を開示しているに過ぎない。従って、特許文献１では、各個人の外出行動を、予め人流データベースに格納しておく必要がある。換言すれば、特許文献１は、そのような予め各個人の外出行動が不明である場合には成り立たない技術的思想である。

　特許文献２は、単に、人口密度分布の時系列データが観測データとして与えられた条件の下で、未観測地点あるいは未来の時空間変数分布を推定・予測する技術的思想を開示しているに過ぎない。また、特許文献２は、時系列データとして人流の速度・向きを対象としたものに適用することができることを記載しているが、具体的にどのように人流の速度・向きに適用するかについては、何ら開示も示唆もしていない。

　特許文献３では、分散処理プラットフォームは、人流全数調査データのうち現在時刻と一致するもの全体を所定数分複製し、各仮説の各人物の位置をランダムに少しずらし、人流仮説データとしている。そのため、各人物ごとに膨大なデータが必要となる。

　特許文献４では、カメラやＧＰＳの情報を活用しているので、個人の行動軌跡を取得することができる。しかしながら、特許文献４の技術では、個人特定につながるため、プラバシーの侵害に当たる可能性が高い。

　非特許文献１に開示されたKernel ABCでは、前述したように、kernel meanを使って再度サンプリングできないという問題がある。すなわち、漸近手法を完成できない。

　非特許文献２では、モバイル空間統計等の人口分布情報を用いているので、プライバシーを保護することが可能である。しかしながら、非特許文献２の開示する方法は、分析粒度が大まかであるため、細かな回遊行動を推定することが困難である。

　非特許文献３に開示された方法では、確率的な行動モデルに基づき、人口分布とアンケート情報を基に回遊行動を分析している。しかしながら、そのような分析方法では、詳細な施設間の関連性を把握することができないため、その分析結果を街作りに活用するには不十分な情報となってしまう。

　本発明の目的は、上述した課題を解決する、パラメータ推定システムおよびパラメータ推定方法を提供することにある。

　本発明の１つの態様として、パラメータ推定システムは、データｙを出力する関数のパラメータθを推定するパターン推定システムであって、前記関数に設定する前記パラメータθを、繰り返し学習により決定する機械学習装置と、を備え、前記機械学習装置は、事前分布Ｐから前記パラメータθのＮ（Ｎは２以上の整数）個のサンプルθ_ｉ（１≦ｉ≦Ｎ）を取得する第１のサンプル取得部と；前記Ｎ個のサンプルθ_ｉを順次前記関数に供給することで、それぞれ、前記関数からＮ個のデータｙ_ｉを取得する関数実行部と；カーネル関数ｋを用いて、前記Ｎ個のデータｙ_ｉと正解データｙ^＊との間のＮ個の類似度ｗ_ｉを算出する類似度算出部と；該Ｎ個の類似度ｗ_ｉと、前記Ｎ個のサンプルθ_ｉと、前記カーネル関数ｋとから、事後分布のカーネル平均μを構築するカーネル平均構築部と；前記カーネル平均μと前記カーネル関数ｋとから、パラメータθのＴ（Ｔは２以上の整数）個の新たなサンプルθ_ｔ（１≦ｔ≦Ｔ）を取得する第２のサンプル取得部と；該Ｔ個の新たなサンプルθ_ｔから選択した１つのサンプルθ_ｋを前記関数に供給することによって前記関数から得られた新たなデータｙ_ｋと前記正解データｙ^＊との間の差分｜｜ｙ_ｋ－ｙ^＊｜｜が、所定の閾値εより小さいか否かを判断するサンプル評価部と；を備え、前記サンプル評価部は、前記差分｜｜ｙ_ｋ－ｙ^＊｜｜が前記所定の閾値εより小さいと判断した場合に、前記選択したサンプルθ_ｋを、前記パラメータθとして推定し、出力する。

　本発明の他の態様として、パターン推定方法は、関数のパラメータθを推定するパラメータ推定システムにおける、前記パラメータθのサンプルから、データｙを出力する前記関数を用いて、前記関数に設定する前記パラメータθを、繰り返し学習により決定するパラメータ推定方法であって、第１のサンプル取得部が、事前分布Ｐから前記パラメータθのＮ（Ｎは２以上の整数）個のサンプルθ_ｉ（１≦ｉ≦Ｎ）を取得し；関数実行部が、前記Ｎ個のサンプルθ_ｉを順次前記関数に供給することで、それぞれ、前記関数からＮ個のデータｙ_ｉを取得し；類似度算出部が、カーネル関数ｋを用いて、前記Ｎ個のデータｙ_ｉとデータを表す正解データｙ^＊との間のＮ個の類似度ｗ_ｉを算出し；カーネル平均構築部が、前記Ｎ個の類似度ｗ_ｉと、前記Ｎ個のサンプルθ_ｉと、前記カーネル関数ｋとから、事後分布のカーネル平均μを構築し；第２のサンプル取得部が、前記カーネル平均μと前記カーネル関数ｋとから、パラメータθのＴ（Ｔは２以上の整数）個の新たなサンプルθ_ｔ（１≦ｔ≦Ｔ）を取得し；サンプル評価部が、前記Ｔ個の新たなサンプルθ_ｔから選択した１つのサンプルθ_ｋを前記関数に供給することによって前記関数から得られた新たなデータｙ_ｋと前記正解データｙ^＊との間の差分｜｜ｙ_ｋ－ｙ^＊｜｜が、所定の閾値εより小さいか否かを判断し；前記サンプル評価部は、前記差分｜｜ｙ_ｋ－ｙ^＊｜｜が前記所定の閾値εより小さいと判断した場合に、前記選択したサンプルθ_ｋを、前記パラメータθとして推定し、出力する。

　本発明によれば、関数の高次元のパラメータ推定が可能になり、計算時間を減少することが可能となる。

本発明の一実施形態に係る人流パターン推定システムの構成を示すブロック図である。正解データを示す人口時間分布データの一例を示す図である。図１に示した人流パターン推定システムに使用される、人流シミュレータに組み込まれる（設定される）、地図情報および施設情報（地図／建物データ）の一例を示す図である。人流パターンの一例を示す図である。図１に示した人流パターン推定システムに使用される、機械学習装置で実現される主な処理部を示すブロック図である。図５に示した機械学習装置に使用される、データ入力部の動作を説明するための図である。図５に示した機械学習装置に使用される、第１のサンプル取得部の動作を説明するための図である。図５に示した機械学習装置に使用される、シミュレーション実行部の動作を説明するための図である。図５に示した機械学習装置に使用される、Kernel ABC 実行部（類似度算出部）の動作を説明するための図である。図５に示した機械学習装置に使用される、カーネル平均構築部の動作を説明するための図である。図５に示した機械学習装置に使用される、第２のサンプル取得部の動作を説明するための図である。図５に示した機械学習装置に使用される、サンプル評価部の動作を説明するための図である。図１に示した人流パターン推定システムの動作の一実施例について説明するためのフローチャートである。図１に示した人流パターン推定システムを実行した場合の実験に使用した、４ｍｅｓｈモデルの一例を示す図である。図１に示した人流パターン推定システムを実行した場合の実験に設定した、５つのルールの一例を示す図である。図１４に示した４ｍｅｓｈモデルの実験結果を示す図である。

　まず、本発明の理解を容易にするために、本明細書中で使用する用語の意味に関して簡略的に説明する。

　近似ベイズ計算法（ＡＢＣ：Approximate Bayesian Computation）とは、モンテカルロ法の一種で、ランダムに選んだパラメータの値から、パフォーマンスの高いパラメータを順次求めていく方法である。例えば、ＡＢＣでは、関数から疑似データを生成し、その擬似データが現実のデータと近ければ、関数に使用されたパラメータを正しいものとして認める。

　換言すれば、ＡＢＣは、観測データと関数から得られたデータの差異について拒絶アルゴリズム法を適用する、ベイズ推定のための尤度無しの方法である。

　Kernel ＡＢＣは、全サンプルを活用して、適合度で重み付けしたカーネル平均（kernel mean）として分布を表現する。

　カーネル法は、高次元空間への写像を構成する方法論である。カーネル法は、データを高次元の再生核ヒルベルト空間（一般には無限次元）へ写像し、解析しやすいデータに変換する。ヒルベルト空間（Hilbert space）とは、ベクトル空間であり、内積が備わっており、完備であるような集合である。

　集合Ω上の関数を要素に持つヒルベルト空間Ｈが再生核ヒルベルト空間（ＲＫＨＲ：reproducing kernel Hilbert space）であるとは、任意のｘ∈Ωに対してφ_Ｘ∈Ｈがあって、任意のｆ∈Ｈに対し、
　　　　　　＜ｆ，φ_Ｘ＞　＝　ｆ（ｘ）
が成り立つことをいう。φ_Ｘのことを再生核という。

　Kernel Herdingは、kernel meanの分布の特徴をうまく利用し、カブリの少ない代表点をスーパサンプルとして選択する。つまり、Kernel Herdingは、少ないサンプル点で、適切な設定の分布を表現する。

　ARDでは、ハイパーパラメータを調整することにより、事前分布からのサンプルでパラメータのほとんどの値が０を取ることで、余分なものを削除している。通常、これはディリクレ分布やベータ分布を事前分布として導入することにより達成される。ディリクレ分布の例で言えば、それに付随するハイパーパラメータの値を小さくとると、分布からのサンプルがほとんど0に近い値を取り、少ないパラメータで現象を説明しようとする力が働く。

　別言すると、パラメータの冗長性をコントロールしたい場合には、パラメータのサンプルをディリクレ分布もしくはベータ分布から取得する。ディリクレ分布、ベータ分布にはパラメータの冗長性をコントロールするハイパーパラメータが存在する。そのハイパーパラメータの値を変更することにより、パラメータのサンプルがもつ冗長性を調整する。例えば、冗長性をなくすようにハイパーパラメータを設定すると、サンプルの要素のごく少数のみが非０の値を取り、多くの要素は０の値を取る。これにより、０の値を取った要素は実質影響を持たない。こうして、得られた事後分布のカーネル平均は冗長性が非常に減少したものとなる。

［関連技術］
　本発明の理解を容易にするために、関連技術とその問題点について説明する。

　問題の事実（真実）が与えられたときに、その事実（真実）に十分に近接したデータを生成する、関数のパラメータの集合を知りたいとする。ベイズ規則を適用したいが、関数は評価されえない。しかしながら、関数の尤度ｐ（ｙ｜θ）から事前確率π（θ）に基づいて、次の式により、サンプリングは可能である。

　そこで、ＡＢＣのフレームワークを使用することが考えられる。表１に、ＡＢＣのアルゴリズムを示す。

　しかしながら、ＡＢＣには、次のような課題がある。

　第１の課題は、データの次元が高いと、パラメータが採択基準(tolerance level)をなかなか満たさず、現実的な時間でアルゴリズムが終了しないことである。

　第２の課題は、パラメータの次元が高いと、ＡＢＣの性能はすぐに低下することである。

　第３の課題は、関数を評価するのに時間がかかるので、実行数を省けないことである。

［本発明に係るアルゴリズムの概要説明］
　本発明では、第１および第２の課題を解決するためにKernel ABCを使用し、第３の課題を解決するためにKernel Herdingを使用している。

　本発明のアルゴリズムの概要は次の通りである。
Ａ１．事前分布からパラメータの組をＮ回取り出す。
Ａ２．Ｎ個のパラメータの組を使用して関数からデータをＮ回生成する。
Ａ３．Kernel ABC
・事後分布のカーネル平均を計算することによって、観測データに近接するデータを生成す関数のパラメータを推察する。
Ａ４．Kernel Herding
・次の繰り返しのために事後分布のカーネル平均から有益なサンプルをＴ回取り出す。
Ａ５．収束するためにＡ２．から繰り返す。

　次に、Kernel ABCについて説明する。
Kernel ABCは、Kernelベイズ規則に基づいており、ベイズ規則のノンパラメトリックな実現である。
Ｂ１．事前分布からパラメータの組をＮ回取り出す。
Ｂ２．パラメータの組を関数にＮ回入力することによって、Ｎ個の関数を評価したデータを取得する。
Ｂ３．真のデータと各関数を評価したデータとの間の類似度（重み）を計算する。
Ｂ４．Ｂ１．でのパラメータの組をそれらの類似度（重み）で平均化することによって、下記数式に従って、パラメータ事後分布のカーネル平均を計算する。

　次に、カーネル平均について説明する。
　Ｘを測度空間Ω上の値を取るランダム変数とし、ｋを測度空間Ω上の正定値カーネルであり、Ｈが正定値カーネルｋによって規定された再生核ヒルベルト空間（ＲＫＨＳ：reproducing kernel Hilbert space）であるとする。このとき、カーネル平均は、次式で求められる。

　特性カーネルによって、確率分布を再生核ヒルベルト空間に一意に埋め込むことができる。特性カーネルは、例えば、ガウスカーネルやラプラスカーネルなどである。

　次に、カーネル平均とカーネルベイズ規則を使用することの利点について説明する。変数間のパラメトリックモデルや、密度関数を想定することなく、ベイズ規則をノンパラメトリックに適用できる。したがって、確率分布を表すカーネル平均は、データからノンパラメトリックに推定される。

　次に、Kernel Herdingについて説明する。
　Kernel ABCの出力は、再生核ヒルベルト空間内の事後分布のカーネル平均であるが、それ自身事後分布ではない。したがって、マルコフ連鎖モンテカルロ法に代表される事後分布からのサンプリング手法を使うことができない。しかし、Kernel Herdingは、事後分布からのサンプルを可能とする。

　更に、Kernel Herdingからのサンプルは、分布の特徴を上手く反映したスーパサンプルと呼ばれる。

　Kernel Herdingは、次の誤差を貪欲に最小化することで、サンプリングを行っていると解釈できる。

　Kernel Herdingは、分布において密度が高い領域からサンプルを取る。つまり、分布のモードを探索すると考えられる。また、一度サンプルが取られた領域からは離れて、次のサンプルが取得される。
　以上を要約すると、次のようになる。

　Kernel ABCは、比較的高い次元のパラメータを処理でき、比較的高い次元のデータ（関数を評価したデータ／問題の真のデータ）を処理でき、扱いやすい計算的な複雑性を持つ。

　Kernel Herdingは、必要とされる関数を評価する回数を削減することができる。

［発明の概要の説明］
　以下では、関数が人流シミュレータである場合を例に挙げて説明するが、本発明は関数が人流シミュレータである場合に限定されないのは勿論である。人流シミュレータへの適用についてそもそも、モバイル空間統計のような、人口分布の時系列に関しての人流シミュレータのデータ同化はこれまで試みられていない。そのときに、以下の課題があり、本発明に係る提案手法ではそれに対処している。

　１）人流シミュレータはマルチエージェントシミュレータであり、多くの場合シミュレーションに非常に計算時間がかかる。

　これに対して、本提案手法では、Kernel herding を用いることにより、より少ない回数のシミュレーションでデータ同化を達成できる。これは、kernel herding が効率的により多くの情報を持つ重要なパラメータのサンプルを取得するからである。

　また、これについて、Ｏ(１／Ｔ) で誤差が収束していく。一方、モンテカルロ法などによるパラメータのサンプル取得では、より遅いＯ(１／√Ｔ) である。

　２）人流シミュレータではないが、他のデータ同化に関して、ＡＢＣなどを使ったベイズ推定では、パラメータの次元が高次元の際に、上手く働かないことが分かっている。

　これに対して、本提案手法では、Kernel ABC が、再生核ヒルベルト空間でベイズ推定を行うことにより、高次元パラメータに対応している。

　３）事前知識が間違っていた場合に、ある地域において存在しない人流パターンを想定して、事前分布を構築してしまうことがある。例えば、銀座において、吉野家で昼食をとり、ロレックスで時計を買い、フェラーリで車を買い、松屋で夜ご飯を食べるなどの人流パターンである。関連技術では、このような存在すると考えにくい人流パターンをアルゴリズムで自動的に削除する仕組みはなかった。

　これに対して、本提案手法では、ＡＲＤにより、そういったパターンを削除するスパース性が導入されている。

　アルゴリズムについて、本提案手法は、関数が評価できない場合に、点推定を行うことを可能としている。

　関連技術では、関数が評価できない場合にベイズ推定を行うことはできた（ＡＢＣなど）。つまり、それによって正しい事後分布を求めることができた。しかしながら、ＡＢＣでは、真値が一点ある状況での点推定を、関数が評価できない場合に行うことはできなかった。

　これに対して、本提案手法では、Kernel ABC + Kernel Herding で取得された事後分布からのサンプルを、事前分布からのサンプルとみて、Kernel ABC+ Kernel herding を繰り返し適用することで、点推定を達成している。

　一度目の適用では以下の数式で表される事後分布がえられる。

　このＰ(θ｜ｙ) を事前分布として、二度目の適用をすると、以下の数式が得られる。

　これは、

ということであり、Ｎ回目の適用では、

となる。つまり、関数ｌ(ｙ｜θ) がＮ乗され、影響が強くなり、一方で事前分布π(θ) の影響が少なくなっている。これを無限回繰り返したあと、Ｐ(θ｜ｙ)を最大化するθを取得することは、π(θ)の影響が消えているため、ｌ(ｙ｜θ)を最大化するθを取得することと一致すると考えられる。つまり、無限回のパラメータに関するベイズ推定を同一のデータに対して繰り返し、Ｐ(θ｜ｙ)を最大化するθを取得することは点推定の最尤法と一致すると考えることができる。これを踏まえ、関数が評価できないときにベイズ推定を行う手法であるABCを使うことで、関数なしの最尤法を実行できる。これをRecursive ABC-MLEと名づけた。このアルゴリズムの内部でkernel ABCを本明細書では用いているが、関数無しに事後分布を求めることができる手法であればどれでも構わない。また、このアルゴリズム内部でＰ(θ｜ｙ)を最大化するθを取得することは、Kernel Herdingのmode-seeking（最頻値追跡）という特性により可能となっている。しかしながら、取得する手法は、必ずしもKernel Herdingである必要はなく、最頻値をサンプリングする特性のある手法であればよい。

　Kernel Herding が、事後分布から効率的により多くの情報を持つ重要なパラメータのサンプルを取得する。これについて、Ｏ(１／Ｔ ) で誤差が収束していく。

　一方、上記のベイズ推定の繰り返しを、関連技術の手法のマルコフ連鎖モンテカルロ法などによるパラメータのサンプル取得で行った場合は、より遅いＯ(１／√Ｔ) である。

　またこのアルゴリズムの本質は、関数が評価できない場合に最尤法を可能にすることにある。このアルゴリズムの適用先としては、後で詳述する人流パターンを推定することのみでなく、次の関数のパラメータを推定するものも考えられる。

１）株価など金融時系列の価格変動の分散（ボラティリティ）を適切にモデル化するためによく使われるアルファ定常分布のパラメータをデータに合わせて推定すること。

２）ハエなど、生物の個体数の時間発展を記述するモデルのパラメータをデータに合わせて推定すること。

３）集団遺伝学において、与えられたデータの遺伝的多様性を実現するようなモデルのパラメータ（突然変異率など）を推定すること。

　尚、このアルゴリズムは、ソフトウェア関連発明であるとしても、計算時間を減少できる効果を奏するので、コンピュータ関連技術を改善（向上）するものであり、抽象概念には該当しない。

［実施の形態］
　以下に述べる実施の形態では、本発明を人流パターン推定システムに適用した場合について説明するが、上述したように、本発明は、一般的に関数のパラメータを推定するパラメータ推定システムにも適用可能である。図１は、本発明の一実施形態に係る人流パターン推定システム１００の構成を示すブロック図である。

　図示の人流パターン推定システム１００は、データを処理するデータ処理装置２００と、後述するプログラムやデータを記憶する記憶装置３００と、データを入力する入力装置４００と、データを出力する出力装置５００と、を備えている。

　出力装置５００は、ＬＣＤ（Liquid Crystal Display）やＰＤＰ（Plasma Display Panel）などの表示装置やプリンタからなる。出力装置５００は、データ処理装置２００からの指示に応じて、操作メニューなどの各種情報を表示したり、最終結果を印字出力する機能を有する。

　記憶装置３００は、ハードディスクやリードオンリメモリ（ＲＯＭ）およびランダムアクセスメモリ（ＲＡＭ）などのメモリからなる。記憶装置３００は、データ処理装置２００における各種処理に必要な処理情報３１０やプログラム３２０を記憶する機能を有する。

　データ処理装置２００は、ＭＰＵ（micro processing unit）などのマイクロプロセッサや中央処理装置（ＣＰＵ）からなる。データ処理装置２００は、記憶装置３００からプログラム３２０を読み込んで、プログラム３２０に従ってデータを処理する各種処理部を実現する機能を有する。

　データ処理装置２００は、人流シミュレータ２１０と機械学習装置２２０とを備える。本実施の形態では、関数が人流シミュレータ２１０である場合を例に挙げて説明するが、前述したように、本発明は、一般的に、他の関数にも適用可能である。

　図示の人流パターン推定システム１００において、アルゴリズムに必要なデータや、装置モバイル空間統計のような、人流を推定したい地域、時刻においての人口分布のデータが与えられているものとする。これを正解データｙ^＊と呼ぶ。

　図２は、正解データｙ^＊を示す人口時間分布データの一例を示す図である。人口時間分布データは、時刻Ｔ＝・・・、ｔ－１，ｔ、ｔ＋１、ｔ＋２、・・・ごとの、人口分布を示している。時刻Ｔは、例えば、３０分毎の時刻である。このような人口時間分布データは、例えば、統計情報やＳＮＳ（Social Networking Service）で取得することが可能である。ここで、ＳＮＳとは、人と人との社会的な繋がりを維持・促進する様々な機能を提供する、会員制のオンラインサービスである。記憶装置３００は、処理情報３１０として正解データｙ^＊を保存している。

　人流パターン推定システム１００は、街中における人流グループ単位での動きである人流パターンを推定するシステムである。上述したように、本発明は、一般的に、関数のパラメータを推定するパラメータ推定システムに適用可能である。

　人流シミュレータ２１０は、人流パターンを、任意の地域、時刻においてシミュレーションすることができる。

　詳述すると、人流シミュレータ２１０には、各地点、各時刻での人口分布の情報および前記街中の地図情報や施設情報が組み込こまれる。

　図３は、人流シミュレータ２１０に組み込まれる（設定される）、地図情報および施設情報（地図／建物データ）の一例を示す図である。記憶装置３００は、処理情報３１０として地図情報および施設情報を保存している。これら地図情報および施設情報は、後述するように、人流シミュレータ２１０に設定される。

　人流シミュレータ２１０は、後述するように機械学習装置２２０から供給される、人流パターンを規定するパラメータθのサンプルから、各地点、各時刻での人口分布データｙを出力する。ここで、人流シミュレータ２１０から得られた人口分布データｙは、通常、正解データｙ^＊には一致しないことに注意されたい。何故なら、正解データｙ^＊を満たす人流パターンを誰も知らないからである。一般的に、本発明に係る関数は、パラメータθのサンプルから、データｙを出力する。

　図４に、人流パターンの一例を示す。図４に示す人流パターンは、地点Ａ、Ｂ、Ｃ、Ｄ、Ｅ、Ｆ、Ｇにおける、時刻９：００以降の人流パターンを示している。

　パラメータθは、各集団の、出発点、出発時刻、経由点、滞在時間、目的地、および人数の情報を含む。

　機械学習装置２２０は、人流シミュレータ２１０に設定するパラメータθを、後述するように、繰り返し学習により決定する。

　図５は、機械学習装置２００で実現される主な処理部を示すブロック図である。

　機械学習装置２００は、データ入力部２２１と、第１のサンプル取得部２２２と、シミュレーション実行部２２３と、Kernel ABC実行部２２４と、カーネル平均構築部２２５と、第２のサンプル取得部２２６と、サンプル評価部２２７とから成る。なお、Kernel ABC実行部２２４は、後述するように類似度算出部として働く。

　次に、図６乃至図１２を参照して、各処理部の動作について説明する。

　図６を参照して、データ入力部２２１の動作について説明する。

　図６に示されるように、データ入力部２２１は、人流シミュレータ２１０に人流パターンを推定したい地域の情報を入力する。詳述すると、記憶装置３００は、処理情報３１０に、前述したように、地域の情報として地図情報や施設情報を保存している。データ入力部２２１は、記憶装置３００から地図情報と施設情報とを読み出して、それら情報を人流シミュレータ２１０に入力（設定）する。

　尚、本例では、データ入力部２２１は、記憶装置３００に保存されている地図情報および施設情報を読み出して、人流シミュレータ２１０に設定しているが、入力装置４００から供給された地図情報および施設情報を、直接、人流シミュレータ２１０に設定してもよい。

　図７を参照して、第１のサンプル取得部２２２の動作について説明する。第１のサンプル取得部２２２は、「事前分布からのサンプル取得部」とも呼ばれる。この第１のサンプル取得部２２には、ＡＲＤが関わっている。

　図７に示されるように、記憶装置３００には、処理情報３１０として、事前分布Ｐと、その事前分布Ｐのハイパーパラメータηと、事前分布Ｐからのサンプルの個数Ｎとが格納されている。事前分布Ｐは、任意のパラメータサンプルの合計が１となるような、mixture モデル、もしくは解の構造がmixture モデルと同じ場合を反映できる適当な事前分布である。

　この事前分布Ｐのハイパーパラメータηを調整することにより、このアルゴリズムの最終的な解にスパース性を入れるかどうかが決定される。

　第１のサンプル取得部２２２は、下記の数式のように、事前分布からサンプルを取り出す。

　したがって、第１のサンプル取得部２２２は、人流シミュレータ２１０へのパラメータのＮ個のサンプル{θ_ｉ}（ｉ＝１～Ｎ）を取得する。この取得したパラメータのＮ個のサンプル{θ_ｉ}は、処理情報３１０として記憶装置３００に保存される。

　図８を参照して、シミュレーション実行部２２３の動作について説明する。

　図８に示されるように、記憶装置３００は、前述したサンプルの個数Ｎと、パラメータのＮ個のサンプル{θ_ｉ}と、を処理情報３１０として保存している。

　シミュレーション実行部２２３は、人流シミュレータ３１０にパラメータのサンプル{θ_ｉ}をＮ回入力して、人流シミュレータ３１０からＮ個の人口分布データ{ｙ_ｉ}（１≦ｉ≦）を得る。この得られたＮ個の人口分布データ{ｙｉ}は、処理情報３１０として記憶装置３００に保存される。

　図９を参照して、Kernel ABC 実行部（類似度算出部）２２４の動作について説明する。

　図９に示されるように、記憶装置３００は、前述したＮ個の人口分布データ{ｙ_ｉ}と、前述した正解データｙ^＊と、カーネル関数ｋと、正則化項δと、を処理情報３１０として保存している。

　Kernel ABC 実行部２２４は、下記の数式で表される計算を行い、正解データと各シミュレーションされた人口分布データ{ｙ_ｉ}との類似度{ｗ_ｉ}（１≦ｉ≦Ｎ）を得る。

ここで、Ｇはグラム行列を示し、Ｉは単位行列を示す。この得られたＮ個の類似度{ｗ_ｉ}は、処理情報３１０として記憶装置３００に保存される。

　なお、本例では、類似度算出部として、Kernel ABCアルゴリズムにより前記類似度ｗ_ｉを算出するKernel ABC 実行部２２４を用いているが、本発明はこれに限られない。すなわち、類似度算出部で使用するアルゴリズムは、後述する事後分布のカーネル平均が構築できるようにＮ個の類似度{ｗ_ｉ}を出力するアルゴリズムなら何でも構わない。

　図１０を参照して、カーネル平均構築部２２５の動作について説明する。

　図１０に示されるように、記憶装置３００は、前述したパラメータのＮ個のサンプル{θ_ｉ}と、前述したＮ個の類似度{ｗ_ｉ}と、前述したカーネル関数ｋと、を処理情報３１０として保存している。

　カーネル平均構築部２２５は、下記の数式で表される計算を行ない、事後分布のカーネル平均μを得る。

この得られたカーネル平均μは、処理情報３１０として記憶装置３００に保存される。

　図１１を参照して、第２のサンプル取得部２２６の動作について説明する。第２のサンプル取得部２２６は、「事後分布からのサンプル取得部」とも呼ばれる。この第２のサンプル取得部２２６は、上述したKernel Herdingを実行する。

　図１１に示されるように、記憶装置３００には、前述したカーネル平均μと、前述したカーネル関数ｋと、取り出すサンプル数Ｔと、を処理情報３１０として保存している。

　第２のサンプル取得部２２６は、以下の数式で現れる誤差を最小化する計算を行って、新たなパラメータのＴ個のサンプル{θ_ｔ}（１≦ｔ≦Ｔ）を得る。

この得られた新たなパラメータのＴ個のサンプル{θ_ｔ}は、処理情報３１０として記憶装置３００に保存される。

　最後に、図１２を参照して、サンプル評価部２２７の動作について説明する。

　図１２に示されるように、記憶装置３００には、前述した新たなパラメータのＴ個のサンプル{θ_ｔ}と、前述した正解データｙ^＊と、閾値εと、を処理情報３１０として保存している。

　サンプル評価部２２７は、先ず、新たなパラメータのＴ個のサンプル{θ_ｔ}から適当なサンプルθ_ｋを選び、その選択したサンプルθ_ｋを人流シミュレータ２１０へ供給する。それにより、サンプル評価部２２７は、人流シミュレータ２１０から新たな人口分布データｙ_ｋを取得する。引き続いて、サンプル評価部２２は、下記の式が満たされているか否かを判断する。

すなわち、サンプル評価部２２７は、新たな人口分布データｙ_ｋと正解データｙ^＊との差分｜｜ｙ_ｋ－ｙ^＊｜｜が、閾値εより小さいか否かを判断する。上記式が満たされていれば、サンプル評価部２２７は、選択したサンプルθ_ｋを、推定された人流パターンを規定する真のパラメータとして推定し出力し、処理を終了する。

　一方、上記式が満たされなかった場合、シミュレーション実行部２２３での動作に戻る。

　もしくは、第１のサンプル取得部２２２での動作に戻る。このとき、スパース性が十分に入っていなければ、第１のサンプル評価部２２２の入力であるハイパーパラメータηも変更されるものとする（ＡＲＤ）。

　次に、図１３のフローチャートを参照して、具体例を挙げながら、人流パターン推定システム１００の動作の一実施例について説明する。

　先ず、データ入力部２２１は、地図情報と施設情報とを人流シミュレータ２１０に設定する（ステップＳ１０１）。

　具体的には、データ入力部２２１は、人流シミュレータ２１０に銀座などのマップを与える。

　次に、第１のサンプル取得部２２２は、ハイパーパラメータηに基づいて事前分布Ｐからサンプルを取り出し、人流シミュレータ２１０へ入力すべきパラメータのＮ個のサンプル{θ_ｉ}を取得する（ステップＳ１０２）。

　具体的には、第１のサンプル取得部２２２は、銀座マップにおいてシミュレーションする人口集団を制御するパラメータθのサンプルを、例えば１０００個取ってくる（Ｎ＝１０００）。

　例えば、一つのパラメータθは、以下のような要素で構成される。銀座を訪れるサラリーマンの集団を考えると、パラメータθは、集団の人数、彼らが銀座を訪れる時刻、彼らが銀座を出ていく時刻、銀座で立ち寄る店舗の位置情報、その店舗での滞在時間などである。ハイパーパラメータηは、それを調整することによって余分なパラメータを削除することができる。例えば、銀座には、コンビニエンスストアしか立ち寄らない集団は来ないといった推定がこれによってできる。

　引き続いて、シミュレーション実行部２２３は、パラメータのサンプル{θ_ｉ}を、人流シミュレータ２１０にＮ回入力して、Ｎ個の人口分布データ{ｙ_ｉ}を得る（ステップＳ１０３）。

　具体的には、シミュレーション実行部２２３は、ステップＳ１０２で得られた一つ一つのパラメータのサンプルに関して、シミュレーションを行う。あるパラメータ集合のサンプルでは、集団の人数が１０００人かもしれず、別のサンプルでは２０００人かもしれないといった違いがある。ステップＳ１０２で、１０００個のパラメータ集合を取ってきたので、シミュレーション実行部２２３は１０００回のシミュレーションを行う。そして、人流シミュレータ２１０から、そのシミュレーションされた人口分布データｙ_ｉを取得する。

　次に、Kernel ABC 実行部（類似度算出部）２２４は、上記数１０の計算を行い、正解データｙ^＊と各シミュレーションされた人口分布データ{ｙ_ｉ}との類似度{ｗ_ｉ}を得る（ステップＳ１０４）。

　具体的には、Kernel ABC 実行部２２４は、銀座で１０００個の違った形でシミュレーションされたデータｙ_ｉを、銀座の正解データｙ^＊と比べる。この比較により、正解データｙ^＊とシミュレーションからのデータｙ_ｉとの類似度ｗ_ｉが計算される。Kernel ABC 実行部２２４は、この類似度ｗ_ｉをもとに、シミュレーションされたデータに紐づいているパラメータのサンプルθ_ｉのもっともらしさを決定する。

　引き続いて、カーネル平均構築部２２５は、上記数１１に従って、事後分布のカーネル平均μを計算する（ステップＳ１０５）。

　具体的には、カーネル平均構築部２２５は、ステップＳ１０４で計算された類似度ｗ_ｉを元に、パラメータのサンプルθ_ｉに関して、事後分布のカーネル平均μを構築する。これは、事前分布からのサンプルθ_ｉを元に、類似度ｗ_ｉを使い重みを再調整し、確率分布を再生核ヒルベルト空間上で構築することに相当する。

　次に、第２のサンプル取得部２２６は、上記数１２の誤差を最小化する計算を行なって、新たなパラメータのＴ個のサンプル{θ_ｔ}を取得する（ステップＳ１０６）。

　具体的には、第２のサンプル取得部２２６は、事後分布のカーネル平均μからのサンプリングを行う。このとき、サンプルされるのは新たなパラメータのサンプルθ_ｔだが、ステップＳ１０２において取得したパラメータのサンプルθ_ｉとは違い、正解データｙ^＊を生み出すようなパラメータのサンプルが求められる。

　引き続いて、サンプル評価部２２７は、新たなパラメータのＴ個のサンプル{θ_ｔ}を評価する（ステップＳ１０７）。

　具体的には、サンプル評価部２２７は、ステップ１０６で求められたパラメータ集合{θ_ｔ}（１≦ｔ≦Ｔ）のなかからパラメータθ_ｋを一つ適当に選び、それを人流シミュレータ２１０に入力し、データｙ_ｋを獲得する。そして、サンプル評価部２２７は、獲得したデータｙ_ｋを正解データｙ^＊と比べることで、パラメータの確からしさを評価する。選択したパラメータθ_ｋが充分に確からしい場合、サンプル評価部２２７は、それを銀座における人流推定値とする（ステップＳ１０８）。

　そのパラメータθ_ｋには、銀座にどういった集団が、どこから来て、どこに滞在し、どこから帰っていくかという情報が含まれている。パラメータθ_ｋの確からしさが満足のいくものでない場合（ステップＳ１０７のＮｏ）、人流パターン推定システム１００は、もう一度、ステップＳ１０３もしくはステップＳ１０２から動作を繰り返す。

　次に、人流パターン推定システム１００を実行した場合の実験結果について説明する。

　ここでは、図１４に示すように、銀座の４ｍｅｓｈモデルの実験を設定した。本例では、（２×２）の４つのメッシュに分割されているが、実際の例では、例えば、５×５＝２５メッシュに分割されることに注意されたい。図示の例では、４つのメッシュは、（０，０）、（０，１）、（１，０）、および（１，１）の領域から成る。時刻は、９００～１７：００で３０分間隔とした。

　図１５に示すような、５つのルールを設定した。そのため、推定パラメータの次元は、（Start, Goal, Transit×2, Stay×2, T, N）×５ルールで、４０次元となる。ここで、Tは集団の出発時刻を示し、Nは集団の人数を示す。なお、この各ルールにおいて、Nはパラメータθのサンプル数Ｎとは異なり、Tは新たなパラメータのサンプル数Ｔとは異なることに注意されたい。

　図１６は、図１４に示した銀座の４ｍｅｓｈモデルの実験結果を示す。図１６において、左側のKernel ABCの図は、シミュレーション実行部２２３およびKernel ABC 実行部２２４のみを１回だけ実行した場合の、各メッシュでの人口分布の変化を示している。真ん中のKernel Herding (1^st)の図は、シミュレーション実行部２２３、Kernel ABC 実行部２２４、カーネル平均構築部２２５、および第２のサンプル取得部２２６を１回だけ実行した場合の、各メッシュでの人口分布の変化を示している。右側のKernel Herding (2nd)の図は、シミュレーション実行部２２３、Kernel ABC 実行部２２４、カーネル平均構築部２２５、および第２のサンプル取得部２２６を２回実行した場合の、各メッシュでの人口分布の変化を示している。

　図１６から、Kernel ABCだけよりも（Kernel ABC + Kernel Herding）の方が、より正解データに近い人口分布を持つ推定結果が得られることが分かる。また、（Kernel ABC + Kernel Herding）を１回実行するよりも（Kernel ABC + Kernel Herding）を２回実行する方が、より正解データに近い人口分布を持つ推定結果が得られることが分かる。

　尚、上記実施の形態において、機械学習装置２２０の各部は、ハードウェアとソフトウェアとの組み合わせを用いて実現可能である。ハードウェアとソフトウェアとを組み合わせた形態では、記憶装置３００にプログラム３２０としてパラメータ推定プログラムが格納され、該パラメータ推定プログラムに基づいて機械学習装置２２０のハードウェアを動作させることによって、各部を各種手段として実現する。また、該パラメータ推定プログラムは、記録媒体に記録されて頒布されても良い。当該記録媒体に記録されたパラメータ推定プログラムは、有線、無線、又は記録媒体そのものを介して、メモリに読込まれ、処理部等を動作させる。尚、記録媒体を例示すれば、オプティカルディスクや磁気ディスク、半導体メモリ装置、ハードディスクなどが挙げられる。

　上記実施の形態を別の表現で説明すれば、機械学習装置２２０として動作させるコンピュータを、記憶装置３００にプログラム３２０として格納されたパラメータ推定プログラムに基づき、データ入力部２２１、第１のサンプル取得部２２２、シミュレーション実行部（関数実行部）２２３、Kernel ABC 実行部（類似度算出部）２２４、カーネル平均構築部２２５、第２のサンプル取得部２２６、およびサンプル評価部２２７として動作させることで実現することが可能である。

　以上、本発明の実施の形態および実施例を、図面を参照しつつ説明してきたが、当業者であれば、他の類似する実施形態および実施例を使用することができること、また、本発明から逸脱することなく適宜形態の変更又は追加を行うことができることに留意すべきである。

　上記の実施形態の一部又は全部は、以下の付記のようにも記載されうるが、以下には限られない。

（付記１）
　データｙを出力する関数のパラメータθを推定するパラメータ推定システムであって、
　前記関数に設定する前記パラメータθを、繰り返し学習により決定する機械学習装置を備え、
　前記機械学習装置は、
　事前分布Ｐから前記パラメータθのＮ（Ｎは２以上の整数）個のサンプルθ_ｉ（１≦ｉ≦Ｎ）を取得する第１のサンプル取得部と、
　前記Ｎ個のサンプルθ_ｉを順次前記関数に供給することで、それぞれ、前記関数からＮ個のデータｙ_ｉを取得する関数実行部と、
　カーネル関数ｋを用いて、前記Ｎ個のデータｙ_ｉとデータを表す正解データｙ^＊との間のＮ個の類似度ｗ_ｉを算出する類似度算出部と、
　該Ｎ個の類似度ｗ_ｉと、前記Ｎ個のサンプルθ_ｉと、前記カーネル関数ｋとから、事後分布のカーネル平均μを構築するカーネル平均構築部と、
　前記カーネル平均μと前記カーネル関数ｋとから、パラメータθのＴ（Ｔは２以上の整数）個の新たなサンプルθ_ｔ（１≦ｔ≦Ｔ）を取得する第２のサンプル取得部と、
　該Ｔ個の新たなサンプルθ_ｔから選択した１つのサンプルθ_ｋを前記関数に供給することによって前記関数から得られた新たなデータｙ_ｋと前記正解データｙ^＊との間の差分｜｜ｙ_ｋ－ｙ^＊｜｜が、所定の閾値εより小さいか否かを判断するサンプル評価部と、
を備え、
　前記サンプル評価部は、前記差分｜｜ｙ_ｋ－ｙ^＊｜｜が前記所定の閾値εより小さいと判断した場合に、前記選択したサンプルθ_ｋを、前記パラメータθとして推定し、出力することを特徴とするパラメータ推定システム。

（付記２）
　前記第１のサンプル取得部は、前記事前分布Ｐのハイパーパラメータηを調整しながら、前記Ｎ個のサンプルθ_ｉを取得する、付記１に記載のパラメータ推定システム。

（付記３）
　前記類似度算出部は、Kernel ABCアルゴリズムにより前記Ｎ個の類似度ｗ_ｉを算出するKernel ABC 実行部から成る、付記１又は２に記載のパラメータ推定システム。

（付記４）
　前記Kernel ABC 実行部は、前記カーネル関数ｋと、グラム行列Ｇと、正規化項δとを用いて、前記Ｎ個のデータｙ_ｉと前記正解データｙ^＊との間の前記Ｎ個の類似度ｗ_ｉを算出する、付記３に記載のパラメータ推定システム。

（付記５）
　前記カーネル平均構築部は、前記Ｎ個のサンプルθ_ｉを元に、前記Ｎ個の類似度ｗ_ｉと前記カーネル関数ｋとを使い重みを再調整し、確率分布を表す前記カーネル平均μを再生核ヒルベルト空間上で構築する、付記１乃至４のいずれか１項に記載のパラメータ推定システム。

（付記６）
　前記第２のサンプル取得部は、前記カーネル関数ｋ、前記新たなサンプルθ_ｔおよび前記取り出すサンプル数Ｔと、前記カーネル平均μとの間で求められる誤差を最小化する計算を行なって、前記新たなサンプルθ_ｔを取得する、付記１乃至５のいずれか１項に記載のパラメータ推定システム。

（付記７）
　前記サンプル評価部において、前記差分｜｜ｙ_ｋ－ｙ^＊｜｜が前記所定の閾値εより大きいと判断された場合に、前記シミュレーション実行部、前記類似度算出部、前記カーネル平均構築部、第２のサンプル取得部、および前記サンプル評価部の動作を繰り返す、付記１乃至６のいずれか１項に記載のパラメータ推定システム。

(付記８)
　前記サンプル評価部において、前記差分｜｜ｙ_ｋ－ｙ^＊｜｜が前記所定の閾値εより大きいと判断された場合に、前記第１のサンプル取得部、前記関数実行部、前記類似度算出部、前記カーネル平均構築部、第２のサンプル取得部、および前記サンプル評価部の動作を繰り返す、付記１乃至６のいずれか１項に記載のパラメータ推定システム。

(付記９)
　前記繰り返す動作を行なう前に、前記パラメータθの冗長性をコントロールしたい場合、前記第１のサンプル取得部は、前記事前分布Ｐのハイパーパラメータηを調整する、付記８に記載のパラメータ推定システム。

（付記１０）
　関数のパラメータθを推定するパラメータ推定システムにおける、前記パラメータθのサンプルから、データｙを出力する前記関数を用いて、前記関数に設定する前記パラメータθを、繰り返し学習により決定するパラメータ推定方法であって、
　第１のサンプル取得部が、事前分布Ｐから前記パラメータθのＮ（Ｎは２以上の整数）個のサンプルθ_ｉ（１≦ｉ≦Ｎ）を取得し、
　関数実行部が、前記Ｎ個のサンプルθ_ｉを順次前記関数に供給することで、それぞれ、前記関数からＮ個のデータｙ_ｉを取得し、
　類似度算出部が、カーネル関数ｋを用いて、前記Ｎ個のデータｙ_ｉとデータを表す正解データｙ^＊との間のＮ個の類似度ｗ_ｉを算出し、
　カーネル平均構築部が、前記Ｎ個の類似度ｗ_ｉと、前記Ｎ個のサンプルθ_ｉと、前記カーネル関数ｋとから、事後分布のカーネル平均μを構築し、
　第２のサンプル取得部が、前記カーネル平均μと前記カーネル関数ｋとから、パラメータθのＴ（Ｔは２以上の整数）個の新たなサンプルθ_ｔ（１≦ｔ≦Ｔ）を取得し、
　サンプル評価部が、前記Ｔ個の新たなサンプルθ_ｔから選択した１つのサンプルθ_ｋを前記関数に供給することによって前記関数から得られた新たなデータｙ_ｋと前記正解データｙ^＊との間の差分｜｜ｙ_ｋ－ｙ^＊｜｜が、所定の閾値εより小さいか否かを判断し、
　前記サンプル評価部は、前記差分｜｜ｙ_ｋ－ｙ^＊｜｜が前記所定の閾値εより小さいと判断した場合に、前記選択したサンプルθ_ｋを、前記パラメータθとして推定し、出力することを特徴とするパラメータ推定方法。

（付記１１）
　関数のパラメータθを推定するパラメータ推定システムにおいて、前記パラメータθのサンプルから、データｙを出力する前記関数を用いて、コンピュータに、前記関数に設定する前記パラメータθを、繰り返し学習により決定させるパラメータ推定プログラムを記録した記録媒体であって、前記パラメータ推定プログラムは、前記コンピュータに、
　事前分布Ｐから前記パラメータθのＮ（Ｎは２以上の整数）個のサンプルθ_ｉ（１≦ｉ≦Ｎ）を取得する第１のサンプル取得手順と、
　前記Ｎ個のサンプルθ_ｉを順次前記関数に供給することで、それぞれ、前記関数からＮ個のある分布からデータｙ_ｉを取得する関数実行手順と、
　カーネル関数ｋを用いて、前記Ｎ個のデータｙ_ｉとデータを表す正解データｙ^＊との間のＮ個の類似度ｗ_ｉを算出する類似度算出手順と、
　該Ｎ個の類似度ｗ_ｉと、前記Ｎ個のサンプルθ_ｉと、前記カーネル関数ｋとから、事後分布のカーネル平均μを構築するカーネル平均構築手順と、
　前記カーネル平均μと前記カーネル関数ｋとから、パラメータθのＴ（Ｔは２以上の整数）個の新たなサンプルθ_ｔ（１≦ｔ≦Ｔ）を取得する第２のサンプル取得手順と、
　該Ｔ個の新たなサンプルθ_ｔから選択した１つのサンプルθ_ｋを前記関数に供給することによって前記関数から得られた新たなデータｙ_ｋと前記正解データｙ^＊との間の差分｜｜ｙ_ｋ－ｙ^＊｜｜が、所定の閾値εより小さいか否かを判断するサンプル評価手順と、
を実行させ、
　前記サンプル評価手順は、前記コンピュータに、前記差分｜｜ｙ_ｋ－ｙ^＊｜｜が前記所定の閾値εより小さいと判断した場合に、前記選択したサンプルθ_ｋを、前記パラメータθとして推定させ、出力させることを特徴とするパラメータ推定プログラム記録媒体。

　本発明は、関数が評価できない場合に最尤法によりパラメータを推定する分野に適用可能である。例えば、本発明は、金融時系列の分散を表すモデルのパラメータや、生物の個体数の時間発展を記述するモデルのパラメータ、遺伝的多様性を実現するようなモデルのパラメータを推定する分野にも適用可能である。

　この出願は、２０１７年１０月３日に出願された日本出願特願２０１７－１９３２４５を基礎とする優先権を主張し、その開示の全てをここに取り込む。

　１００　　人流パターン推定システム（パラメータ推定システム）
　２００　　データ処理装置
　２１０　　人流シミュレータ（関数）
　２２０　　機械学習装置
　２２１　　データ入力部
　２２２　　第１のサンプル取得部（事前分布からのサンプル取得部）
　２２３　　シミュレーション実行部（関数実行部）
　２２４　　Kernel ABC 実行部（類似度算出部）
　２２５　　カーネル平均構築部
　２２６　　第２のサンプル取得部（事後分布からのサンプル取得部）
　２２７　　サンプル評価部
　３００　　記憶装置
　３１０　　処理情報
　３２０　　プログラム
　４００　　入力装置
　５００　　出力装置

Claims

　データｙを出力する関数のパラメータθを推定するパラメータ推定システムであって、
　前記関数に設定する前記パラメータθを、繰り返し学習により決定する機械学習装置を備え、
　前記機械学習装置は、
　事前分布Ｐから前記パラメータθのＮ（Ｎは２以上の整数）個のサンプルθ_ｉ（１≦ｉ≦Ｎ）を取得する第１のサンプル取得部と、
　前記Ｎ個のサンプルθ_ｉを順次前記関数に供給することで、それぞれ、前記関数からＮ個のデータｙ_ｉを取得する関数実行部と、
　カーネル関数ｋを用いて、前記Ｎ個のデータｙ_ｉとデータを表す正解データｙ^＊との間のＮ個の類似度ｗ_ｉを算出する類似度算出部と、
　該Ｎ個の類似度ｗ_ｉと、前記Ｎ個のサンプルθ_ｉと、前記カーネル関数ｋとから、事後分布のカーネル平均μを構築するカーネル平均構築部と、
　前記カーネル平均μと前記カーネル関数ｋとから、パラメータθのＴ（Ｔは２以上の整数）個の新たなサンプルθ_ｔ（１≦ｔ≦Ｔ）を取得する第２のサンプル取得部と、
　該Ｔ個の新たなサンプルθ_ｔから選択した１つのサンプルθ_ｋを前記関数に供給することによって前記関数から得られた新たなデータｙ_ｋと前記正解データｙ^＊との間の差分｜｜ｙ_ｋ－ｙ^＊｜｜が、所定の閾値εより小さいか否かを判断するサンプル評価部と、
を備え、
　前記サンプル評価部は、前記差分｜｜ｙ_ｋ－ｙ^＊｜｜が前記所定の閾値εより小さいと判断した場合に、前記選択したサンプルθ_ｋを、前記パラメータθとして推定し、出力することを特徴とするパラメータ推定システム。
　前記第１のサンプル取得部は、前記事前分布Ｐのハイパーパラメータηを調整しながら、前記Ｎ個のサンプルθ_ｉを取得する、請求項１に記載のパラメータ推定システム。
　前記類似度算出部は、Kernel ABCアルゴリズムにより前記Ｎ個の類似度ｗ_ｉを算出するKernel ABC 実行部から成る、請求項１又は２に記載のパラメータ推定システム。
　前記Kernel ABC 実行部は、前記カーネル関数ｋと、グラム行列Ｇと、正規化項δとを用いて、前記Ｎ個のデータｙ_ｉと前記正解データｙ^＊との間の前記Ｎ個の類似度ｗ_ｉを算出する、請求項３に記載のパラメータ推定システム。
　前記カーネル平均構築部は、前記Ｎ個のサンプルθ_ｉを元に、前記Ｎ個の類似度ｗ_ｉと前記カーネル関数ｋとを使い重みを再調整し、確率分布を表す前記カーネル平均μを再生核ヒルベルト空間上で構築する、請求項１乃至４のいずれか１項に記載のパラメータ推定システム。
　前記第２のサンプル取得部は、前記カーネル関数ｋ、前記新たなサンプルθ_ｔおよび前記取り出すサンプル数Ｔと、前記カーネル平均μとの間で求められる誤差を最小化する計算を行なって、前記新たなサンプルθ_ｔを取得する、請求項１乃至５のいずれか１項に記載のパラメータ推定システム。
　前記サンプル評価部において、前記差分｜｜ｙ_ｋ－ｙ^＊｜｜が前記所定の閾値εより大きいと判断された場合に、前記シミュレーション実行部、前記類似度算出部、前記カーネル平均構築部、第２のサンプル取得部、および前記サンプル評価部の動作を繰り返す、請求項１乃至６のいずれか１項に記載のパラメータ推定システム。
　前記サンプル評価部において、前記差分｜｜ｙ_ｋ－ｙ^＊｜｜が前記所定の閾値εより大きいと判断された場合に、前記第１のサンプル取得部、前記関数実行部、前記類似度算出部、前記カーネル平均構築部、第２のサンプル取得部、および前記サンプル評価部の動作を繰り返す、請求項１乃至６のいずれか１項に記載のパラメータ推定システム。
　前記繰り返す動作を行なう前に、前記パラメータθの冗長性をコントロールしたい場合、前記第１のサンプル取得部は、前記事前分布Ｐのハイパーパラメータηを調整する、請求項８に記載のパラメータ推定システム。
　関数のパラメータθを推定するパラメータ推定システムにおける、前記パラメータθのサンプルから、データｙを出力する前記関数を用いて、前記関数に設定する前記パラメータθを、繰り返し学習により決定するパラメータ推定方法であって、
　第１のサンプル取得部が、事前分布Ｐから前記パラメータθのＮ（Ｎは２以上の整数）個のサンプルθ_ｉ（１≦ｉ≦Ｎ）を取得し、
　関数実行部が、前記Ｎ個のサンプルθ_ｉを順次前記関数に供給することで、それぞれ、前記関数からＮ個のデータｙ_ｉを取得し、
　類似度算出部が、カーネル関数ｋを用いて、前記Ｎ個のデータｙ_ｉとデータを表す正解データｙ^＊との間のＮ個の類似度ｗ_ｉを算出し、
　カーネル平均構築部が、前記Ｎ個の類似度ｗ_ｉと、前記Ｎ個のサンプルθ_ｉと、前記カーネル関数ｋとから、事後分布のカーネル平均μを構築し、
　第２のサンプル取得部が、前記カーネル平均μと前記カーネル関数ｋとから、パラメータθのＴ（Ｔは２以上の整数）個の新たなサンプルθ_ｔ（１≦ｔ≦Ｔ）を取得し、
　サンプル評価部が、前記Ｔ個の新たなサンプルθ_ｔから選択した１つのサンプルθ_ｋを前記関数に供給することによって前記関数から得られた新たなデータｙ_ｋと前記正解データｙ^＊との間の差分｜｜ｙ_ｋ－ｙ^＊｜｜が、所定の閾値εより小さいか否かを判断し、
　前記サンプル評価部は、前記差分｜｜ｙ_ｋ－ｙ^＊｜｜が前記所定の閾値εより小さいと判断した場合に、前記選択したサンプルθ_ｋを、前記パラメータθとして推定し、出力することを特徴とするパラメータ推定方法。
　関数のパラメータθを推定するパラメータ推定システムにおいて、前記パラメータθのサンプルから、データｙを出力する前記関数を用いて、コンピュータに、前記関数に設定する前記パラメータθを、繰り返し学習により決定させるパラメータ推定プログラムを記録した記録媒体であって、前記パラメータ推定プログラムは、前記コンピュータに、
　事前分布Ｐから前記パラメータθのＮ（Ｎは２以上の整数）個のサンプルθ_ｉ（１≦ｉ≦Ｎ）を取得する第１のサンプル取得手順と、
　前記Ｎ個のサンプルθ_ｉを順次前記関数に供給することで、それぞれ、前記関数からＮ個のある分布からデータｙ_ｉを取得する関数実行手順と、
　カーネル関数ｋを用いて、前記Ｎ個のデータｙ_ｉとデータを表す正解データｙ^＊との間のＮ個の類似度ｗ_ｉを算出する類似度算出手順と、
　該Ｎ個の類似度ｗ_ｉと、前記Ｎ個のサンプルθ_ｉと、前記カーネル関数ｋとから、事後分布のカーネル平均μを構築するカーネル平均構築手順と、
　前記カーネル平均μと前記カーネル関数ｋとから、パラメータθのＴ（Ｔは２以上の整数）個の新たなサンプルθ_ｔ（１≦ｔ≦Ｔ）を取得する第２のサンプル取得手順と、
　該Ｔ個の新たなサンプルθ_ｔから選択した１つのサンプルθ_ｋを前記関数に供給することによって前記関数から得られた新たなデータｙ_ｋと前記正解データｙ^＊との間の差分｜｜ｙ_ｋ－ｙ^＊｜｜が、所定の閾値εより小さいか否かを判断するサンプル評価手順と、
を実行させ、
　前記サンプル評価手順は、前記コンピュータに、前記差分｜｜ｙ_ｋ－ｙ^＊｜｜が前記所定の閾値εより小さいと判断した場合に、前記選択したサンプルθ_ｋを、前記パラメータθとして推定させ、出力させることを特徴とするパラメータ推定プログラム記録媒体。