JP2019126425A - 推定装置、その方法、及びプログラム - Google Patents

Abstract

【課題】アイトラッカーによって計測された注視点の時系列から注視点の生成モデルのパラメータを推定する技術を提供する。【解決手段】時刻ｔの人の真の注目点は、時刻ｔ近傍にサッカードが発生する場合は時刻t-1の真の注目点を当該サッカードの方向及び大きさに応じて移動させた位置とし、時刻ｔ近傍にサッカードが発生しない場合は時刻t-1の真の注目点をランダムな方向及び大きさに応じて移動させた位置とし、時刻ｔの人の注視点の遅れは、AR(2)モデルに従うものとし、時刻ｔの人の真の注視点は、時刻ｔの真の注目点を時刻ｔの注視点の遅れにより補正した位置とし、実際に計測される時刻ｔの注視点を、上記時刻ｔの真の注視点にノイズが加わったものとしてモデル化したものを注視点の生成モデルとして、対象者の眼の動きを計測して得た注視点の時系列から、注視点の生成モデルのモデルパラメータを推定する。【選択図】図５

Description

本発明は、ヒトの注視点の時系列をモデル化した注視点モデルのパラメータを推定する技術に関する。

サリエンシーマップは、視覚刺激に対するボトムアップ性の注意を表現するための計算論的なモデルであり、元画像の色や輝度といった特徴から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に向かって発生したタイミングでそのニューロンが発火するイメージを示す図。 ニューロンの配置と受容野の中心位置との対応関係を示す図。 真の注視点を、真の注目点a_rと注視点の遅れs_rとに分解のイメージを示す図。 第一実施形態に係る推定装置の機能ブロック図。 第一実施形態に係る推定装置の処理フローの例を示す図。 事後期待値とその解釈の例を示す図。 被験者のデータの１反復目の推定結果を示す図。 被験者のデータの１１反復目の推定結果を示す図。 被験者のデータから計算した注視点の頻度分布と、サリエンシーマップにサポートベクター回帰を適用して得られた初期値とを表す図。 被験者のデータの１１反復目の推定結果を示す図。 第一実施形態の変形例に係る推定装置の機能ブロック図。

以下、本発明の実施形態について、説明する。なお、以下の説明に用いる図面では、同じ機能を持つ構成部や同じ処理を行うステップには同一の符号を記し、重複説明を省略する。以下の説明において、テキスト中で使用する記号「^」「~」「^-」等は、本来直後の文字の真上に記載されるべきものであるが、テキスト記法の制限により、当該文字の直前に記載する。式中においてはこれらの記号は本来の位置に記述している。また、ベクトルや行列の各要素単位で行われる処理は、特に断りが無い限り、そのベクトルやその行列の全ての要素に対して適用されるものとする。

＜サリエンシーマップのモデル＞
まず、サリエンシーマップのモデルについて説明する。

人がある一枚の画像を見ており、その間の注視点の時系列がアイトラッカーによって計測されている状況を想定する。本実施形態では、アイトラッカーによって計測された眼球が向いている位置を「注視点」、人が真に注目している位置を「注目点」と呼び分けることにする。本実施形態の目的は、画像情報を用いることで「注視点」の軌跡から眼球運動による遅れの成分やノイズ成分を除去し、「注目点」の軌跡を推定することである。この推定の際に、「注目点」の軌跡に影響を与える潜在変数としてサリエンシーマップが同時に計算される。

本実施形態では、注視点の時系列を生成モデルで表現し、その事後分布を推定する。

＜サリエンシーマップの条件付き強度関数による定義＞
従来のサリエンシーマップは、サッカードの行き先になりやすい地点で値が大きくなるように定義されている(参考文献１参照)。
（参考文献１）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.

そこで本実施形態では、このサリエンシーマップの考え方に則り、各地点のサリエンシーの値をその地点に向かうサッカードの発生頻度によって定義することにする。

サッカードは、眼球運動に含まれる跳躍性のジャンプのことであり、瞬間的に発生する運動である。そのため、その時系列は瞬間的に発生するイベントを並べたような時系列となる。統計的モデリングの枠組みにおいては、このようなランダムに生起するイベントの時系列を記述するための確率過程として点過程がよく用いられている。とくに各イベントが何らかの特徴量をもつ場合はマーク付き点過程が用いられる。マーク付き点過程とは、各イベントにマークとよばれる特徴量が付加された系列の確率的構造を記述するための確率過程であり、モデリングの対象としては地震や為替取引等が挙げられる。図１にマーク付き点過程から生成された標本の例を示す。図１では、横軸が時刻を表し、縦軸が所定のマークを表す。サッカードの時系列をマーク付き点過程とみなす場合は、サッカードの発生時刻をイベントの発生時刻、サッカードの振幅・方向・到着地点などの特徴量をマークとみなすことができる。

