(57)【特許請求の範囲】
【請求項１】
個々を識別可能な第１の個体群に含まれる個体ｉ（１≦ｉ≦Ｉ，Ｉは１以上の整数）とオブジェクトｊ（１≦ｊ≦Ｊ、Ｊは１以上の整数）との関連度を表すと共に、平均パラメータμｉｊ及び分散パラメータを持つ対数正規分布に従う要素ｘｉｊｋからなるベクトルｘｉｊを要素として持つＩ×Ｊのベクトル値行列Ｘから、
前記個体ｉが、クラスタｒ（１≦ｒ≦Ｒ、Ｒは１以上の整数）に所属することを表す非負値の要素ａｉｒを持つＩ×Ｒの第１の特徴行列Ａと、前記オブジェクトｊが、前記クラスタｒに所属することを表す非負値の要素ｂｊｒを持つＪ×Ｒの第２の特徴行列Ｂと、を抽出する分布推定装置であって、
前記ベクトル値行列Ｘの各ベクトルｘｉｊの各要素ｘｉｊｋの対数値、前記第１の特徴行列Ａ、前記第２の特徴行列Ｂ、及び前記分散パラメータを用いて表される目的関数を最適化するように、前記第１の特徴行列Ａ及び前記第２の特徴行列Ｂを推定する特徴行列推定部と、
予め定めた反復終了条件を満足するまで、前記特徴行列推定部による推定を繰り返す反復判定部と、
を備えた分布推定装置。
【請求項２】
前記ベクトルｘｉｊの要素ｘｉｊｋは、平均パラメータμｉｊ及び個体ｉに依存する分散パラメータσｉを持つ対数正規分布、又は、平均パラメータμｉｊ及びオブジェクトｊに依存する分散パラメータσｊを持つ対数正規分布に従い、
前記特徴行列推定部は、前記ベクトル値行列Ｘの各ベクトルｘｉｊの各要素ｘｉｊｋの対数値、前記第１の特徴行列Ａ、前記第２の特徴行列Ｂ、及び前記分散パラメータσｉ又は前記分散パラメータσｊを用いて表される目的関数を最適化するように、前記第１の特徴行列Ａ及び前記第２の特徴行列Ｂを推定する
請求項１記載の分布推定装置。
【請求項３】
前記ベクトル値行列Ｘは、時間帯ｉにおける各ユーザの場所ｊの滞在時間を表す要素ｘｉｊｋからなるベクトルｘｉｊを要素として持ち、
前記第１の特徴行列Ａは、前記時間帯ｉが、前記クラスタｒに所属することを表す非負値の要素ａｉｒを持ち、
前記第２の特徴行列Ｂは、前記場所ｊが、前記クラスタｒに所属することを表す非負値の要素ｂｊｒを持つ
請求項１記載の分布推定装置。
【請求項４】
前記ベクトルｘｉｊの要素ｘｉｊｋは、平均パラメータμｉｊ及び場所ｊに依存する分散パラメータσｊを持つ対数正規分布に従い、
前記特徴行列推定部は、下記の式で表される前記目的関数を最小化するように、前記第１の特徴行列Ａ及び前記第２の特徴行列Ｂを推定する
【数１】
0006469033.tif 000036
請求項３記載の分布推定装置。
【請求項５】
前記ベクトルｘｉｊの要素ｘｉｊｋは、平均パラメータμｉｊ、及び時間帯ｉに依存する分散パラメータσｉを持つ対数正規分布に従い、
前記特徴行列推定部は、下記の式で表される前記目的関数を最小化するように、前記第１の特徴行列Ａ及び前記第２の特徴行列Ｂを推定する
【数２】
0006469033.tif 000037
請求項３記載の分布推定装置。
【請求項６】
個々を識別可能な第１の個体群に含まれる個体ｉ（１≦ｉ≦Ｉ，Ｉは１以上の整数）とオブジェクトｊ（１≦ｊ≦Ｊ、Ｊは１以上の整数）との関連度を表すと共に、平均パラメータμｉｊ及び分散パラメータを持つ対数正規分布に従う要素ｘｉｊｋからなるベクトルｘｉｊを要素として持つＩ×Ｊのベクトル値行列Ｘから、
前記個体ｉが、クラスタｒ（１≦ｒ≦Ｒ、Ｒは１以上の整数）に所属することを表す非負値の要素ａｉｒを持つＩ×Ｒの第１の特徴行列Ａと、前記オブジェクトｊが、前記クラスタｒに所属することを表す非負値の要素ｂｊｒを持つＪ×Ｒの第２の特徴行列Ｂと、を抽出する分布推定装置における分布推定方法であって、
特徴行列推定部が、前記ベクトル値行列Ｘの各ベクトルｘｉｊの各要素ｘｉｊｋの対数値、前記第１の特徴行列Ａ、前記第２の特徴行列Ｂ、及び前記分散パラメータを用いて表される目的関数を最適化するように、前記第１の特徴行列Ａ及び前記第２の特徴行列Ｂを推定する特徴行列推定ステップと、
反復判定部が、予め定めた反復終了条件を満足するまで、前記特徴行列推定ステップによる推定を繰り返す反復判定ステップと、
を行う分布推定方法。
【請求項７】
前記ベクトル値行列Ｘは、時間帯ｉにおける各ユーザの場所ｊの滞在時間を表す要素ｘｉｊｋからなるベクトルｘｉｊを要素として持ち、
前記第１の特徴行列Ａは、前記時間帯ｉが、前記クラスタｒに所属することを表す非負値の要素ａｉｒを持ち、
前記第２の特徴行列Ｂは、前記場所ｊが、前記クラスタｒに所属することを表す非負値の要素ｂｊｒを持つ
請求項６記載の分布推定方法。
【請求項８】
コンピュータを、請求項１〜５の何れか１項記載の分布推定装置の各部として機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、行列形式で与えられるデータから、データが従う分布を推定する分布推定装置、分布推定方法、及び分布推定プログラムに係り、特に、ユーザの特定の場所における滞在時間のデータから、データが従う分布を推定する分布推定装置、分布推定方法、及び分布推定プログラムに関する。
【背景技術】
【０００２】
ＰＯＳ（Point of Sales）データに代表される購買履歴等の構造化されたデータ、テキストデータ・画像データ等の構造化されていないデータ等の多くは、前処理によって実数値を要素に持つ行列形式で表現できることが知られている。このように行列表現されたデータ中に存在するクラスタを発見する手法として、非負値行列分解（Non-negative Matrix Factorization, NMF）と呼ばれる手法の有用性がこれまで示されている（例えば、非特許文献１を参照）。
【０００３】
ＮＭＦを適用する際に入力される行列データは、当該行列データより低次の階数の行列の積に分解される。この各々の低次の階数の行列が、それぞれ各行、各列に対応する事物のクラスタへの寄与度を表しており、この寄与度によりクラスタ発見が可能となる。従って、例えば、ユーザの訪問履歴に関するデータを行列形式により表現することにより、おすすめのショップのリストを作成すること等が可能となる。
【０００４】
図１０に、ユーザの訪問履歴Ｘに関するデータをＮＭＦに適用した適応例を示す。図１０に示すように、ユーザの場所毎の訪問回数を表すユーザ訪問場所行列Ｘは、行列中の各行ｉがユーザを表し、各列ｊが訪問場所を表し、各値が訪問回数を表すＩ行Ｊ列の行列である。図１０に示す例では、１行目に対応するユーザがユーザ１であり、１列目に対応する訪問場所が場所１であり、ユーザ１の場所１への訪問回数が４である。
【０００５】
ユーザ訪問場所行列Ｘにおいては、訪問回数の値が大きいほど、その場所が人気のスポットであることを表す。このように、ユーザ訪問場所行列ＸをＮＭＦに適用することで、
【数１】
0006469033.tif 000002
となるＩ行Ｒ列のユーザ特徴行列Ａ＝｛ａｉｒ｝と、Ｊ行Ｒ列の訪問場所特徴行列Ｂ＝｛ｂｊｒ｝が求まる。ただし、記号の上付きの記号Ｔは行列の転置を表す。
【０００６】
また、ここで、記号
【数２】
0006469033.tif 000003
で表現した類似の尺度について説明する。上記非特許文献１にも記述されているように、行列の類似の尺度には、ユークリッド距離に基づくものや一般化カルバックライブラーダイバージェンス（ＫＬ距離）により定義される距離尺度が用いられ、値が小さいほど両者が類似していることを表す。
【０００７】
図１０において、ユーザ特徴行列Ａのクラスタ１に対応する１列目に着目すると、ユーザ１に対応する１行目、ユーザ２に対応する２行目、及びユーザ３に対応する３行目の値がそれぞれ０より大きい値となっている。これは、ユーザ１、ユーザ２、及びユーザ３がクラスタ１に所属することを示している。
【０００８】
また、図１０において、訪問場所特徴行列Ｂに着目すると、場所１に対応する１列目、場所２に対応する２列目、及び場所３に対応する３列目では、クラスタ１に対応する１行目の値がクラスタ２に対応する２行目の値より大きい値となっている。これは、場所１、場所２、及び場所３は、ユーザ１、ユーザ２、及びユーザ３が訪れやすい場所であるというクラスタ１が持つ特徴を表している。
【０００９】
これらを踏まえ、図１１に示すように、場所１、場所２、及び場所３をまとめてクラスタ１の場所特徴とする。同様に、クラスタ１に所属するユーザ１、ユーザ２、及びユーザ３をクラスタ１のユーザ特徴とする。以下、クラスタ１の場所特徴及びユーザ特徴をまとめてクラスタ１の特徴ともいう。
【００１０】
このように、ＮＭＦの適用によって得られたユーザ特徴行列Ａ、及び訪問場所特徴行列Ｂに基づき、図１１に示すようなクラスタ抽出が可能となる。
【００１１】
なお、クラスタの総数に相当する訪問場所特徴行列Ｂの階数は、解析する前に予め決定しておくものとする。
【００１２】
ＮＭＦでは、行列を分解する際に、行列に対して類似の尺度を定義しており、定義した尺度のもとで値が小さいほど両者が類似していることを表す。従って、ＮＭＦでは、採用した類似の尺度を最小化する行列（例えば、ユーザ特徴行列Ａ、及び訪問場所特徴行列Ｂを求める手法として定式化される。
【００１３】
また、ＮＭＦでは、データが持つ性質を考慮して、利用する距離が決定される。例えば、距離としてＫＬダイバージェンスが採用される場合は、行列の各要素ｘｉｊは平均μ＝ｘ＾ｉｊのポアソン分布Ｐｏ（ｘｉｊ｜μｉｊ）に従っていると仮定していることに相当する。ポアソン分布Ｐｏ（ｘｉｊ｜μｉｊ）に従っているデータとしては、上述したユーザの場所毎の訪問回数等が例示される。
【００１４】
なお、行列Ｘと行列Ｘ＾のＫＬダイバージェンスＤＫＬは、下記（１）式で定義される。
【００１５】
【数３】
0006469033.tif 000004
【００１６】
実数値が従う確率分布としては正規分布が有用であり、頻度を表す離散値が従う確率分布としてはポアソン分布が有用であることは広く認識された事実である。
【先行技術文献】
【非特許文献】
【００１７】
【非特許文献１】澤田宏, “非負値行列因子分解NMFの基礎とデータ／信号解析への応用”, 電子情報通信学会誌, Vol. 95, No. 9, pp. 829-833, 2012.
【非特許文献２】Kyosuke Nishida, Hiroyuki Toda, Takeshi Kurashima, and Yoshihiko Suhara, "Probabilistic Identification of Visited Point-of-Interest for Personalized Automatic Check-in", pp.631-642, in Proceedings of the 2014 ACM International Joint Conference on Pervasive and Ubiquitous Computing (UbiComp 2014).
【発明の概要】
【発明が解決しようとする課題】
【００１８】
従来、ある時間帯にある場所を訪れたユーザがどのくらいの時間その場所に滞在したかを表す入力データから、滞在時間の従う確率分布の推定を行う技術が用いられている。この場合の入力データは、図１２の左図に示すように、各行ｉが時間帯を表し、各列ｊが場所を表し、要素としてベクトルｘｉｊ＝（ｘｉｊ１，ｘｉｊ２，…，ｘｉｊＫｉｊ）を持つＩ行Ｊ列のベクトル値行列Ｘとして表現される。ただし、Ｋｉｊは時間帯ｉに場所ｊを訪れたユーザの総数であり、ｘｉｊｋが時間帯ｉに場所ｊにｋ番目に訪れた人の滞在時間を表す。以後、このベクトル値行列Ｘを時間帯場所ベクトル値行列という。
【００１９】
なお、ここでいう時間帯ｉは、例えばｉ＝１が［７：００−８：００］の間の１時間、ｉ＝２が［８：００−９：００］の１時間のように、２４時間をなんらかの基準で分割したものであってもよいし、例えばｉ＝１が休日の午前中、ｉ＝２が休日の午後、ｉ＝３が平日の午前中、ｉ＝４が平日の午後のように、平日及び休日といった日に関する情報と、時間に関する情報とを組み合わせて作成したものであってもよい。
【００２０】
また、場所ｊは、喫茶店、レストラン等の店舗に対応させても良く、地図を何らかの基準でグリッド化したグリッドのＩＤに対応させても良い。
【００２１】
このような入力データから、データが従う分布の推定を行う場合、時間帯場所ベクトル値行列Ｘを、通常の要素がスカラである行列に変換する平均化処理を適用した上で、既存のＮＭＦを適用することが考えられる。より正確に述べると、時間帯場所ベクトル値行列Ｘにおいてユーザ毎に平均を取った時間帯場所平均行列
【数４】
0006469033.tif 000005
を下記（２）式に従って作成し、時間帯場所平均行列ＸにＮＭＦを適用する。
【００２２】
【数５】
0006469033.tif 000006
【００２３】
しかし、この手法では、各要素の分散の情報を潰しており、時間帯ｉに場所ｊに滞在した滞在時間の平均値のみに基づいて分布の推定を行うことになる。すなわち、この手法では、以下の（Ａ）及び（Ｂ）により、滞在時間が従う確率分布を推定したい場合には不適切だといえる。
【００２４】
（Ａ）上記手法は、正規分布から得られたと仮定されるデータに対して推定を行うためである。上記非特許文献２でも示されるように、ユーザがある場所に滞在した滞在時間の分布は対数正規分布に従うことが知られている。つまり、上記手法により分布の推定を行うと、実際には図１２の右図に示すように、ユーザがある場所に滞在した滞在時間の分布は正規分布とは限らない。しかし、図１３に示すように、ユーザがある場所に滞在した滞在時間の分布を正規分布として推定を行ってしまい、滞在時間が従う分布の推定を精度良く行うことができない。
【００２５】
（Ｂ）上記手法では、平均値を推定しているため、分散値も考慮した分布の形を推定する場合には不適切である。
【００２６】
本発明は、以上のような事情に鑑みてなされたものであり、要素の分散を考慮して、分布を推定することができる分布推定装置、分布推定方法、及び分布推定プログラムを提供することを目的とする。
【課題を解決するための手段】
【００２７】
上記目的を達成するために、本発明の分布推定装置は、個々を識別可能な第１の個体群に含まれる個体ｉ（１≦ｉ≦Ｉ，Ｉは１以上の整数）とオブジェクトｊ（１≦ｊ≦Ｊ、Ｊは１以上の整数）との関連度を表すと共に、平均パラメータμｉｊ及び分散パラメータを持つ対数正規分布に従う要素ｘｉｊｋからなるベクトルｘｉｊを要素として持つＩ×Ｊのベクトル値行列Ｘから、前記個体ｉが、クラスタｒ（１≦ｒ≦Ｒ、Ｒは１以上の整数）に所属することを表す非負値の要素ａｉｒを持つＩ×Ｒの第１の特徴行列Ａと、前記オブジェクトｊが、前記クラスタｒに所属することを表す非負値の要素ｂｊｒを持つＪ×Ｒの第２の特徴行列Ｂと、を抽出する分布推定装置であって、前記ベクトル値行列Ｘの各ベクトルｘｉｊの各要素ｘｉｊｋの対数値、前記第１の特徴行列Ａ、前記第２の特徴行列Ｂ、及び前記分散パラメータを用いて表される目的関数を最適化するように、前記第１の特徴行列Ａ及び前記第２の特徴行列Ｂを推定する特徴行列推定部と、予め定めた反復終了条件を満足するまで、前記特徴行列推定部による推定を繰り返す反復判定部と、を備える。
【００２８】
なお、本発明のｖ装置において、前記ベクトルｘｉｊの要素ｘｉｊｋは、平均パラメータμｉｊ及び個体ｉに依存する分散パラメータσｉを持つ対数正規分布、又は、平均パラメータμｉｊ及びオブジェクトｊに依存する分散パラメータσｊを持つ対数正規分布に従い、前記特徴行列推定部は、前記ベクトル値行列Ｘの各ベクトルｘｉｊの各要素ｘｉｊｋの対数値、前記第１の特徴行列Ａ、前記第２の特徴行列Ｂ、及び前記分散パラメータσｉ又は前記分散パラメータσｊを用いて表される目的関数を最適化するように、前記第１の特徴行列Ａ及び前記第２の特徴行列Ｂを推定するようにしても良い。
【００２９】
また、本発明の分布推定装置において、前記ベクトル値行列Ｘは、時間帯ｉにおける各ユーザの場所ｊの滞在時間を表す要素ｘｉｊｋからなるベクトルｘｉｊを要素として持ち、前記第１の特徴行列Ａは、前記時間帯ｉが、前記クラスタｒに所属することを表す非負値の要素ａｉｒを持ち、前記第２の特徴行列Ｂは、前記場所ｊが、前記クラスタｒに所属することを表す非負値の要素ｂｊｒを持つようにしても良い。
【００３０】
また、本発明の分布推定装置において、前記ベクトルｘｉｊの要素ｘｉｊｋは、平均パラメータμｉｊ及び場所ｊに依存する分散パラメータσｊを持つ対数正規分布に従い、前記特徴行列推定部は、前記ベクトル値行列Ｘの各ベクトルｘｉｊの各要素ｘｉｊｋの対数値、前記第１の特徴行列Ａ、前記第２の特徴行列Ｂ、及び前記分散パラメータσｊを用いて表される目的関数を最小化するように、前記第１の特徴行列Ａ及び前記第２の特徴行列Ｂを推定するようにしても良い。
【００３１】
また、本発明の分布推定装置において、前記ベクトルｘｉｊの要素ｘｉｊｋは、平均パラメータμｉｊ、及び時間帯ｉに依存する分散パラメータσｉを持つ対数正規分布に従い、前記特徴行列推定部は、前記ベクトル値行列Ｘの各ベクトルｘｉｊの各要素ｘｉｊｋの対数値、前記第１の特徴行列Ａ、前記第２の特徴行列Ｂ、及び前記分散パラメータσｊを用いて表される目的関数を最小化するように、前記第１の特徴行列Ａ及び前記第２の特徴行列Ｂを推定するようにしても良い。
【００３２】
上記目的を達成するために、本発明の分布推定方法は、個々を識別可能な第１の個体群に含まれる個体ｉ（１≦ｉ≦Ｉ，Ｉは１以上の整数）とオブジェクトｊ（１≦ｊ≦Ｊ、Ｊは１以上の整数）との関連度を表すと共に、平均パラメータμｉｊ及び分散パラメータを持つ対数正規分布に従う要素ｘｉｊｋからなるベクトルｘｉｊを要素として持つＩ×Ｊのベクトル値行列Ｘから、前記個体ｉが、クラスタｒ（１≦ｒ≦Ｒ、Ｒは１以上の整数）に所属することを表す非負値の要素ａｉｒを持つＩ×Ｒの第１の特徴行列Ａと、前記オブジェクトｊが、前記クラスタｒに所属することを表す非負値の要素ｂｊｒを持つＪ×Ｒの第２の特徴行列Ｂと、を抽出する分布推定装置における分布推定方法であって、特徴行列推定部が、前記ベクトル値行列Ｘの各ベクトルｘｉｊの各要素ｘｉｊｋの対数値、前記第１の特徴行列Ａ、前記第２の特徴行列Ｂ、及び前記分散パラメータを用いて表される目的関数を最適化するように、前記第１の特徴行列Ａ及び前記第２の特徴行列Ｂを推定する特徴行列推定ステップと、反復判定部が、予め定めた反復終了条件を満足するまで、前記特徴行列推定ステップによる推定を繰り返す反復判定ステップと、を行う。
【００３３】
また、本発明の分布推定方法において、前記ベクトル値行列Ｘは、時間帯ｉにおける各ユーザの場所ｊの滞在時間を表す要素ｘｉｊｋからなるベクトルｘｉｊを要素として持ち、前記第１の特徴行列Ａは、前記時間帯ｉが、前記クラスタｒに所属することを表す非負値の要素ａｉｒを持ち、前記第２の特徴行列Ｂは、前記場所ｊが、前記クラスタｒに所属することを表す非負値の要素ｂｊｒを持つようにしても良い。
【００３４】
上記目的を達成するために、本発明の文書分類プログラムは、コンピュータを、上記分布推定装置の各部として機能させるためのプログラムである。
【発明の効果】
【００３５】
本発明によれば、要素の分散を考慮して、分布を推定することができる。
【図面の簡単な説明】
【００３６】
【図１】実施形態に係る分布推定装置における場所及び滞在時間の分布を推定する推定方法の一例を示す機能ブロック図である。
【図２】対数正規分布の一例を示す模式図である。
【図３】実施形態に係る分布推定装置における時間帯場所ベクトル値行列を時間帯場所ベクトル値対行列に変換する変換方法の一例を示す模式図である。
【図４】実施形態に係る分布推定装置における時間帯場所ベクトル値対数行列を時間帯場所対数平均値ベクトルに変換する変換方法の一例を示す模式図である。
【図５】実施形態に係る分布推定装置における時間帯場所行列及び場所特徴行列を推定する推定方法の一例を示す模式図である。
【図６】実施形態に係る分布推定装置における時間帯場所ベクトル値行列の平均値を推定する推定方法の一例を示す模式図である。
【図７】実施形態に係る分布推定装置の構成の一例を示す機能ブロック図である。
【図８】実施形態に係る分布推定装置により実行される全体処理の流れを示すフローチャートである。
【図９】実施形態に係る分布推定装置により実行される特徴行列推定処理の流れを示すフローチャートである。
【図１０】ユーザの訪問履歴に関するデータをＮＭＦに適用した行列分解の一例を示す模式図である。
【図１１】ユーザの訪問履歴に関するデータをＮＭＦに適用した行列分解によりクラスタを抽出方法の一例を示す模式図である。
【図１２】時間帯場所ベクトル値行列の一例を示す模式図である。
【図１３】場所及び対外時間の分布を推定する推定方法の一例である。
【発明を実施するための形態】
【００３７】
以下、本発明の実施形態について図面を用いて説明する。
【００３８】
本実施形態に係る分布推定装置は、時間帯場所ベクトル値行列Ｘにおける滞在時間が従う分布の推定を行う非負値行列分解を用いる。特に、本実施形態では、図１に示すように、要素である滞在時間を対数正規分布に従って得られると仮定して、ユーザの滞在時間の分布の推定を行う。
【００３９】
対数正規分布は、ユーザの滞在時間等のモデリングにおいて広く利用される確率分布である（上記非特許文献２を参照）。対数正規分布は、正規分布とは異なり、図２に示すように、右に裾の長いデータを表現することが可能な確率分布である。対数正規分布は、この性質によって、例えばカフェ等における『話し込んでしまい長居をしてしまう』（滞在時間が長くなる）ことを表現することが可能となる。
【００４０】
本実施形態では、時間帯場所ベクトル値行列Ｘにおける各要素ｘｉｊｋが、平均パラメータμｉｊ、分散パラメータσｉｊを持つ対数正規分布ＬＮ（ｘｉｊｋ｜μｉｊ，σｉｊ）に従い、平均パラメータμｉｊは、時間帯特徴行列Ａ、場所特徴行列Ｂの要素を用いて下記（３）式のように表現されるとしてモデリングを行う。
【００４１】
【数６】
0006469033.tif 000007
【００４２】
時間帯特徴行列Ａ、場所特徴行列Ｂ、及び分散パラメータσｉｊは、下記（４）式に示す目的関数を用いて最小化を行うことにより推定できる。
【００４３】
【数７】
0006469033.tif 000008
【００４４】
このように、本実施形態では、平均パラメータμｉｊが、時間帯特徴行列Ａと場所特徴行列Ｂとの積で表現される。これは一見、入力行列である時間帯場所ベクトル値行列Ｘの行列分解を行う通常のＮＭＦとは全く異なる手法に見える。しかしながら、本実施形態の技術とＮＭＦには強いつながりがある。それは、対数正規分布と正規分布との関係、すなわち「対数正規分布ＬＮ（ｘ｜μ，σ）に従う変数Ｘの対数関数をとった場合、新たな確率変数Ｙ＝ｌｏｇ（Ｘ）が正規分布Ｎ（ｙ｜μ，σ）に従う」という関係に由来するものである。
【００４５】
まず、図３に示すように、時間帯場所ベクトル値行列Ｘを対数変換した時間帯場所ベクトル値対数行列Ｙを算出する。時間帯場所ベクトル値対数行列Ｙのｙｉｊの要素には、訪問したＫｉｊ人のユーザの訪問時間情報（対数値）が、時間帯場所ベクトル値行列Ｘと同様にＫｉｊ次元ベクトル要素として格納されている。
【００４６】
また、図４に示すように、時間帯場所ベクトル値対数行列Ｙを、時間帯場所平均対数値行列
【数８】
0006469033.tif 000009
に変換する。時間帯場所平均対数値行列
【数９】
0006469033.tif 000010
の
【数１０】
0006469033.tif 000011
の要素には、訪問したＫｉｊ人のユーザの平均訪問時間情報（対数平均値）が格納されている。なお、この場合の
【数１１】
0006469033.tif 000012
の要素は、ベクトルではなくスカラである。
【００４７】
本実施形態によるＮＭＦは、図５に示すように、ｙｉｊとの誤差を上記（４）式に示す目的関数に従って最小化するように、ＮＭＦによって因子行列である時間帯特徴行列Ａ及び場所特徴行列Ｂを導出する。これらの因子行列によって再構築されるμｉｊの要素は、上記（４）式に従って、分散パラメータσｉｊも小さくするようにパラメータが抽出される。なお、分散パラメータσｉｊは、ｙｉｊｋが平均パラメータμｉｊからどれだけばらつくかという分散を表したパラメータとなっている。このように、本実施形態では、ＮＭＦの拡張技術であるとみなすことができる。
【００４８】
さらに、上記（４）式では、データ数Ｋｉｊの和を取っていることがわかる。このことから、データ数Ｋｉｊが多ければ多いほど、目的関数の値が大きくなってしまう。そのため、目的関数を小さくするためには分散パラメータσｉｊを小さくすることになる。すなわち、データ数Ｋｉｊが多い場所に関する分布を推定する場合ほど、分散パラメータσｉｊが小さくなり、正確な推定を行えるような目的関数となっている。
【００４９】
このように、本実施形態では、時間帯場所ベクトル値行列Ｘを入力行列とした場合の、各場所及び各時間に対するユーザの滞在時間が従う分布の推定を行うことで、滞在時間の推定が行える。また、対数正規分布という滞在時間のモデリングに適した分布を利用していることから、ユーザの滞在時間の推定精度の向上も期待できる。これは、図６に示すように、時間帯ｉに場所ｊに滞在する時間を表すベクトルｘｉｊを生成する確率が最も高い対数正規分布のパラメータを推定していることに等しい。パラメータ推定結果を用いることで、時間帯及び場所毎の平均滞在時間を対数正規分布の平均値
【数１２】
0006469033.tif 000013
として推定できる。なお、
【数１３】
0006469033.tif 000014
の要素の平均値は、下記（５）式で表される。
【００５０】
【数１４】
0006469033.tif 000015
【００５１】
また、分散パラメータσｉｊを時間帯及び場所によって異なるように取得した場合、例えば、「テーマパークＤでは、午前中早くに来場する客は長時間楽しむことを目的としているが、夕方に来場する客は、より安価に購入できるナイトパスを利用して、短時間で少数のアトラクションを楽しむ傾向にある」という分析及びパターンの抽出が可能になる。
【００５２】
本実施形態では、さらに、特殊な場合として、時間帯場所ベクトル値対数行列Ｙのｊ列目の任意のベクトルｘｉｊに対して分散パラメータσｉｊが共通である場合についても考えることもできる。この場合には分散パラメータσｉｊはｊにしか依存しないので、分散パラメータσｊとして、上述した方法を適用することができる。この場合には、最小化の目的関数は、下記（６）式のように表される。
【００５３】
【数１５】
0006469033.tif 000016
【００５４】
分散パラメータσｉｊをｊ成分について共通化するということは、ｉ成分に時間帯の情報を持ち、ｊ成分に場所の情報を持つデータにおいて、分散パラメータσｉｊは時間帯に依存しないと考えていることとなる。データの密度によって偏りがあるデータ、つまり場所及び時間帯によって偏りがあるデータ等では、このようにパラメータを共通化することは非常に有効である。
【００５５】
これにより、分散パラメータσｉｊを、時間帯依存がなく、場所のみに特徴があると仮定して抽出した場合、例えば「サービスエリアＡでは、どの時間帯でも平均的に１時間滞在する傾向があるので、データ数は少ないが、同じような利用をされている場所が、道の駅、国道沿いのコンビニエンスストア等でも抽出できた。』という分析及びパターンの抽出が可能になる。
【００５６】
また、同様に、ｉ行の任意ベクトル値に対して分散パラメータσｉｊが共通である場合を考えた場合、分散パラメータσｉｊを分散パラメータσｉと書き直し、最小化する目的関数を下記（７）式のように設定することで、上述した方法を適応することができる。
【００５７】
【数１６】
0006469033.tif 000017
【００５８】
これにより、分散パラメータσｉｊを、場所依存がなく、時間帯のみに特徴があると仮定して抽出した場合、例えば「ある因子に属する場所では、昼の時間帯は３０分ほどなのに、夜の利用客は２時間以上の利用がされていて、これは多くのレストランでは、昼は客の回転数を上げる努力がなされている一方、夜は客の単価を上げるように努力している」という分析及びパターンの抽出が可能になる。
【００５９】
次に、本実施形態に係る分布推定装置１の機能について説明する。
【００６０】
図７に示すように、本実施形態に係る分布推定装置１は、時間帯場所情報処理部１０、特徴行列推定部２０、反復判定部２０ａ、特徴行列処理部３０、記憶部４０、及び、入出力部５０を有する。入出力部５０には、入力装置、表示装置等の外部装置２に接続されており、入出力部５０は、外部装置２に対して情報の入出力を行う。
【００６１】
記憶部４０は、時間帯場所情報テーブル４１、時間帯特徴テーブル４２、場所特徴テーブル４３、及び、場所分散テーブル４４を有する。
【００６２】
以下に各テーブルについて説明する。なお、テーブル形式のデータは行列形式にて表現できることから、以下では、各テーブルと各特徴行列を同一視し、区別せずに用いる。また、時間帯場所行列の各成分はベクトル値であり、各特徴行列の各成分は実数値であるが、同じ「行列」という表現を用い、文脈によって使い分ける。
【００６３】
＜時間帯場所情報テーブル４１＞
【００６４】
時間帯場所情報テーブル４１は、時間帯フィールド、場所ＩＤフィールド、及び、滞在時間ベクトルフィールドを有する。時間帯フィールドは、時間帯場所情報処理部１０により追加された時間帯を特定する識別子が設定されるフィールドである。場所ＩＤフィールドには、時間帯場所情報処理部１０により追加された場所を特定する識別子が設定されるフィールドである。滞在時間ベクトルフィールドは、時間帯場所情報処理部１０により各々の時間帯に各々の場所に滞在した複数人のユーザの滞在時間の値（以下、滞在時間ベクトルという。）が設定されるフィールドである。なお、各ユーザの滞在時間の値には、０又は正の実数を設定できるが、負の実数又は虚数を設定することはできない。
【００６５】
＜時間帯特徴テーブル４２＞
【００６６】
時間帯特徴テーブル４２は、時間帯フィールド、平均値基底ベクトルＩＤフィールド、及び、時間帯特徴値フィールドを有する。時間帯フィールドは、特徴行列推定部２０により時間帯を特定する識別子が設定されるフィールドである。平均値基底ベクトルＩＤフィールドは、特徴行列推定部２０により各データが従う分布の平均値を表現する基底を特定する識別子が設定されるフィールドである。時間帯特徴値フィールドは、特徴行列推定部２０により算出された当該時間帯の平均値を表現する基底ベクトルの特徴値が設定されるフィールドである。
【００６７】
＜場所特徴テーブル４３＞
【００６８】
場所特徴テーブル４３は、場所ＩＤフィールド、係数ＩＤフィールド、及び、場所特徴値フィールドを有する。場所ＩＤフィールドは、特徴行列推定部２０により場所を特定する識別子が設定されるフィールドである。係数ＩＤフィールドは、特徴行列推定部２０によりどの平均値基底ベクトルに対する係数かを特定する識別子が設定されるフィールドである。場所特徴値フィールドは、特徴行列推定部２０により算出された当該場所の平均値基底ベクトルに対する重みを特定する特徴値が設定されるフィールドである。
【００６９】
＜場所分散テーブル４４＞
【００７０】
場所分散テーブル４４は、時間帯フィールド、場所ＩＤフィールド、及び、場所分散値フィールドを有する。時間帯フィールドは、特徴行列推定部２０により時間帯を特定する識別子が設定されるフィールドである。場所ＩＤフィールドは、特徴行列推定部２０により場所を特定する識別子が設定されるフィールドである。場所分散値フィールドは、特徴行列推定部２０により算出された当該場所の場所分散パラメータの値が設定されるフィールドである。
【００７１】
時間帯場所情報処理部１０は、入出力部５０から、入力行列として、時間帯場所行列を入力する。また、時間帯場所情報処理部１０は、入力した時間帯場所行列に基づき、時間帯場所情報テーブル４１に、追加された場所、滞在時間に応じて、時間帯フィールド、場所ＩＤフィールド、滞在時間の値を設定した行を挿入する。
【００７２】
特徴行列推定部２０は、時間帯場所情報テーブル４１に格納された情報を取得し、後述する特徴行列推定処理を行うことにより特徴行列を推定し、得られた時間帯の特徴値を時間帯特徴テーブル４２に格納し、場所の特徴値を場所特徴テーブル４３に格納し、場所分散パラメータの値を場所分散テーブル４４に格納する。特徴行列推定部２０が特徴行列を推定する際には、時間帯（個体）ｉと場所（オブジェクト）ｊとの関連度を表すと共に、平均パラメータμｉｊ及び分散パラメータσｉｊを持つ対数正規分布に従う要素ｘｉｊｋからなるベクトルｘｉｊを要素として持つＩ×Ｊの時間帯場所ベクトル値行列Ｘの各ベクトルｘｉｊの各要素ｘｉｊｋの対数値、時間帯ｉがクラスタｒに所属することを表す非負値の要素ａｉｒを持つＩ×Ｒの時間帯特徴行列（第１の特徴行列）Ａ、場所ｊがクラスタｒに所属することを表す非負値の要素ｂｊｋを持つＪ×Ｒの場所特徴行列（第２の特徴行列）Ｂ、及び分散パラメータσｉｊを用いて表される上記（４）式の目的関数を最小化するように、時間帯特徴行列Ａ及び場所特徴行列Ｂを推定する。
【００７３】
反復判定部２０ａは、予め定めた反復終了条件を満足するまで、特徴行列推定部２０による推定を繰り返すように制御する。
【００７４】
特徴行列処理部３０は、時間帯特徴テーブル４２、及び場所特徴テーブル４３を参照し、リクエストの引数に対応する特徴行列を出力する。本実施形態では、例えば、外部装置２から特徴出力のリクエストが入力された場合に、特徴行列の出力する実行するが、特徴行列の出力を実行するタイミングはこれに限らず、予め定めた時間が経過する毎に、特徴行列を出力しても良い。また、出力する特徴は全ての特徴であっても一部の特徴であっても良く、全ての特徴を出力する場合には、時間帯特徴テーブル４２、及び場所特徴テーブル４３の全ての行を出力すればよい。
【００７５】
入出力部５０は、外部装置２から、入力行列として時間帯場所行列を入力し、入力した時間帯場所行列を時間帯場所情報処理部１０に引き渡す。また、入出力部５０は、時間帯特徴テーブル４２に格納された特徴値、及び場所特徴テーブル４３に格納された特徴値を外部装置２に対して出力する。
【００７６】
なお、本実施形態に係る分布推定装置１は、例えば、ＣＰＵ（Central Processing Unit）、ＲＡＭ（Random Access Memory）、各種プログラムを記憶するＲＯＭ（Read Only Memory）を備えたコンピュータ装置で構成される。また分布推定装置１を構成するコンピュータは、ハードディスクドライブ、不揮発性メモリ等の記憶部を備えていても良い。本実施形態では、ＣＰＵがＲＯＭ、ハードディスク等の記憶部に記憶されているプログラムを読み出して実行することにより、上記のハードウェア資源とプログラムとが協働し、上述した機能が実現される。
【００７７】
本実施形態に係る分布推定装置１は、時間帯場所行列を入力として特徴行列を推定し、特徴行列を出力する全体処理を行う。
【００７８】
まず、本実施形態に係る分布推定装置１による全体処理の流れを、図８に示すフローチャートを用いて説明する。
【００７９】
本実施形態では、時間帯場所情報処理部１０が時間帯場所情報を更新するタイミングは、例えば、システム管理者が外部装置２から供給されるデータをもとに手動で管理できるようにしてもよいし、あるユーザが新たな場所に滞在した場合に外部装置２が自動的に処理を起動するようにしてもよい。
【００８０】
ステップＳ１０１では、時間帯場所情報処理部１０が、入出力部５０から、入力行列として、時間帯場所行列を入力する。また、時間帯場所情報処理部１０が、入力した時間帯場所行列に基づき、時間帯場所情報テーブル４１に、追加された場所、滞在時間に応じて、時間帯フィールド、場所ＩＤフィールド、滞在時間の値を設定した行を挿入する。
【００８１】
ステップＳ１０３では、特徴行列推定部２０が、後述する特徴行列推定処理を行うことにより、各特徴行列、及び場所分散を推定する。
【００８２】
ステップＳ１０５では、特徴行列処理部３０が、時間帯特徴テーブル４２、及び場所特徴テーブル４３を参照し、リクエストの引数に対応する特徴行列を出力する。本実施形態では、例えば、外部装置２から特徴出力のリクエストが入力された場合に、特徴行列の出力する実行するが、特徴行列の出力を実行するタイミングはこれに限らず、予め定めた時間が経過する毎に、特徴行列を出力しても良い。また、出力する特徴は全ての特徴であっても一部の特徴であっても良く、全ての特徴を出力する場合には、時間帯特徴テーブル４２、及び場所特徴テーブル４３の全ての行を出力すればよい。
【００８３】
ここで、特徴行列推定部２０は、上述した特徴行列推定処理として、以下の方法で各特徴行列、場所分散を推定し、記憶部４０の時間帯特徴テーブル４２、場所特徴テーブル４３、及び場所分散テーブル４４に格納する処理を行う。
【００８４】
本実施形態に係る分布推定装置１による特徴行列推定処理の流れを、図９に示すフローチャートを用いて説明する。
【００８５】
なお、本実施形態では、時間帯場所情報テーブル４１中に存在する全データを下記（８）式のように表す。また、時間帯特徴行列Ａ、及び場所特徴行列Ｂをそれぞれ下記（９）式及び（１０）式のように表す。また、場所分散パラメータσを下記（１１）式のように表す。各式におけるＩは全時間帯数であり、Ｊは全場所数である。また、各式におけるｉは時間帯を特定する識別子であり、ｊは場所を特定する識別子であり、ｒは平均値を表現する基底ベクトルを特定する識別子である。
【００８６】
【数１７】
0006469033.tif 000018
【００８７】
【数１８】
0006469033.tif 000019
【００８８】
【数１９】
0006469033.tif 000020
【００８９】
【数２０】
0006469033.tif 000021
【００９０】
ステップＳ２０１では、特徴行列推定部２０が、時間帯特徴行列Ａ、及び場所特徴行列Ｂをそれぞれ初期化する。同様に終了条件の閾値ε、最大繰り返し回数を設定する。
【００９１】
ステップＳ２０３では、特徴行列推定部２０が、終了条件に用いる変数として特徴更新の最大変化幅を示す変数δをδ＝０として初期化する。
【００９２】
ステップＳ２０５では、特徴行列推定部２０が、下記（１２）式に従い、時間帯特徴行列Ａを更新する。その後、更新前の時間帯特徴行列Ａの要素の値と更新後の時間帯特徴行列Ａの要素の値の差の絶対値の最大値である
【数２１】
0006469033.tif 000022
がδより大きければ、
【数２２】
0006469033.tif 000023
と更新する。
【００９３】
【数２３】
0006469033.tif 000024
【００９４】
ただし、
【数２４】
0006469033.tif 000025
と書いた、μｉｊは、時間帯特徴行列Ａ、及び場所特徴行列Ｂによるｌｏｇ（ｘｉｊｋ）の推定値であると見なせる。
【００９５】
なお、上記（１２）式における記号「←」は、右辺の計算結果を左辺の変数に代入する処理を意味する。また、各式において、代入処理前のユーザ特徴行列Ａの要素の値をａｏｌｄｉｒとし、代入処理後の値をａｎｅｗｉｒとして記述している。また、上記（１２）式におけるμｉｊは、時間帯特徴行列Ａ及び場所特徴行列Ｂによるｌｏｇ（ｘｉｊｋ）の推定値であると見なせる。
【００９６】
ステップＳ２０７では、特徴行列推定部２０が、下記（１３）式に従い、場所特徴行列Ｂを更新する。その後、更新前の場所特徴行列Ｂの要素の値と更新後の場所特徴行列Ｂの要素の値の差の絶対値の最大値
【数２５】
0006469033.tif 000026
がδより大きければ、
【数２６】
0006469033.tif 000027
と更新する。なお、各式において、代入処理前の場所特徴行列Ｂの要素の値をｂｏｌｄｉｒとし、代入処理後の値をｂｎｅｗｉｒとして記述している。
【００９７】
【数２７】
0006469033.tif 000028
【００９８】
ステップＳ２０９では、特徴行列推定部２０が、下記（１４）式に従い、場所分散パラメータσｉｊを更新する。
【００９９】
【数２８】
0006469033.tif 000029
【０１００】
ステップＳ２１１では、特徴行列推定部２０が、計算繰り返し回数を更新する。
【０１０１】
ステップＳ２１３では、反復判定部２０ａが、計算繰り返し回数が予め定めた最大繰り返し数を超えたか否か、又は、特徴行列を更新することにより、上述した変数δが予め定めた閾値εより小さいか否かを判定する。
【０１０２】
ステップＳ２１３で計算繰り返し回数が予め定めた最大繰り返し数を超えていないと判定した場合、又は、上述した変数δが予め定めた閾値ε以上であると判定した場合（Ｓ２１３，Ｎ）は、ステップＳ２０５に戻る。また、ステップＳ２１３で計算繰り返し回数が予め定めた最大繰り返し数を超えたと判定した場合、又は、上述した変数δが予め定めた閾値εより小さいと判定した場合（Ｓ２１３，Ｙ）は、本特徴行列推定処理のプログラムの実行を終了する。
【０１０３】
このようにして、分布推定装置１により、滞在時間が従う対数正規分布が推定される。
【０１０４】
ここで、本実施形態の特殊な例として、データＸの場所分散パラメータσを、ｉに関して共通化したものを考える場合、すなわち目的関数が上記（６）式に従う場合には、特徴行列を更新するための式である上記（１２）乃至（１４）式の代わりに、下記（１５）乃至（１７）式を用いて、後述する手法で特徴行列を更新する。
【０１０５】
【数２９】
0006469033.tif 000030
【０１０６】
【数３０】
0006469033.tif 000031
【０１０７】
【数３１】
0006469033.tif 000032
【０１０８】
また、目的関数が上記（７）式に従う場合は、特徴行列を更新するための式である更新式である上記（１２）乃至（１４）式の代わりに、下記（１８）乃至（２０）式を用いて下記のように更新する。
【０１０９】
【数３２】
0006469033.tif 000033
【０１１０】
【数３３】
0006469033.tif 000034
【０１１１】
【数３４】
0006469033.tif 000035
【０１１２】
上記（１２）乃至（２０）式の各式において、全ての時間帯ｉ、場所ｊ、ユーザｋについてμｉｊ＝ｌｏｇ（ｘｉｊｋ）が成立する場合、各式の左辺と右辺が一致し、更新の最大変化幅を示す変数δの値が閾値ε以下となるため、更新が停止する。
【０１１３】
なお、上記の実施の形態では、時間帯場所行列を表現した行列から分布を推定する例を示しているが、この例に限定されることはない。例えば、地域と年齢とに応じてユーザの年収を表現する行列など、ユーザ、年齢、地域等のように１つ１つにＩＤ番号を付与して識別可能であり行列形式としてデータを表現することが可能な事物であるならば、あらゆるものが本装置による分布推定が可能である。また、本実施形態の図１１に示す分布推定装置の各構成要素の動作をプログラムとして構築し、分布推定装置として利用されるコンピュータにインストールして実行させる、または、ネットワークを介して流通させることが可能である。
【０１１４】
また、本実施形態では、図１に示す機能の構成要素の動作をプログラムとして構築し、分布推定装置１として利用されるコンピュータにインストールして実行させるが、これに限らず、ネットワークを介して流通させても良い。
【０１１５】
また、構築されたプログラムをハードディスクやフレキシブルディスク、ＣＤ−ＲＯＭ等の可搬記憶媒体に格納し、コンピュータにインストールしたり、配布したりしても良い。
【符号の説明】
【０１１６】
１  分布推定装置
２  外部装置
１０  時間帯場所情報処理部
２０  特徴行列推定部
２０ａ  反復判定部
３０  特徴行列処理部
４０  記憶部
４１  時間帯場所情報テーブル
４２  時間帯特徴テーブル
４３  場所特徴テーブル
４４  場所分散テーブル
５０  入出力部
【図１】
0006469033.tif 000038
【図２】
0006469033.tif 000039
【図３】
0006469033.tif 000040
【図４】
0006469033.tif 000041
【図５】
0006469033.tif 000042
【図６】
0006469033.tif 000043
【図７】
0006469033.tif 000044
【図８】
0006469033.tif 000045
【図９】
0006469033.tif 000046
【図１０】
0006469033.tif 000047
【図１１】
0006469033.tif 000048
【図１２】
0006469033.tif 000049
【図１３】
0006469033.tif 000050