(57)【要約】【課題】アイトラッカーによって計測された注視点の時系列から注視点の生成モデルのパラメータを推定する技術を提供する。【解決手段】時刻ｔの人の真の注目点は、時刻ｔ近傍にサッカードが発生する場合は時刻t-1の真の注目点を当該サッカードの方向及び大きさに応じて移動させた位置とし、時刻ｔ近傍にサッカードが発生しない場合は時刻t-1の真の注目点をランダムな方向及び大きさに応じて移動させた位置とし、時刻ｔの人の注視点の遅れは、AR(2)モデルに従うものとし、時刻ｔの人の真の注視点は、時刻ｔの真の注目点を時刻ｔの注視点の遅れにより補正した位置とし、実際に計測される時刻ｔの注視点を、上記時刻ｔの真の注視点にノイズが加わったものとしてモデル化したものを注視点の生成モデルとして、対象者の眼の動きを計測して得た注視点の時系列から、注視点の生成モデルのモデルパラメータを推定する。【選択図】図５
2019126425.tif 000002
【特許請求の範囲】
【請求項１】
サッカードの時系列は、各サッカードの特徴量をマークとするマーク付き点過程で生成されるものとし、
時刻ｔの人の真の注目点は、時刻ｔ近傍にサッカードが発生する場合は時刻t-1の真の注目点を当該サッカードの方向及び大きさに応じて移動させた位置とし、時刻ｔ近傍にサッカードが発生しない場合は時刻t-1の真の注目点をランダムな方向及び大きさに応じて移動させた位置とし、
時刻ｔの人の注視点の遅れは、AR(2)モデルに従うものとし、
時刻ｔの人の真の注視点は、前記時刻ｔの真の注目点を前記時刻ｔの注視点の遅れにより補正した位置とし、
実際に計測される時刻ｔの注視点を、上記時刻ｔの真の注視点にノイズが加わったものとしてモデル化したものを注視点の生成モデルとして、
対象者の眼の動きを計測して得た注視点の時系列から、前記注視点の生成モデルのモデルパラメータを推定するモデル推定部を含む、
推定装置。
【請求項２】
請求項１記載の推定装置であって、
前記モデル推定部で学習したモデルパラメータに基づいて、サッカードの発生時刻を推定するサッカード発生時刻推定部と、
前記モデル推定部で学習したモデルパラメータに基づいて、注視点の時系列の推定結果を求める注視点系列推定部と、
前記モデル推定部で学習したモデルパラメータに基づいて、注目点の時系列の推定結果を求める注目点系列推定部と、
前記モデル推定部で学習したモデルパラメータに基づいて、注目範囲または注目対象の推定結果を求める注目範囲推定部と、
前記モデル推定部で学習したモデルパラメータに基づいて、サリエンシーマップを生成するサリエンシーマップ生成部と、の少なくともいずれかをさらに含む、
推定装置。
【請求項３】
請求項１記載の推定装置であって、
κrを、時刻rにサッカードが発生した場合はそのサッカードのベクトル、発生しなかった場合は空集合をとる確率変数とし、jrを時刻rにおけるサッカードがサリエンシーマップのどの混合成分に向かって発生したものかを示す潜在変数とし、Hrを時刻rまでの履歴とし、Sをサリエンシーマップ画像とし、r*を時刻r以前で発生した最後のサッカードの発生時刻とし、時刻r以前に発生したサッカードがない場合はr*=0とし、arを時刻rにおける真の注目点とし、srを時刻rにおける注視点の遅れとし、hNBをサッカードが短い時間間隔で発生しないようにするための修正項とし、i*を時刻r以前で発生した最後のサッカードがどの混合成分に向かって発生したものかを示すものとし、
前記生成モデルは、
【数６２】
2019126425.tif 000064
により与えられる、
ことを特徴とする推定装置。
【請求項４】
請求項３記載の推定装置であって、
j(r+1)^*を時刻r+1以前に発生した最後のサッカードが、どの混合成分に向かって発生したものかを示すものとし、
前記生成モデルは、
【数６３】
2019126425.tif 000065
により与えられる、
ことを特徴とする推定装置。
【請求項５】
請求項４記載の推定装置であって、
【数６４】
2019126425.tif 000066
とし、q(b,z,φ)を任意の事後分布とし、p(b,z,φ|x,ψ)を真の事後分布とし、前記モデル推定部は、次式により与えられるL(q)が大きくなるように、事後分布q(b),q(z),q(φ),ψを繰り返し更新することにより前記モデルパラメータを学習する、
【数６５】
2019126425.tif 000067
ことを特徴とする推定装置。
【請求項６】
サッカードの時系列は、各サッカードの特徴量をマークとするマーク付き点過程で生成されるものとし、
時刻ｔの人の真の注目点は、時刻ｔ近傍にサッカードが発生する場合は時刻t-1の真の注目点を当該サッカードの方向及び大きさに応じて移動させた位置とし、時刻ｔ近傍にサッカードが発生しない場合は時刻t-1の真の注目点をランダムな方向及び大きさに応じて移動させた位置とし、
時刻ｔの人の注視点の遅れは、AR(2)モデルに従うものとし、
時刻ｔの人の真の注視点は、前記時刻ｔの真の注目点を前記時刻ｔの注視点の遅れにより補正した位置とし、
実際に計測される時刻ｔの注視点を、上記時刻ｔの真の注視点にノイズが加わったものとしてモデル化したものを注視点の生成モデルとして、
対象者の眼の動きを計測して得た注視点の時系列から、前記注視点の生成モデルのモデルパラメータを推定するモデル推定ステップを含む、
推定方法。
【請求項７】
請求項１から請求項５の何れかの推定装置としてコンピュータを機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、ヒトの注視点の時系列をモデル化した注視点モデルのパラメータを推定する技術に関する。
【背景技術】
【０００２】
サリエンシーマップは、視覚刺激に対するボトムアップ性の注意を表現するための計算論的なモデルであり、元画像の色や輝度といった特徴から2 次元のグレースケール画像として計算される。サリエンシーマップによって画像の中で人が視線を向けやすい箇所を知ることができるため、計算論的神経科学に限らない様々な分野で研究されている。
【０００３】
非特許文献１には、サリエンシーマップという概念について詳しく記載されている。非特許文献１では、サリエンシーマップの中で最も値が大きい箇所に視線が向くと仮定されている。この仮定はWinner take allルールと呼ばれており、このルールに従うと我々の視線は画像情報のみから決定論的に決まることとなる。
【先行技術文献】
【非特許文献】
【０００４】
【非特許文献１】C. Koch and S. Ullman, "Shifts in selective visual attention: towards the underlying neural circuitry",  Human Neurobiology, vol. 4, pp. 219-27, 1985.
【発明の概要】
【発明が解決しようとする課題】
【０００５】
しかしながら、非特許文献１で用いられているWinner take allルールは、我々が画像を見たときにいつも同じ箇所に着目するとは限らないという点で、現実に即していない。さらに、Winner take all ルールではサリエンシーマップの最大値のみに意味を見出すため、各画素に割り当てられた絶対値や異なる画素間の値の比に解釈を与えることができない。
【０００６】
ここで挙げた問題は、従来のサリエンシーマップが「画像の情報のみ」から計算されるあくまで計算論的なモデルであることに起因する。これらの問題を解決するためには、画像を見ているときの注意状態を代替するような教師データを用意し、教師付き学習を行う必要がある。
【０００７】
そこで本発明では、サリエンシーマップがサッケード生成に関連の深い上丘に表象されているという知見をもとに、サリエンシーマップを点過程モデルを用いてモデル化し、アイトラッカーによって計測された注視点の時系列からサリエンシーマップのモデルパラメータを簡便に推定する推定装置、その方法、及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【０００８】
上記の課題を解決するために、本発明の一態様によれば、推定装置は、サッカードの時系列は、各サッカードの特徴量をマークとするマーク付き点過程で生成されるものとし、時刻ｔの人の真の注目点は、時刻ｔ近傍にサッカードが発生する場合は時刻t-1の真の注目点を当該サッカードの方向及び大きさに応じて移動させた位置とし、時刻ｔ近傍にサッカードが発生しない場合は時刻t-1の真の注目点をランダムな方向及び大きさに応じて移動させた位置とし、時刻ｔの人の注視点の遅れは、AR(2)モデルに従うものとし、時刻ｔの人の真の注視点は、時刻ｔの真の注目点を時刻ｔの注視点の遅れにより補正した位置とし、実際に計測される時刻ｔの注視点を、上記時刻ｔの真の注視点にノイズが加わったものとしてモデル化したものを注視点の生成モデルとして、対象者の眼の動きを計測して得た注視点の時系列から、注視点の生成モデルのモデルパラメータを推定するモデル推定部を含む。
【発明の効果】
【０００９】
本発明によれば、アイトラッカーによって計測された注視点の時系列から注視点モデルのモデルパラメータを簡便に推定することができるという効果を奏する。
【図面の簡単な説明】
【００１０】
【図１】マーク付き点過程から生成された標本の例を示す図。
【図２】サッカードがmovement fieldに向かって発生したタイミングでそのニューロンが発火するイメージを示す図。
【図３】ニューロンの配置と受容野の中心位置との対応関係を示す図。
【図４】真の注視点を、真の注目点arと注視点の遅れsrとに分解のイメージを示す図。
【図５】第一実施形態に係る推定装置の機能ブロック図。
【図６】第一実施形態に係る推定装置の処理フローの例を示す図。
【図７】事後期待値とその解釈の例を示す図。
【図８】被験者のデータの１反復目の推定結果を示す図。
【図９】被験者のデータの１１反復目の推定結果を示す図。
【図１０】被験者のデータから計算した注視点の頻度分布と、サリエンシーマップにサポートベクター回帰を適用して得られた初期値とを表す図。
【図１１】被験者のデータの１１反復目の推定結果を示す図。
【図１２】第一実施形態の変形例に係る推定装置の機能ブロック図。
【発明を実施するための形態】
【００１１】
以下、本発明の実施形態について、説明する。なお、以下の説明に用いる図面では、同じ機能を持つ構成部や同じ処理を行うステップには同一の符号を記し、重複説明を省略する。以下の説明において、テキスト中で使用する記号「^」「~」「-」等は、本来直後の文字の真上に記載されるべきものであるが、テキスト記法の制限により、当該文字の直前に記載する。式中においてはこれらの記号は本来の位置に記述している。また、ベクトルや行列の各要素単位で行われる処理は、特に断りが無い限り、そのベクトルやその行列の全ての要素に対して適用されるものとする。
【００１２】
＜サリエンシーマップのモデル＞
まず、サリエンシーマップのモデルについて説明する。
【００１３】
人がある一枚の画像を見ており、その間の注視点の時系列がアイトラッカーによって計測されている状況を想定する。本実施形態では、アイトラッカーによって計測された眼球が向いている位置を「注視点」、人が真に注目している位置を「注目点」と呼び分けることにする。本実施形態の目的は、画像情報を用いることで「注視点」の軌跡から眼球運動による遅れの成分やノイズ成分を除去し、「注目点」の軌跡を推定することである。この推定の際に、「注目点」の軌跡に影響を与える潜在変数としてサリエンシーマップが同時に計算される。
【００１４】
本実施形態では、注視点の時系列を生成モデルで表現し、その事後分布を推定する。
【００１５】
＜サリエンシーマップの条件付き強度関数による定義＞
従来のサリエンシーマップは、サッカードの行き先になりやすい地点で値が大きくなるように定義されている(参考文献１参照)。
（参考文献１）L. Itti, C. Koch and E. Niebur, "A Model of Saliency-Based Visual Attention for Rapid Scene Analysis", IEEE Transactions on Pattern Analysis and Machie Intelligence, vol. 20, pp. 1254-1259, 1998.
【００１６】
そこで本実施形態では、このサリエンシーマップの考え方に則り、各地点のサリエンシーの値をその地点に向かうサッカードの発生頻度によって定義することにする。
【００１７】
サッカードは、眼球運動に含まれる跳躍性のジャンプのことであり、瞬間的に発生する運動である。そのため、その時系列は瞬間的に発生するイベントを並べたような時系列となる。統計的モデリングの枠組みにおいては、このようなランダムに生起するイベントの時系列を記述するための確率過程として点過程がよく用いられている。とくに各イベントが何らかの特徴量をもつ場合はマーク付き点過程が用いられる。マーク付き点過程とは、各イベントにマークとよばれる特徴量が付加された系列の確率的構造を記述するための確率過程であり、モデリングの対象としては地震や為替取引等が挙げられる。図１にマーク付き点過程から生成された標本の例を示す。図１では、横軸が時刻を表し、縦軸が所定のマークを表す。サッカードの時系列をマーク付き点過程とみなす場合は、サッカードの発生時刻をイベントの発生時刻、サッカードの振幅・方向・到着地点などの特徴量をマークとみなすことができる。
【００１８】
点過程モデリングを行う場合、データから点過程の確率的構造を推定することが目標となる。しかしながら、その確率密度関数を直接推定するのは難しい。これは、点過程が定義される距離空間の構造が複雑であり、尤度の計算が困難であることに起因する。したがって、点過程によって時系列データのモデリングを行う場合は、確率密度関数を推定するかわりに条件付き強度関数と呼ばれる関数を推定することが多い。
【００１９】
条件付き強度関数は、過去の履歴があたえられたもとでの単位時間あたりに発生するイベントの数の期待値に相当する関数である。具体的には
【００２０】
【数１】
2019126425.tif 000003
【００２１】
と定義される。ただし、Nは計数測度、κ∈Κは各イベントに付加されたマーク、Htは時刻tまでの履歴である。このHtは、t以前に発生したイベントの情報や同時に観測された共変量の値など、全ての情報を含んでもよい。
【００２２】
この条件付き強度関数によって、ほとんどの点過程の確率的構造は一意に決定されるため（例えば参考文献２参照）、点過程の推定を行う際は条件付き強度関数を推定することを考えれば十分である。
（参考文献２）D. J. Daley and D. Vere-Jones, "An Introduction to The Theory of Point Processes", Springer Science and Business Media, New York, 2003.
【００２３】
本実施形態では、サッカードのベクトルをマークとするようなマーク付き点過程からサッカードの系列が生成されていると仮定する。今回の問題設定では画像という時間的に定常な刺激を呈示しているため、サッカードの発生も時間的に定常であるとする。すなわち、条件付き強度関数が
【００２４】
【数２】
2019126425.tif 000004
【００２５】
と分解できるとし、サッカードの発生を示す関数λ(t|Ht)が定常であるとする。ただし、条件付き強度関数それ自体が確率過程である。このとき、λS(κ|Ht)はサッカードが画像平面上のどの地点に向かって発生しやすいかを説明する項であり、サリエンシーマップに対応する関数となる。
【００２６】
次に、λS(κ|Ht)としてどのような関数が適しているかを考察する。本実施形態では、サッカード生成に関連の深い中脳の上丘を参考にλS(κ|Ht)の形を決めることにする。上丘にはある特定の振幅・方向のサッカードが発生したときのみに発火するという性質をもつニューロンが存在する。すなわち、各ニューロンは視野上でmovement field(図２中、MFと表す)と呼ばれるある特定の領域をもち、サッカードがその領域に向かって発生したタイミングでそのニューロンが発火する。図２にそのイメージ図を示す。図中、spike trainはニューロンの発火状態を示し、scan pathは注視点を示す。また、それらのニューロンの配置(motor map on SC)は各々が反応する振幅・方向に空間的に対応したマップをなしていることが知られている。図３にその対応関係を示す。図３では上丘におけるニューロンの配置と受容野の中心位置は空間的に対応している。このような対応関係を根拠に、サリエンシーマップが上丘に表象されているといわれることもある。
【００２７】
サリエンシーマップが上丘によって表象されている場合、サリエンシーマップは上丘ニューロンの受容野(ニューロンの神経応答に変化を生じるような刺激が提示される空間の領域)を足し合わせた形で表現されていると考えられる。本実施形態では、この生理学的な知見を参考に、λS(κ|Ht)を次のようなガウシアンカーネルの非負結合によって表すことにする。
【００２８】
【数３】
2019126425.tif 000005
【００２９】
ただし、
【００３０】
【数４】
2019126425.tif 000006
【００３１】
であり、cjは発火頻度、μjは受容野の中心位置、Λjは受容野の広がり具合を表すパラメータである。とくに、cjについては過去のサッカードの履歴Htに依存してよいとする。また、xt-は時刻t直前での注視点である。時刻tでサッカードが発生した場合は、xt-はサッカードの発生地点、κ+xt-はサッカードの到着地点に対応する。
【００３２】
以上が本実施形態におけるサリエンシーマップの定義である。ここで、サリエンシーマップを条件付き強度関数によって定義することの利点を補足しておく。条件付き強度関数によってサリエンシーマップを定義する最大の利点は、サリエンシーの値がサッカードの頻度になっているという点である。したがって、Winner take allルールで説明できないような注視点の確率的な挙動を説明することができる。さらに、異なる画素間の値の比や、異なる画像における値の比に意味を与えることができる。例えば、A地点のサリエンシー値がB 地点のサリエンシー値の2倍である場合、A地点の方がB地点よりも2倍注視されやすいと解釈できる。さらに、サリエンシーマップを上丘ニューロンの発火活動の強度関数と直接対応付けることも可能である。
【００３３】
画像の情報については、次のように利用する。事前に、従来手法で計算されたサリエンシーマップ画像Sに対してガウシアンカーネルを用いたサポートベクター回帰により、{(cj,μj,Λj)}Jj=1を推定しておく。ただし、条件付き強度関数は非負である必要があるため、サポートベクター回帰を適用する際は重み係数-cjについて非負制約を導入する。事後分布の推定の段階では、事前に推定しておいたパラメータのうち{(μj,Λj)}Jj=1を固定したハイパーパラメータとしてそのまま利用する。また、{cj}Jj=1は、cj(Ht)を推定する際の事前分布として利用する。
【００３４】
＜注視点の生成モデル＞
観測時間をT、アイトラッカーのサンプリング回数をRとし、サンプリング間隔をΔt = T/Rとする。アイトラッカーによって計測された注視点をxr∈R2,r=1,…,Rとし、人が真に着目している点の位置(つまり、注目点)をar∈R2,r=1,…,Rとする。この注目点arは、普段は連続的に変化しつつ、注意の対象が遷移した際に大きくジャンプするとする。本実施形態では、この瞬間的な注目点のジャンプをサッカードと呼ぶことにする。
【００３５】
ここで、κrを、時刻rにサッカードが発生した場合はそのサッカードのベクトル、発生しなかった場合は空集合{}(ただし、式中では、空集合を
【００３６】
【数５】
2019126425.tif 000007
【００３７】
とも表記する）をとる確率変数とする。このκrを用いて、注目点arが次のようなジャンプを伴うランダムウォークに従って生成されているとする。
【００３８】
【数６】
2019126425.tif 000008
【００３９】
このような仮定のもとで、注目点arは不連続な時系列となる。つまり、時刻r毎の人の真の注目点arは、１時刻前の真の注目点ar-1を基準として、時刻rの近傍(時刻r-1〜時刻rの間)にサッカードが発生する場合はサッカードの方向及び大きさ(つまり、サッカードのベクトルκr)に応じて移動させた位置とし、時刻rの近傍にサッカードが発生しない場合はランダムな方向及び大きさへ移動させた位置となるようモデル化される。言い換えれば、時刻rの近傍にサッカードが発生しない場合はランダムウォークur〜N(0,U)により移動先の位置が生成されるものとしてモデル化される。
【００４０】
しかしながら、眼球の物理的な制約上、注視点は瞬間的に移動することはない。なぜならば、我々がある点を注目しようと思い立ったのちに注視点がその点に移動するまでには、筋肉が眼球を動かす分の遅れが生じるからである。そこで、この眼球運動の制約によって生じる注目点arからの注視点の遅れをsrで表すことにし、時刻毎の人の注視点の遅れsrは、次のようなAR(2)モデルによって生成されると仮定され、モデル化される。
【００４１】
【数７】
2019126425.tif 000009
【００４２】
ただし、F1、F2は眼球運動の特性を決めるパラメータである。
【００４３】
このもとで、真の注視点の軌跡はar+srによって表される。つまり、真の注視点は、真の注目点arを注視点の遅れsrにより補正した位置としてモデル化される。図４に分解のイメージを示す。そして、実際に計測される注視点は、真の注視点ar+srにノイズwrが加わたものとして次のようにモデル化される。
【００４４】
【数８】
2019126425.tif 000010
【００４５】
次に、サッカード{κr}がどのように生成されるかを考える。本実施形態では、＜サリエンシーマップの条件付き強度関数による定義＞で説明したとおり、サッカードが(2)式および(3)式の条件付き強度関数によって定義されるマーク付き点過程から生成されていると仮定する。ただし、観測が離散的であるため、条件付き強度関数は離散時間で定義されることに注意する。いま、サッカードが発生した直後は、新たなサッカードは発生しにくく、かつそのサッカードの行き先は発生直前の注視点に依存すると考えられる。そこで、サッカードがサリエンシーマップのどの混合成分に向かって発生したものかを示す潜在変数
【００４６】
【数９】
2019126425.tif 000011
【００４７】
を導入し、{(κr,jr)}が次のようなマーク付き点過程から生成されていると想定する。
【００４８】
【数１０】
2019126425.tif 000012
【００４９】
ただし、r*は時刻r以前で発生した最後(直近)のサッカードの発生時刻のインデックスであり、i*はそのサッカードがどの混合成分に向かって発生したものかを示すインデックスである。r以前で発生したサッカードがない場合はr*=0とする。A={Aij}は確率推移行列であり、サッカードが発生したもとでそのサッカードが混合成分iから混合成分jに向かうものである確率がAijに対応する。また、πjは最初に発生したサッカードが向かう混合成分がjである確率である。h(・)はサッカードが短い時間間隔で発生しないようにするための修正項であり、本実施形態では負の二項分布のハザード関数
【００５０】
【数１１】
2019126425.tif 000013
【００５１】
を用いる。ただし、
【００５２】
【数１２】
2019126425.tif 000014
【００５３】
であり、θは発生間隔の平均を、mはサッカード直後の発生しにくさを操作するパラメータである。とくに、m=1のとき負の二項分布は幾何分布と一致しハザードは定数となる。
【００５４】
さらに、A,πおよびθについては次のような事前分布を設定する。
【００５５】
【数１３】
2019126425.tif 000015
【００５６】
ただし、
【００５７】
【数１４】
2019126425.tif 000016
【００５８】
であり、{(cj,μj,Λj)}は事前にサポートベクター回帰によって計算しておいたパラメータである。また、α0は画像情報の事前分布への反映度合いを調節するパラメータであり、事前に適切な値を設定しておく。
【００５９】
以上が本実施形態における生成モデルの定義である。最後に、モデルをまとめておく。
【００６０】
【数１５】
2019126425.tif 000017
【００６１】
＜事後分布の推定＞
推定のステップでは、(10)式のモデルのもとでの潜在変数の事後分布およびハイパーパラメータを推定することが目標となる。その推定の方法を順を追って説明する。まずはじめに、(10)式の生成モデルを扱いが容易なモデルに書き換える。次に、書き換えたモデルのもとで、事後分布の近似を求める方法を説明する。最後に得られた事後分布をどう解釈すればよいかを述べる。
【００６２】
＜スイッチング線形ガウス状態空間モデルへの帰着＞
上述の(10)式の生成モデルは、等価な隠れセミマルコフ・スイッチング線形ガウス状態空間モデルに変形することができる。このようなモデルに書き換えることで、変分ベイズの枠組みで潜在変数の事後分布を推定することが可能となる。
【００６３】
サッカードがどの混合成分に向けて発生したかを示す潜在変数{jr}Rr=1のかわりとして、新たな潜在変数{zr}Rr=1を
【００６４】
【数１６】
2019126425.tif 000018
【００６５】
と定義する。ただし、j(r+1)^*は時刻r+1以前に発生した最後のサッカードが、どの混合成分に向かって発生したものかを示すインデックスである。なお、上付き添え字及び下付き添え字におけるA^BはABを意味し、A_BはABを意味するものとする。また、{zr}Rr=1に関連する事象について
【００６６】
【数１７】
2019126425.tif 000019
【００６７】
と表記することにする。以上の定義のもとで、{zr}Rr=1は推移確率が
【００６８】
【数１８】
2019126425.tif 000020
【００６９】
であり、初期確率が
【００７０】
【数１９】
2019126425.tif 000021
【００７１】
であるセミマルコフ過程に従う。
【００７２】
次に、{(xr,ar,sr,κr)}Rr=1の生成モデルを書き換える。書き換えの際のポイントはκrを消すことである。
【００７３】
ここで、新たな変数を以下のように定める。
【００７４】
【数２０】
2019126425.tif 000022
【００７５】
このように変数を定義すると、{(xr,br)}Rr=1の生成モデルを次のような線形状態空間モデルに書き換えることができる。
【００７６】
【数２１】
2019126425.tif 000023
【００７７】
したがって、(10)式の生成モデルは、背後にセミマルコフ過程に従う潜在変数が存在し、その潜在変数に従って局所的に線形ガウス状態空間モデルをつなげた形で表すことができる。このようなモデルをスイッチング線形状態空間モデルとよぶ。
【００７８】
最後に、(10)式を書き換えたモデルをまとめておく。
【００７９】
【数２２】
2019126425.tif 000024
【００８０】
＜変分ベイズ＞
以降、表記の簡単のため
【００８１】
【数２３】
2019126425.tif 000025
【００８２】
とする。
【００８３】
推定のステップでは、xが与えられたもとでのb,z,φの事後分布を計算することが目標となる。しかしながら、モデルの構造が複雑であるため、真の事後分布を解析的に計算することができない。そこで、本稿では変分ベイズの枠組みを用いて、真の事後分布の近似を計算することにする。
【００８４】
変分ベイズとは、計算が簡単になるようにあらかじめ指定しておいた分布族の中から、真の事後分布とのKLダイバージェンスが最小になるような分布を求め、その分布を真の事後分布の近似とする方法である。いま、xの周辺対数尤度log p(x|ψ)は、任意の分布q(b,z,φ)を用いて
【００８５】
【数２４】
2019126425.tif 000026
【００８６】
と分解できる。ただし、KL(q||p(b,z,φ|x,ψ))はq(b,z,φ)と真の事後分布p(b,z,φ|x,ψ)との間のKLダイバージェンスである。(19)式の一行目の左辺がqに依存しないことより、L(q)の最大化とKLダイバージェンスの最小化は一致する。とくに、qが真の事後分布p(b,z,φ|x,ψ)と一致する場合にL(q)は最大値をとる。
【００８７】
したがって、事後分布の近似qを求める際は、真の事後分布とのKLダイバージェンスの最小化を考えるかわりに、L(q)の最大化を考えれば良い。しかしながら、任意の分布qについてL(q)を計算するのは困難である。そこで、L(q)の計算が簡単な分布の族Qを指定し、この分布族の中でL(q)が最大となる分布を求めることにする。すなわち、事後分布の推定を次のような最適化問題を解くことに帰着させる。
【００８８】
【数２５】
2019126425.tif 000027
【００８９】
本実施形態ではこの分布族Qとして、独立性制約
【００９０】
【数２６】
2019126425.tif 000028
【００９１】
を満たす分布の族を指定する。このような制約を導入することにより、座標降下法によってL(q)の最適化を効率的に行うことができる。具体的には、q(b)の更新にカルマンスムーザ、q(z)の更新にforward-backwardアルゴリズムが利用できる。さらに、q(φ)の更新も解析的な最適解を導出できる。また、ψについてL(q)を最大化することにより、これらのハイパーパラメータをもデータから決定することが可能である。
【００９２】
以下、上述の処理を実現する推定装置について説明する。
【００９３】
＜第一実施形態に係る推定装置＞
注視点の時系列を入力として、注視点モデルのパラメータを推定する推定装置について説明する。
【００９４】
推定装置は、例えば、中央演算処理装置（CPU: Central Processing Unit）、主記憶装置（RAM: Random Access Memory）などを有する公知又は専用のコンピュータに特別なプログラムが読み込まれて構成された特別な装置である。推定装置は、例えば、中央演算処理装置の制御のもとで各処理を実行する。推定装置に入力されたデータや各処理で得られたデータは、例えば、主記憶装置に格納され、主記憶装置に格納されたデータは必要に応じて中央演算処理装置へ読み出されて他の処理に利用される。推定装置の各処理部は、少なくとも一部が集積回路等のハードウェアによって構成されていてもよい。推定装置が備える各記憶部は、例えば、RAM（Random Access Memory）などの主記憶装置、ハードディスクや光ディスクもしくはフラッシュメモリ（Flash Memory）のような半導体メモリ素子により構成される補助記憶装置、またはリレーショナルデータベースやキーバリューストアなどのミドルウェアにより構成することができる。
【００９５】
図５は第一実施形態に係る推定装置の機能ブロック図を、図６はその処理フローの例を示す図である。
【００９６】
推定装置は、モデル推定部１０１と出力部１７０とを含む。
【００９７】
モデル推定部１０１は、初期設定部１１０と、第１更新部１２０と、第２更新部１３０と、第３更新部１４０と、第４更新部１５０と、制御部１６０とを含む。
【００９８】
<<モデル推定部１０１>>
[入力]：アイトラッカーにより計測された注視点の時系列xr（r=1,2,…,R)
[出力]：モデルパラメータの推定結果q(b),q(z),q(φ),ψ
[処理]：モデル推定部１０１は、アイトラッカーにより計測された注視点(対象者の眼の動きを計測して得た注視点)の時系列xr（r=1,2,…,R)を入力として、(17)式により表される注視点の時系列モデルの各パラメータb,z,φ,ψの事後分布を推定することにより、注視点の時系列モデル（学習済みモデル）を求める。
【００９９】
例えば、独立性制約
【０１００】
【数２７】
2019126425.tif 000029
【０１０１】
を満たす分布族の下で、
【０１０２】
【数２８】
2019126425.tif 000030
【０１０３】
のL(q)が最大となる事後分布を求める。
【０１０４】
以下、モデル推定部１０１の具体的な処理について説明する。
【０１０５】
〔初期設定部１１０〕
初期設定部１１０は、q(b,z,φ),ψに適当な初期値を設定し(Ｓ１１０)、出力する。このとき、q(b,z,φ)については、(21)式
【０１０６】
【数２９】
2019126425.tif 000031
【０１０７】
を満たすようにq(b),q(z),q(φ)を与える。ここで、q(b),q(z),q(φ)は、以後の処理で更新の際に用いる各変数の事後期待値が計算できるような分布を与える。例えば、各分布を更新する際に更新先となる分布と、同じ形の分布となるようにパラメータを決定する。本実施形態では、
q(b):正規分布
q(φ):q(A)とq(π)はディリクレ分布、q(θ)はベータ分布
となるようにパラメータを決定する。なお、q(z)は、離散確率変数なので、どのような値でもよい。
【０１０８】
〔第１更新部１２０〕
第１更新部１２０は、q(z),q(φ),ψを入力とし、q(z),q(φ),ψが与えられているもとで、(19)式のL(q)を最大にするq(b)を求め、q(b)の値を求めた値で更新する(Ｓ１２０)。例えば、初回の処理では初期設定部１１０で設定したq(z),q(φ),ψを用い、2回目以降の処理ではそれぞれ第２更新部１３０、第３更新部１４０、第４更新部１５０で更新した最新のq(z),q(φ),ψを用いる。以下に例を示す。
【０１０９】
q(z),q(φ)およびψが与えられているもとで、L(q)を最大にするq(b)は
【０１１０】
【数３０】
2019126425.tif 000032
【０１１１】
となる。ここで、b以外の変数についての完全対数尤度の事後期待値は、次のような線形ガウス状態空間モデルの対数尤度と一致する。
【０１１２】
【数３１】
2019126425.tif 000033
【０１１３】
ただし、~Fr,~or,~Prは、
【０１１４】
【数３２】
2019126425.tif 000034
【０１１５】
であり、{ξrj}r=1,…,R,j=1,…,Jは、
【０１１６】
【数３３】
2019126425.tif 000035
【０１１７】
である。したがって、通常の線形ガウス状態空間モデルと同様、カルマンスムーザを適用することによってq(b)が計算可能である。既存のカルマンスムーザを適用することができるため、その詳細は省略する。
【０１１８】
〔第２更新部１３０〕
第２更新部１３０は、q(b),q(φ),ψを入力とし、q(b),q(φ),ψが与えられているもとで,(19)式のL(q)を最大にするq(z)を求め、q(z)の値を求めた値で更新する(Ｓ１３０)。例えば、初回の処理では初期設定部１１０で設定したq(φ),ψと第１更新部１２０で更新したq(b)を用い、2回目以降の処理ではそれぞれ第１更新部１２０、第３更新部１４０、第４更新部１５０で更新した最新のq(b),q(φ),ψを用いる。以下に例を示す。
【０１１９】
q(b),q(φ)およびψが与えられているもとで,L(q)を最大にするq(z)は
【０１２０】
【数３４】
2019126425.tif 000036
【０１２１】
となる。ここで、z以外の変数についての完全対数尤度の事後期待値は、次のようなセミマルコフ過程の対数尤度と一致する。
【０１２２】
【数３５】
2019126425.tif 000037
【０１２３】
ただし、
【０１２４】
【数３６】
2019126425.tif 000038
【０１２５】
であり、これらの事後期待値は解析的に計算可能である。
【０１２６】
したがって、隠れセミマルコフモデルの推定に用いられるforward-backwardアルゴリズムを適用することによってq(z)が計算可能である。forward-backwardアルゴリズムでは、次式で定義されるforward-message(α,α*)およびbackward-message(β,β*)を再起的に計算する。
【０１２７】
【数３７】
2019126425.tif 000039
【０１２８】
ただし、orは時刻rでの観測に対応する変数である。計算されたメッセージを用いることで、zrの事後分布を得ることができる。具体的な更新式の詳細については、参考文献３等を参照されたい。
(参考文献３)S. Z. Yu, "Hidden semi-Markov models. Articial intelligence", vol. 174, pp. 215-243, 2010.
【０１２９】
しかしながら、通常のforward-backwardアルゴリズムに必要な計算量はO(J2R2)であり、時系列の長さRが長いときは望ましくない。そこで、本実施形態では参考文献４で提案された手法を用いて計算量を削減する。
(参考文献４)M. J. Johnson and A. S. Willsky, "Stochastic Variational Inference for Bayesian Time Series Models", International Conference on Machine Learning, 2014.
【０１３０】
この手法の主なアイディアは、(27)式の隠れセミマルコフモデルを等価な隠れマルコフモデルへ変換するというものである。
【０１３１】
ここで、新たな潜在変数-zrを
【０１３２】
【数３８】
2019126425.tif 000040
【０１３３】
とし、
【０１３４】
【数３９】
2019126425.tif 000041
【０１３５】
とする。ただし、
【０１３６】
【数４０】
2019126425.tif 000042
【０１３７】
とする。以上の設定のもとで、対数尤度が
【０１３８】
【数４１】
2019126425.tif 000043
【０１３９】
と表されるマルコフ過程を考えると、
【０１４０】
【数４２】
2019126425.tif 000044
【０１４１】
となる（参考文献５参照）。
(参考文献５)M. J. Johnson, "Bayesian time series models and scalable inference", PhD thesis, Massachusetts Institute of Technology, 2014.
【０１４２】
すなわち、{-zr}Rr=1の分布が得られれば、{zr}Rr=1の分布も得られるということになる。とくに、￣zrについてのforward-message -αおよびbackward-message -βを
【０１４３】
【数４３】
2019126425.tif 000045
【０１４４】
とすれば、(α,α*,β,β*)および(-α,-β)には
【０１４５】
【数４４】
2019126425.tif 000046
【０１４６】
という関係が成り立つ。
【０１４７】
以上のような関係性を利用してq(z)の更新を行う。まず、隠れマルコフモデルにおけるforward-backwardアルゴリズムを用いて-zrについてのメッセージ(-α,-β)を計算する。そして、その結果を用いてzrについてのメッセージ(α,α*,β,β*)を計算し、q(z)をそのメッセージに対応する分布に更新する。隠れマルコフモデルのforward-backwardアルゴリズムに必要な計算量はO((mJ)2R)であるため、この更新は時系列の長さRについて高々線形時間で済む。
【０１４８】
〔第３更新部１４０〕
第３更新部１４０は、q(b),q(z),ψが与えられているもとで、(19)式のL(q)を最大にするq(φ)を求め、q(φ)の値を求めた値で更新する(Ｓ１４０)。例えば、初回の処理では初期設定部１１０で設定したψと、それぞれ第１更新部１２０、第２更新部１３０で更新したq(b)、q(z)を用い、2回目以降の処理ではそれぞれ第１更新部１２０、第３更新部１４０、第４更新部１５０で更新した最新のq(b),q(φ),ψを用いる。
【０１４９】
q(b),q(z)およびψが与えられているもとで、L(q)を最大にするq(φ)は
【０１５０】
【数４５】
2019126425.tif 000047
【０１５１】
となる。
【０１５２】
以降、それぞれの項についての更新式を示す。まず、q(Ai),i= 1,…,Jを
【０１５３】
【数４６】
2019126425.tif 000048
【０１５４】
と更新する。ただし、
【０１５５】
【数４７】
2019126425.tif 000049
【０１５６】
である。
【０１５７】
また、q(π)を
【０１５８】
【数４８】
2019126425.tif 000050
【０１５９】
と更新する。ただし、
【０１６０】
【数４９】
2019126425.tif 000051
【０１６１】
である。
【０１６２】
さらに、q(θ)を
【０１６３】
【数５０】
2019126425.tif 000052
【０１６４】
と更新する。ただし、
【０１６５】
【数５１】
2019126425.tif 000053
【０１６６】
である。このようにして更新したq(Ai)、q(π)、q(θ)を用いて、(36)式によりq(φ)を更新する。
【０１６７】
〔第４更新部１５０〕
第４更新部１５０は、q(b),q(z),q(φ)を入力とし、q(b),q(z),q(φ)に基づいて、ハイパーパラメータU,V,Wを更新する(Ｓ１５０)。例えば、それぞれ第１更新部１２０、第２更新部１３０、第３更新部１４０で更新した最新のq(b),q(z),q(φ)を用いる。
【０１６８】
ハイパーパラメータそれぞれについての更新式を示す。まず、Uを
【０１６９】
【数５２】
2019126425.tif 000054
【０１７０】
により更新する。
【０１７１】
また、Vを
【０１７２】
【数５３】
2019126425.tif 000055
【０１７３】
により更新する。
【０１７４】
さらに、Wを
【０１７５】
【数５４】
2019126425.tif 000056
【０１７６】
により更新する。
【０１７７】
更新したU,V,Wを用いて、ΨとΩを
【０１７８】
【数５５】
2019126425.tif 000057
【０１７９】
と設定し、このΨとΩを用いた行列方程式
【０１８０】
【数５６】
2019126425.tif 000058
【０１８１】
の解fを求める。ただし、vecはベクトル化作用素であり、vec(Ω)はΩの各行ベクトルを並べたベクトルに対応する。この解fを用いて、F1,F2を
【０１８２】
【数５７】
2019126425.tif 000059
【０１８３】
により、更新する。
【０１８４】
以上の更新式は各パラメータについてL(q)を最大にするものとなっているが、外れ値に強いロバストな手法によって置き換えてもよい。
【０１８５】
以上の処理により、ハイパーパラメータψ=(U,V,W,F1,F2)を更新することができる。
【０１８６】
〔制御部１６０〕
制御部１６０は、所定の終了条件を満たすまで第１更新部１２０〜第４更新部１５０を繰り返し実行させる（Ｓ１６０）。例えば、予め定めた繰り返し回に達したことを終了条件とし、所定の繰り返し回数に到達するまで第１更新部１２０〜第４更新部１５０を繰り返し実行させるよう制御する。
【０１８７】
或いは、第１更新部１２０によりq(b)を更新する前のq(b),q(z),q(φ),ψに基づいて計算されるL(q)と、第１更新部１２０〜第４更新部１５０により更新された後のq(b),q(z),q(φ),ψに基づいて計算されるL(q)の差が所定の閾値以下となることを終了条件とし、それまで第１更新部１２０〜第４更新部１５０を繰り返し実行させるよう制御する。
【０１８８】
要するに、十分L(q)が大きくなるまで（L(q)が最大化に近づくまで）第１更新部１２０〜第４更新部１５０を繰り返し実行させればよい。
【０１８９】
〔出力部１７０〕
出力部１７０は、所定の終了条件を満たした時点のパラメータをモデルパラメータの推定結果q(b),q(z),q(φ),ψとして出力する（Ｓ１７０）。
【０１９０】
＜効果＞
以上の構成により、アイトラッカーによって計測された注視点の時系列からサリエンシーマップのモデルパラメータを簡便に推定することができる。
【０１９１】
また、変分事後分布qについての潜在変数の事後期待値や推定されたハイパーパラメータを観察することで、注視点の時系列に含まれる様々な情報を得ることができる。例えば、ar+srの事後期待値を計算することで注視点からノイズを除去した時系列が得られる。とくに、背後にある注視点のジャンプの性質を加味した雑音除去がなされているという点で、通常の平滑化に比べて有用である。また、Aの事後期待値から、被験者の注意遷移の振る舞いを観察できる。具体的には、Aの事後期待値についてマルコフクラスタリングアルゴリズムを適用することで、条件付き強度関数を表現するために用いたガウシアンカーネルをいくつかのクラスタに分割できる。分割されたカーネルのクラスタそれぞれをオブジェクトだとみなすことで、画像内に存在する注意を引く対象の数を同定することも可能である。図７に代表的なものについての解釈を示す。
【０１９２】
＜シミュレーション結果＞
参考文献６のデータセットに第一実施形態を適用した結果を示す。このデータは、被験者が画像を見ている間の3秒間の注視点の軌跡をアイトラッカーによって計測したものである。
（参考文献６）T. Judd, K. Ehinger, F. Durand and A. Torralba, "Learning to predict where humans look", IEEE International Conference on Computer Vision, 2009
【０１９３】
このデータセットに含まれるひとつの画像のデータについて、第一実施形態を適用した。  15人の被験者のデータのうち、4人の被験者の注視点の時系列データに対して本実施形態を適用した。図８にある被験者のデータの１反復目の推定結果を、図９に１１反復目の推定結果を示す。「Saccade delay」の縦方向の破線は、推定したサッカードの発生時刻を示す。最初の反復では、ノイズによって多くの偽のサッカードが推定されてしまっているものの、反復を繰り返すことによって真のサッカードのみを分離することができているのが見て取れる。また、図１０は推定に用いた4人を含む15人の被験者のデータから計算した注視点の頻度分布と、参考文献１のサリエンシーマップにサポートベクター回帰を適用して得られた初期値とを表す。図１１は１１反復目の推定結果を示す。本実施形態を適用することで、１１反復目の推定結果が初期値よりも注視点の頻度分布に近い画像となっているのが見て取れる。
【０１９４】
＜変形例＞
本実施形態では、q(b),q(z),q(φ),ψの順で、パラメータを更新しているが、更新の順番は変更してもよい。初期設定部１１０で設定した初期値、または、更新した最新のパラメータを用いて、第１更新部１２０、第２更新部１３０、第３更新部１４０、第４更新部１５０においてパラメータを更新すればよい。
【０１９５】
また、推定装置は、サッカード発生時刻推定部１８１、注視点系列推定部１８２、注目点系列推定部１８３、注目範囲推定部１８４、サリエンシーマップ生成部１８５と、の少なくともいずれかをさらに含む構成としてもよい（図１２参照）。
【０１９６】
例えば、サッカード発生時刻推定部１８１は、モデル推定部１０１で学習したモデルパラメータの推定結果q(b),q(z),q(φ),ψを入力とし、これらの値に基づいて、
【０１９７】
【数５８】
2019126425.tif 000060
【０１９８】
により、サッカードの発生時刻を推定し、出力する。
【０１９９】
注視点系列推定部１８２は、モデル推定部１０１で学習したモデルパラメータの推定結果q(b),q(z),q(φ),ψを入力とし、これらの値に基づいて、
【０２００】
【数５９】
2019126425.tif 000061
【０２０１】
により、注視点の時系列の推定結果を求め、出力する。
【０２０２】
注目点系列推定部１８３は、モデル推定部１０１で学習したモデルパラメータの推定結果q(b),q(z),q(φ),ψを入力とし、これらの値に基づいて、
【０２０３】
【数６０】
2019126425.tif 000062
【０２０４】
により、注目点の時系列の推定結果を求め、出力する。
【０２０５】
注目範囲推定部１８４は、モデル推定部１０１で学習したｚの事後分布q(z)に基づいて、注目範囲または注目対象の推定結果を求め、出力する。例えば、どの混合成分に向かってサッカードが発生したかを示すzの事後分布と、受容野の中心位置{μj} j=1Jと、受容野の広がり具合{Λj} j=1Jと、から、注目範囲または注目対象の推定結果を求めことができる。
【０２０６】
サリエンシーマップ生成部１８５は、モデル推定部１０１で学習したモデルパラメータの推定結果q(b),q(z),q(φ),ψを入力とし、これらの値に基づいて、
【０２０７】
【数６１】
2019126425.tif 000063
【０２０８】
により、サリエンシーマップを生成し、出力する。
【０２０９】
出力部１７０は、サリエンシーマップのモデルパラメータの推定結果、サッカードの発生時刻、注視点の時系列の推定結果、注目点の時系列の推定結果、注目範囲または注目対象の推定結果、サリエンシーマップ、の少なくとも何れかを入力とし、推定装置及び出力部１７０は、少なくとも何れかを出力する。
【０２１０】
＜その他の変形例＞
本発明は上記の実施形態及び変形例に限定されるものではない。例えば、上述の各種の処理は、記載に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。その他、本発明の趣旨を逸脱しない範囲で適宜変更が可能である。
【０２１１】
＜プログラム及び記録媒体＞
また、上記の実施形態及び変形例で説明した各装置における各種の処理機能をコンピュータによって実現してもよい。その場合、各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記各装置における各種の処理機能がコンピュータ上で実現される。
【０２１２】
この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。
【０２１３】
また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させてもよい。
【０２１４】
このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶部に格納する。そして、処理の実行時、このコンピュータは、自己の記憶部に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。また、このプログラムの別の実施形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよい。さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。なお、プログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。
【０２１５】
また、コンピュータ上で所定のプログラムを実行させることにより、各装置を構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。
【図１】
2019126425.tif 000068
【図２】
2019126425.tif 000069
【図３】
2019126425.tif 000070
【図４】
2019126425.tif 000071
【図５】
2019126425.tif 000072
【図６】
2019126425.tif 000073
【図７】
2019126425.tif 000074
【図８】
2019126425.tif 000075
【図９】
2019126425.tif 000076
【図１０】
2019126425.tif 000077
【図１１】
2019126425.tif 000078
【図１２】
2019126425.tif 000079