点過程モデリングを行う場合、データから点過程の確率的構造を推定することが目標となる。しかしながら、その確率密度関数を直接推定するのは難しい。これは、点過程が定義される距離空間の構造が複雑であり、尤度の計算が困難であることに起因する。したがって、点過程によって時系列データのモデリングを行う場合は、確率密度関数を推定するかわりに条件付き強度関数と呼ばれる関数を推定することが多い。

条件付き強度関数は、過去の履歴があたえられたもとでの単位時間あたりに発生するイベントの数の期待値に相当する関数である。具体的には

と定義される。ただし、Nは計数測度、κ∈Κは各イベントに付加されたマーク、H_tは時刻tまでの履歴である。このH_tは、t以前に発生したイベントの情報や同時に観測された共変量の値など、全ての情報を含んでもよい。

この条件付き強度関数によって、ほとんどの点過程の確率的構造は一意に決定されるため（例えば参考文献２参照）、点過程の推定を行う際は条件付き強度関数を推定することを考えれば十分である。
（参考文献２）D. J. Daley and D. Vere-Jones, "An Introduction to The Theory of Point Processes", Springer Science and Business Media, New York, 2003.

本実施形態では、サッカードのベクトルをマークとするようなマーク付き点過程からサッカードの系列が生成されていると仮定する。今回の問題設定では画像という時間的に定常な刺激を呈示しているため、サッカードの発生も時間的に定常であるとする。すなわち、条件付き強度関数が

と分解できるとし、サッカードの発生を示す関数λ(t|H_t)が定常であるとする。ただし、条件付き強度関数それ自体が確率過程である。このとき、λ_S(κ|H_t)はサッカードが画像平面上のどの地点に向かって発生しやすいかを説明する項であり、サリエンシーマップに対応する関数となる。

次に、λ_S(κ|H_t)としてどのような関数が適しているかを考察する。本実施形態では、サッカード生成に関連の深い中脳の上丘を参考にλ_S(κ|H_t)の形を決めることにする。上丘にはある特定の振幅・方向のサッカードが発生したときのみに発火するという性質をもつニューロンが存在する。すなわち、各ニューロンは視野上でmovement field(図２中、MFと表す)と呼ばれるある特定の領域をもち、サッカードがその領域に向かって発生したタイミングでそのニューロンが発火する。図２にそのイメージ図を示す。図中、spike trainはニューロンの発火状態を示し、scan pathは注視点を示す。また、それらのニューロンの配置(motor map on SC)は各々が反応する振幅・方向に空間的に対応したマップをなしていることが知られている。図３にその対応関係を示す。図３では上丘におけるニューロンの配置と受容野の中心位置は空間的に対応している。このような対応関係を根拠に、サリエンシーマップが上丘に表象されているといわれることもある。

サリエンシーマップが上丘によって表象されている場合、サリエンシーマップは上丘ニューロンの受容野(ニューロンの神経応答に変化を生じるような刺激が提示される空間の領域)を足し合わせた形で表現されていると考えられる。本実施形態では、この生理学的な知見を参考に、λ_S(κ|H_t)を次のようなガウシアンカーネルの非負結合によって表すことにする。

ただし、

であり、c_jは発火頻度、μ_jは受容野の中心位置、Λ_jは受容野の広がり具合を表すパラメータである。とくに、c_jについては過去のサッカードの履歴H_tに依存してよいとする。また、x_t-は時刻t直前での注視点である。時刻tでサッカードが発生した場合は、x_t-はサッカードの発生地点、κ+x_t-はサッカードの到着地点に対応する。

以上が本実施形態におけるサリエンシーマップの定義である。ここで、サリエンシーマップを条件付き強度関数によって定義することの利点を補足しておく。条件付き強度関数によってサリエンシーマップを定義する最大の利点は、サリエンシーの値がサッカードの頻度になっているという点である。したがって、Winner take allルールで説明できないような注視点の確率的な挙動を説明することができる。さらに、異なる画素間の値の比や、異なる画像における値の比に意味を与えることができる。例えば、A地点のサリエンシー値がB 地点のサリエンシー値の2倍である場合、A地点の方がB地点よりも2倍注視されやすいと解釈できる。さらに、サリエンシーマップを上丘ニューロンの発火活動の強度関数と直接対応付けることも可能である。

画像の情報については、次のように利用する。事前に、従来手法で計算されたサリエンシーマップ画像Sに対してガウシアンカーネルを用いたサポートベクター回帰により、{(c_j,μ_j,Λ_j)}^J _j=1を推定しておく。ただし、条件付き強度関数は非負である必要があるため、サポートベクター回帰を適用する際は重み係数^-c_jについて非負制約を導入する。事後分布の推定の段階では、事前に推定しておいたパラメータのうち{(μ_j,Λ_j)}^J _j=1を固定したハイパーパラメータとしてそのまま利用する。また、{c_j}^J _j=1は、c_j(H_t)を推定する際の事前分布として利用する。

＜注視点の生成モデル＞
観測時間をT、アイトラッカーのサンプリング回数をRとし、サンプリング間隔をΔt = T/Rとする。アイトラッカーによって計測された注視点をx_r∈R²,r=1,…,Rとし、人が真に着目している点の位置(つまり、注目点)をa_r∈R²,r=1,…,Rとする。この注目点a_rは、普段は連続的に変化しつつ、注意の対象が遷移した際に大きくジャンプするとする。本実施形態では、この瞬間的な注目点のジャンプをサッカードと呼ぶことにする。

ここで、κ_rを、時刻rにサッカードが発生した場合はそのサッカードのベクトル、発生しなかった場合は空集合{}(ただし、式中では、空集合を

とも表記する）をとる確率変数とする。このκ_rを用いて、注目点a_rが次のようなジャンプを伴うランダムウォークに従って生成されているとする。

このような仮定のもとで、注目点a_rは不連続な時系列となる。つまり、時刻r毎の人の真の注目点a_rは、１時刻前の真の注目点a_r-1を基準として、時刻rの近傍(時刻r-1〜時刻rの間)にサッカードが発生する場合はサッカードの方向及び大きさ(つまり、サッカードのベクトルκ_r)に応じて移動させた位置とし、時刻rの近傍にサッカードが発生しない場合はランダムな方向及び大きさへ移動させた位置となるようモデル化される。言い換えれば、時刻rの近傍にサッカードが発生しない場合はランダムウォークu_r〜N(0,U)により移動先の位置が生成されるものとしてモデル化される。

しかしながら、眼球の物理的な制約上、注視点は瞬間的に移動することはない。なぜならば、我々がある点を注目しようと思い立ったのちに注視点がその点に移動するまでには、筋肉が眼球を動かす分の遅れが生じるからである。そこで、この眼球運動の制約によって生じる注目点a_rからの注視点の遅れをs_rで表すことにし、時刻毎の人の注視点の遅れs_rは、次のようなAR(2)モデルによって生成されると仮定され、モデル化される。

ただし、F₁、F₂は眼球運動の特性を決めるパラメータである。

このもとで、真の注視点の軌跡はa_r+s_rによって表される。つまり、真の注視点は、真の注目点a_rを注視点の遅れs_rにより補正した位置としてモデル化される。図４に分解のイメージを示す。そして、実際に計測される注視点は、真の注視点a_r+s_rにノイズw_rが加わたものとして次のようにモデル化される。

次に、サッカード{κ_r}がどのように生成されるかを考える。本実施形態では、＜サリエンシーマップの条件付き強度関数による定義＞で説明したとおり、サッカードが(2)式および(3)式の条件付き強度関数によって定義されるマーク付き点過程から生成されていると仮定する。ただし、観測が離散的であるため、条件付き強度関数は離散時間で定義されることに注意する。いま、サッカードが発生した直後は、新たなサッカードは発生しにくく、かつそのサッカードの行き先は発生直前の注視点に依存すると考えられる。そこで、サッカードがサリエンシーマップのどの混合成分に向かって発生したものかを示す潜在変数

を導入し、{(κ_r,j_r)}が次のようなマーク付き点過程から生成されていると想定する。

ただし、r^*は時刻r以前で発生した最後(直近)のサッカードの発生時刻のインデックスであり、i^*はそのサッカードがどの混合成分に向かって発生したものかを示すインデックスである。r以前で発生したサッカードがない場合はr^*=0とする。A={A_ij}は確率推移行列であり、サッカードが発生したもとでそのサッカードが混合成分iから混合成分jに向かうものである確率がA_ijに対応する。また、π_jは最初に発生したサッカードが向かう混合成分がjである確率である。h(・)はサッカードが短い時間間隔で発生しないようにするための修正項であり、本実施形態では負の二項分布のハザード関数

を用いる。ただし、

であり、θは発生間隔の平均を、mはサッカード直後の発生しにくさを操作するパラメータである。とくに、m=1のとき負の二項分布は幾何分布と一致しハザードは定数となる。

さらに、A,πおよびθについては次のような事前分布を設定する。

ただし、

であり、{(c_j,μ_j,Λ_j)}は事前にサポートベクター回帰によって計算しておいたパラメータである。また、α₀は画像情報の事前分布への反映度合いを調節するパラメータであり、事前に適切な値を設定しておく。

以上が本実施形態における生成モデルの定義である。最後に、モデルをまとめておく。

＜事後分布の推定＞
推定のステップでは、(10)式のモデルのもとでの潜在変数の事後分布およびハイパーパラメータを推定することが目標となる。その推定の方法を順を追って説明する。まずはじめに、(10)式の生成モデルを扱いが容易なモデルに書き換える。次に、書き換えたモデルのもとで、事後分布の近似を求める方法を説明する。最後に得られた事後分布をどう解釈すればよいかを述べる。

＜スイッチング線形ガウス状態空間モデルへの帰着＞
上述の(10)式の生成モデルは、等価な隠れセミマルコフ・スイッチング線形ガウス状態空間モデルに変形することができる。このようなモデルに書き換えることで、変分ベイズの枠組みで潜在変数の事後分布を推定することが可能となる。

サッカードがどの混合成分に向けて発生したかを示す潜在変数{j_r}^R _r=1のかわりとして、新たな潜在変数{z_r}^R _r=1を

と定義する。ただし、j_(r+1)^*は時刻r+1以前に発生した最後のサッカードが、どの混合成分に向かって発生したものかを示すインデックスである。なお、上付き添え字及び下付き添え字におけるA^BはA^Bを意味し、A_BはA_Bを意味するものとする。また、{z_r}^R _r=1に関連する事象について

と表記することにする。以上の定義のもとで、{z_r}^R _r=1は推移確率が

であり、初期確率が

であるセミマルコフ過程に従う。

次に、{(x_r,a_r,s_r,κ_r)}^R _r=1の生成モデルを書き換える。書き換えの際のポイントはκ_rを消すことである。

ここで、新たな変数を以下のように定める。

このように変数を定義すると、{(x_r,b_r)}^R _r=1の生成モデルを次のような線形状態空間モデルに書き換えることができる。

したがって、(10)式の生成モデルは、背後にセミマルコフ過程に従う潜在変数が存在し、その潜在変数に従って局所的に線形ガウス状態空間モデルをつなげた形で表すことができる。このようなモデルをスイッチング線形状態空間モデルとよぶ。

最後に、(10)式を書き換えたモデルをまとめておく。

＜変分ベイズ＞
以降、表記の簡単のため

とする。

推定のステップでは、xが与えられたもとでのb,z,φの事後分布を計算することが目標となる。しかしながら、モデルの構造が複雑であるため、真の事後分布を解析的に計算することができない。そこで、本稿では変分ベイズの枠組みを用いて、真の事後分布の近似を計算することにする。

変分ベイズとは、計算が簡単になるようにあらかじめ指定しておいた分布族の中から、真の事後分布とのKLダイバージェンスが最小になるような分布を求め、その分布を真の事後分布の近似とする方法である。いま、xの周辺対数尤度log p(x|ψ)は、任意の分布q(b,z,φ)を用いて

と分解できる。ただし、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)が最大となる分布を求めることにする。すなわち、事後分布の推定を次のような最適化問題を解くことに帰着させる。

本実施形態ではこの分布族Qとして、独立性制約

を満たす分布の族を指定する。このような制約を導入することにより、座標降下法によって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）のような半導体メモリ素子により構成される補助記憶装置、またはリレーショナルデータベースやキーバリューストアなどのミドルウェアにより構成することができる。

図５は第一実施形態に係る推定装置の機能ブロック図を、図６はその処理フローの例を示す図である。

推定装置は、モデル推定部１０１と出力部１７０とを含む。

モデル推定部１０１は、初期設定部１１０と、第１更新部１２０と、第２更新部１３０と、第３更新部１４０と、第４更新部１５０と、制御部１６０とを含む。

<<モデル推定部１０１>>
[入力]：アイトラッカーにより計測された注視点の時系列x_r（r=1,2,…,R)
[出力]：モデルパラメータの推定結果q(b),q(z),q(φ),ψ
[処理]：モデル推定部１０１は、アイトラッカーにより計測された注視点(対象者の眼の動きを計測して得た注視点)の時系列x_r（r=1,2,…,R)を入力として、(17)式により表される注視点の時系列モデルの各パラメータb,z,φ,ψの事後分布を推定することにより、注視点の時系列モデル（学習済みモデル）を求める。

例えば、独立性制約

を満たす分布族の下で、

のL(q)が最大となる事後分布を求める。

以下、モデル推定部１０１の具体的な処理について説明する。

〔初期設定部１１０〕
初期設定部１１０は、q(b,z,φ),ψに適当な初期値を設定し(Ｓ１１０)、出力する。このとき、q(b,z,φ)については、(21)式

を満たすように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)は

となる。ここで、b以外の変数についての完全対数尤度の事後期待値は、次のような線形ガウス状態空間モデルの対数尤度と一致する。

ただし、~F_r,~o_r,~P_rは、

であり、{ξ_rj}_{r=1,…,R,j=1,…,J}は、

である。したがって、通常の線形ガウス状態空間モデルと同様、カルマンスムーザを適用することによって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)は

となる。ここで、z以外の変数についての完全対数尤度の事後期待値は、次のようなセミマルコフ過程の対数尤度と一致する。

ただし、

であり、これらの事後期待値は解析的に計算可能である。

したがって、隠れセミマルコフモデルの推定に用いられるforward-backwardアルゴリズムを適用することによってq(z)が計算可能である。forward-backwardアルゴリズムでは、次式で定義されるforward-message(α,α^*)およびbackward-message(β,β^*)を再起的に計算する。

ただし、o_rは時刻rでの観測に対応する変数である。計算されたメッセージを用いることで、z_rの事後分布を得ることができる。具体的な更新式の詳細については、参考文献３等を参照されたい。
(参考文献３)S. Z. Yu, "Hidden semi-Markov models. Articial intelligence", vol. 174, pp. 215-243, 2010.

しかしながら、通常のforward-backwardアルゴリズムに必要な計算量はO(J²R²)であり、時系列の長さRが長いときは望ましくない。そこで、本実施形態では参考文献４で提案された手法を用いて計算量を削減する。
(参考文献４)M. J. Johnson and A. S. Willsky, "Stochastic Variational Inference for Bayesian Time Series Models", International Conference on Machine Learning, 2014.

この手法の主なアイディアは、(27)式の隠れセミマルコフモデルを等価な隠れマルコフモデルへ変換するというものである。

ここで、新たな潜在変数^-z_rを

とし、

とする。ただし、

とする。以上の設定のもとで、対数尤度が

と表されるマルコフ過程を考えると、

となる（参考文献５参照）。
(参考文献５)M. J. Johnson, "Bayesian time series models and scalable inference", PhD thesis, Massachusetts Institute of Technology, 2014.

すなわち、{^-z_r}^R _r=1の分布が得られれば、{z_r}^R _r=1の分布も得られるということになる。とくに、￣z_rについてのforward-message ^-αおよびbackward-message ^-βを

とすれば、(α,α^*,β,β^*)および(^-α,^-β)には

という関係が成り立つ。

以上のような関係性を利用してq(z)の更新を行う。まず、隠れマルコフモデルにおけるforward-backwardアルゴリズムを用いて^-z_rについてのメッセージ(^-α,^-β)を計算する。そして、その結果を用いてz_rについてのメッセージ(α,α^*,β,β^*)を計算し、q(z)をそのメッセージに対応する分布に更新する。隠れマルコフモデルのforward-backwardアルゴリズムに必要な計算量はO((mJ)²R)であるため、この更新は時系列の長さ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(φ)は

となる。

以降、それぞれの項についての更新式を示す。まず、q(A_i),i= 1,…,Jを

と更新する。ただし、

である。

また、q(π)を

と更新する。ただし、

である。

さらに、q(θ)を

と更新する。ただし、

である。このようにして更新したq(A_i)、q(π)、q(θ)を用いて、(36)式によりq(φ)を更新する。

〔第４更新部１５０〕
第４更新部１５０は、q(b),q(z),q(φ)を入力とし、q(b),q(z),q(φ)に基づいて、ハイパーパラメータU,V,Wを更新する(Ｓ１５０)。例えば、それぞれ第１更新部１２０、第２更新部１３０、第３更新部１４０で更新した最新のq(b),q(z),q(φ)を用いる。

ハイパーパラメータそれぞれについての更新式を示す。まず、Uを

により更新する。

また、Vを

により更新する。

さらに、Wを

により更新する。

更新したU,V,Wを用いて、ΨとΩを

と設定し、このΨとΩを用いた行列方程式

の解fを求める。ただし、vecはベクトル化作用素であり、vec(Ω)はΩの各行ベクトルを並べたベクトルに対応する。この解fを用いて、F₁,F₂を

により、更新する。

以上の更新式は各パラメータについてL(q)を最大にするものとなっているが、外れ値に強いロバストな手法によって置き換えてもよい。

以上の処理により、ハイパーパラメータψ=(U,V,W,F₁,F₂)を更新することができる。

〔制御部１６０〕
制御部１６０は、所定の終了条件を満たすまで第１更新部１２０〜第４更新部１５０を繰り返し実行させる（Ｓ１６０）。例えば、予め定めた繰り返し回に達したことを終了条件とし、所定の繰り返し回数に到達するまで第１更新部１２０〜第４更新部１５０を繰り返し実行させるよう制御する。

或いは、第１更新部１２０により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についての潜在変数の事後期待値や推定されたハイパーパラメータを観察することで、注視点の時系列に含まれる様々な情報を得ることができる。例えば、a_r+s_rの事後期待値を計算することで注視点からノイズを除去した時系列が得られる。とくに、背後にある注視点のジャンプの性質を加味した雑音除去がなされているという点で、通常の平滑化に比べて有用である。また、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(φ),ψを入力とし、これらの値に基づいて、

により、サッカードの発生時刻を推定し、出力する。

注視点系列推定部１８２は、モデル推定部１０１で学習したモデルパラメータの推定結果q(b),q(z),q(φ),ψを入力とし、これらの値に基づいて、

により、注視点の時系列の推定結果を求め、出力する。

注目点系列推定部１８３は、モデル推定部１０１で学習したモデルパラメータの推定結果q(b),q(z),q(φ),ψを入力とし、これらの値に基づいて、

により、注目点の時系列の推定結果を求め、出力する。

注目範囲推定部１８４は、モデル推定部１０１で学習したｚの事後分布q(z)に基づいて、注目範囲または注目対象の推定結果を求め、出力する。例えば、どの混合成分に向かってサッカードが発生したかを示すzの事後分布と、受容野の中心位置{μ_j} _j=1 ^Jと、受容野の広がり具合{Λ_j} _j=1 ^Jと、から、注目範囲または注目対象の推定結果を求めことができる。

サリエンシーマップ生成部１８５は、モデル推定部１０１で学習したモデルパラメータの推定結果q(b),q(z),q(φ),ψを入力とし、これらの値に基づいて、

により、サリエンシーマップを生成し、出力する。

出力部１７０は、サリエンシーマップのモデルパラメータの推定結果、サッカードの発生時刻、注視点の時系列の推定結果、注目点の時系列の推定結果、注目範囲または注目対象の推定結果、サリエンシーマップ、の少なくとも何れかを入力とし、推定装置及び出力部１７０は、少なくとも何れかを出力する。

＜その他の変形例＞
本発明は上記の実施形態及び変形例に限定されるものではない。例えば、上述の各種の処理は、記載に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。その他、本発明の趣旨を逸脱しない範囲で適宜変更が可能である。

＜プログラム及び記録媒体＞
また、上記の実施形態及び変形例で説明した各装置における各種の処理機能をコンピュータによって実現してもよい。その場合、各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記各装置における各種の処理機能がコンピュータ上で実現される。

この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。

また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させてもよい。

このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶部に格納する。そして、処理の実行時、このコンピュータは、自己の記憶部に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。また、このプログラムの別の実施形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよい。さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。なお、プログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。

また、コンピュータ上で所定のプログラムを実行させることにより、各装置を構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。

Claims (7)

1. サッカードの時系列は、各サッカードの特徴量をマークとするマーク付き点過程で生成されるものとし、
   時刻ｔの人の真の注目点は、時刻ｔ近傍にサッカードが発生する場合は時刻t-1の真の注目点を当該サッカードの方向及び大きさに応じて移動させた位置とし、時刻ｔ近傍にサッカードが発生しない場合は時刻t-1の真の注目点をランダムな方向及び大きさに応じて移動させた位置とし、
   時刻ｔの人の注視点の遅れは、AR(2)モデルに従うものとし、
   時刻ｔの人の真の注視点は、前記時刻ｔの真の注目点を前記時刻ｔの注視点の遅れにより補正した位置とし、
   実際に計測される時刻ｔの注視点を、上記時刻ｔの真の注視点にノイズが加わったものとしてモデル化したものを注視点の生成モデルとして、
   対象者の眼の動きを計測して得た注視点の時系列から、前記注視点の生成モデルのモデルパラメータを推定するモデル推定部を含む、
   推定装置。
2. 請求項１記載の推定装置であって、
   前記モデル推定部で学習したモデルパラメータに基づいて、サッカードの発生時刻を推定するサッカード発生時刻推定部と、
   前記モデル推定部で学習したモデルパラメータに基づいて、注視点の時系列の推定結果を求める注視点系列推定部と、
   前記モデル推定部で学習したモデルパラメータに基づいて、注目点の時系列の推定結果を求める注目点系列推定部と、
   前記モデル推定部で学習したモデルパラメータに基づいて、注目範囲または注目対象の推定結果を求める注目範囲推定部と、
   前記モデル推定部で学習したモデルパラメータに基づいて、サリエンシーマップを生成するサリエンシーマップ生成部と、の少なくともいずれかをさらに含む、
   推定装置。
3. 請求項１記載の推定装置であって、
   κ_rを、時刻rにサッカードが発生した場合はそのサッカードのベクトル、発生しなかった場合は空集合をとる確率変数とし、j_rを時刻rにおけるサッカードがサリエンシーマップのどの混合成分に向かって発生したものかを示す潜在変数とし、H_rを時刻rまでの履歴とし、Sをサリエンシーマップ画像とし、r^*を時刻r以前で発生した最後のサッカードの発生時刻とし、時刻r以前に発生したサッカードがない場合はr^*=0とし、a_rを時刻rにおける真の注目点とし、s_rを時刻rにおける注視点の遅れとし、h_NBをサッカードが短い時間間隔で発生しないようにするための修正項とし、i^*を時刻r以前で発生した最後のサッカードがどの混合成分に向かって発生したものかを示すものとし、
   前記生成モデルは、
   
   により与えられる、
   ことを特徴とする推定装置。
4. 請求項３記載の推定装置であって、
   j_(r+1)^*を時刻r+1以前に発生した最後のサッカードが、どの混合成分に向かって発生したものかを示すものとし、
   前記生成モデルは、
   
   により与えられる、
   ことを特徴とする推定装置。
5. 請求項４記載の推定装置であって、
   
   とし、q(b,z,φ)を任意の事後分布とし、p(b,z,φ|x,ψ)を真の事後分布とし、前記モデル推定部は、次式により与えられるL(q)が大きくなるように、事後分布q(b),q(z),q(φ),ψを繰り返し更新することにより前記モデルパラメータを学習する、
   
   ことを特徴とする推定装置。
6. サッカードの時系列は、各サッカードの特徴量をマークとするマーク付き点過程で生成されるものとし、
   時刻ｔの人の真の注目点は、時刻ｔ近傍にサッカードが発生する場合は時刻t-1の真の注目点を当該サッカードの方向及び大きさに応じて移動させた位置とし、時刻ｔ近傍にサッカードが発生しない場合は時刻t-1の真の注目点をランダムな方向及び大きさに応じて移動させた位置とし、
   時刻ｔの人の注視点の遅れは、AR(2)モデルに従うものとし、
   時刻ｔの人の真の注視点は、前記時刻ｔの真の注目点を前記時刻ｔの注視点の遅れにより補正した位置とし、
   実際に計測される時刻ｔの注視点を、上記時刻ｔの真の注視点にノイズが加わったものとしてモデル化したものを注視点の生成モデルとして、
   対象者の眼の動きを計測して得た注視点の時系列から、前記注視点の生成モデルのモデルパラメータを推定するモデル推定ステップを含む、
   推定方法。
7. 請求項１から請求項５の何れかの推定装置としてコンピュータを機能させるためのプログラム。