社団法人 電子情報通信学会
THE INSTITUTE OF ELECTRONICS,
INFORMATION AND COMMUNICATION ENGINEERS

信学技報
TECHNICAL REPORT OF IEICE.

局所的変分法による非補償型時系列 IRT
玉野 浩嗣†,†† 持橋 大地††† NEC データサイエンス研究所 〒211–8666 神奈川県川崎市中原区下沼部 1753†† 統計数理研究所 〒190–8562 東京都立川市緑町 10–3
E-mail: †h-tamano@nec.com, ††daichi@ism.ac.jp
あらまし

項目応答理論 (IRT) の時間発展版にあたるナレッジトレーシングは，学習者の潜在スキル状態の時間変化を推定し，各問題に正解できるかどうかを予測するモデルである．教育ドメインの特性上，ナレッジトレーシングは予測精度以外に，予測モデルの解釈性や，予測理由の説明性が重要となる．これには人間が理解可能なスキル名を使った説明が必要不可欠である．スキルをある問題に関連付ける場合，それらスキルがすべて満たされることでその問題が解けるとすることが一般的である．このようなモデルは多次元項目応答理論において非補償型モデルと呼ばれ，非補償型モデルを使った予測理由の説明は自然である．そこで，本研究では非補償型の項目応答モデルを持つ確率モデ
ルによる，解釈可能なナレッジトレーシングを提案する．提案モデルは線形動的システムをベースとし，出力確率に非補償型の項目応答モデルを用いる．状態潜在変数の事後確率が解析的に得られなくなる問題に対し，局所的変分法
により事後確率をガウス近似する方法を提案する．人工データにより提案方法が真の事後分布を近似できていることを確認し，ASSISTment Data により深層学習ベースの手法に対しても高い予測性能を示すことを確認する．
キーワード 教育工学，項目応答理論，ナレッジトレーシング，カルマンフィルター，変分近似

Non-compensatory Temporal IRT with Local Variational Approximation
Hiroshi TAMANO†,†† and Daichi MOCHIHASHI††† NEC Data Science Laboratories 1753 Shimonumabe, Nakahara-ku, Kawasaki-shi 211–8666 Japan†† The Institute of Statistical Mathematics 10–3 Midori-cho, Tachikawa-shi 190–8562 Japan
E-mail: †h-tamano@nec.com, ††daichi@ism.ac.jp
Abstract Knowledge tracing, the time varying extension of item response theory (IRT), is a task to trace learner’s
latent skill states to predict whether the learner can answer a new question correctly. Due to its educational domain,knowledge tracing needs high interpretability of its prediction and explainability of its result. As for explainability,
explanation with human interpretable skills is necessary, which are usually given as conjunctive conditions. Such amodel is called a non-compensatory model in multi-dimensional IRT and explanation using non-compensatory item
response model is desired. To realize an interpretable and explainable knowledge tracing, we propose a probabilistic
model based on non-compensatory item response model combined with a linear dynamical system. Since it results
in a complicated posterior on the skill states of the learners, we approximate it using a local variational distribution.
We also show that our posterior adequately approximates the true posterior in artificial data, and our prediction
performance is better than two popular deep learning based knowledge tracing in ASSISTment Data.
Key words Educational Technology, IRT, Knowledge Tracing, Kalman Filter, Variational Approximation

1.

まえがき

用したデジタル教材や，Coursera や Udacity などの Massive

Open Online Courses (MOOCs) によるオンライン学習によ教育のデジタル化により学習履歴データが蓄積され，データ

り，学習者がいつ，どのコンテンツをどれくらい学習したかなに基づき個々の学習者に適した学びを提供するアダプティブ

どの様々なログが蓄積され，テスト結果との関連の分析などがラーニングシステムの利用が広がっている [1] [2]．教育のデジ進められている [3]．また，アダプティブラーニングとしては，

タル化では，PC，タブレット PC，スマートフォンなどを利米 Knewton 社 [1] に代表されるように，学習履歴から個々に最—1—

適な課題のレコメンデーションが提示され，学習効果を上げている．

そこで，本研究では非補償型の項目応答モデルを持つ確率モデルによる解釈可能なナレッジトレーシングを提案する．提案

アダプティブラーニングシステムにおいて，学習者のスキルモデルは，線形動的システムをベースとして出力確率が非補償状態を追跡し，個々に適した学習コンテンツを推薦するために型の項目応答モデルである．線形動的システムにおける状態の中心的な役割を果たす技術がナレッジトレーシングである．ナ

事後分布はガウス分布であるため解析的に計算可能であるが，レッジトレーシングは，学習者が問題を解いた時系列正誤ログ

非補償型の項目応答モデルを導入することで，ガウス分布ではを入力として，次のタイムステップにおいて，ある問題を正解

なくなり解析的に計算ができなくなる．そこで，非補償型の尤

できるかどうかを予測する．項目応答理論の時間発展版と見る度関数を局所的変分法によりガウス近似することで，事後分布

こともできる．代表的な手法としては，二値分類問題としてモを解析的にガウス近似する．事後分布がガウス分布で近似でき

デル化する手法 [4] [5]，隠れマルコフモデルをベースとした時ることにより，モデルパラメータの学習は，線形動的システム

系列の確率モデルによる手法 [6]，Recurrent Neural Networksとほぼ同様に行うことができ，非補償型項目応答モデルのパラ(RNNs) をベースとした時系列の深層学習モデル [7] [8] が提案

メータの最適化のみモンテカルロ EM [11] を用いる．

されている．

以下，2 章で提案手法の基礎となる確率モデルとして，多次またナレッジトレーシングでは，教育ドメインの特性上，予元項目応答理論と線形動的システムについて説明する．次に，測精度以外に予測モデルの解釈性、予測理由の説明性が重要と3 章にて，提案手法である非補償型の時系列 IRT について説明なる．予測モデルの解釈性は，予測モデルが人間に理解可能か

し，4 章にて，提案手法について事後分布の近似に関する実験

どうかである．深層学習モデルは非常に複雑で人間が理解する

結果，ASSISTments データによる予測精度の実験結果についことは難しく，解釈性は低いとされる．こうした挙動が完全にて示す．5 章にて，関連研究を述べ，6 章にてまとめる．

はわからない予測モデルは，問題が発生しても対処できないた
め学校教育に利用することは難しい．また，予測理由の説明性
は，問題の正誤予測においてその理由を学習者または教師が理解可能な形でフィードバックできるかどうかであり，予測理由を提示することで，学習者は自分に何が足らないものを認識で

2.

準

備

まず，提案手法のベースとなる多次元項目応答理論と線形動的システムについて説明する．

2. 1.

多次元項目応答理論

き，または教師がその子に何が足りないかを認識でき，その後

項目応答理論は，各問題を各学習者が正解したかどうかのテの学習や指導に役立てることができる．さらに，予測には不確

スト結果のデータから，各問題の難易度，各学習者の能力の推実性が伴うため，予測の不確かさについても学習者または教師

定などに用いられる．項目応答理論では通常，学習者の能力をにフィードバックすることで，どの程度フィードバックを信じ

一次元の値として扱うが，多次元項目応答理論では多次元の値るべきか判断することができるため有用である．

として扱う．例えば，複数のスキルを必要とする問題で構成さ

予測理由の説明には，人間が理解可能なスキル名を使った説
明が必要不可欠であり，多次元項目応答理論における非補償型モ

れるテストの結果から，学習者の各スキルの習熟度を測定する
場合などで多次元項目応答理論が利用される．

3
デルを使った説明が適している．例えば，方程式 51 x + 10
= 2x

多次元項目項目応答理論は，補償型 (compensatory) と非補を解くためには，分数のスキルと方程式のスキルの両方が満た

償型 (non-compensatory) のモデルに大別することができる．

されないと解くことはできない．このように問題を解くために補償型モデルでは，問題を解くためには複数スキルの総和が高必要なスキルは，多くの場合 AND 条件で記述されており，こければ正解できるという仮定をおく．一方で，非補償型モデルのようなモデルは多次元項目応答理論において非補償型モデル

では，問題を解くためにはすべてのスキルが高くなければ正解と呼ばれている．学習者の潜在スキルの時間変化がない，項目できないという仮定をおく．どちらのモデルを利用するかは，応答データの確率モデルについては，[9] [10] が提案されている

そのデータの生成がどちらの仮定の方が適切かで分析者が選択

が，学習者の潜在スキルが時間変化するナレッジトレーシング

する．

において，非補償型モデルを使った確率モデルと，時間変化する学習者の潜在スキルの事後分布を推定する方法は我々の知る補償型モデルでは，問題 i を学習者 j が正解するかどうかを確率変数 yi,j ∈ {0, 1} で表し，正解する確率を限り提案されていない．これは，ナレッジトレーシング研究がこれまで予測にフォーカスしてきたためと考えられる．学習者
へ予測理由を説明することを考えなければ，学習者の潜在ス
キル状態を低次元空間に埋め込みロジスティック回帰する方法 [6] [7] [8] でも十分である．また，[9] [10] はサンプリングを用いた推定方法であるため，学習者の潜在スキルをリアルタイムに推定できることが望ましいナレッジトレーシングにおいては直接応用することは適さず，また解析的に事後分布を求めよう
としても非補償型モデルの尤度関数が複雑であるため解析解はp(yi,j = 1|zj , ai , bi ) =1
1 + exp(−(aTi zj − bi ))(1)のように定義する．ここで，zj は学習者 i の各スキルの習熟度を表すベクトル，ai は問題 i における各スキルの識別力を表すベクトル，bi は問題 i の難易度である．補償型では，各スキルの習熟度が足し合わされるため，各スキルはお互いに補完できるとするモデルになっている．
非補償型モデルでは，問題 i を学習者 j が正解する確率を求まらない．

—2—

z (t−1)…

z (t)式（7）から α̂(z (t−1) ) がガウス分布であれば，α̂(z (t) ) もガウス

z (t+1)分布になることが分かる．t = 1 の場合，α̂(z (1) ) = p(z (1) |y (1) )…

がガウス分布であることから，すべてのフォワード再帰はガウス分布となる．バックワード再帰についても，同様にガウス分
布で計算することができ，こちらは γ メッセージと呼ばれる．

y (t−1) y (t)LDS のパラメータの学習について述べる．LDS のパラメー

y (t+1)

タは，θ = {D, Γ, C, Σ, µ0 , P0 } であり，EM アルゴリズムを用
いて学習することができる．E ステップでは，フォワード再帰

図 1: LDS のグラフィカルモデル．

p(yi,j = 1|zj , ai , bi ) =∏kとバックワード再帰を用いて現在のパラメータにおける潜在

1
1 + exp(−(ai,k zj,k − bi,k ))(2)のように定義する．ここで，k はスキルのインデックスを表

変数の事後分布 p(z (t) |y (1) , . . . , y (T ) ; θ old ) を算出する．M ス
テップでは，E ステップで求めた事後分布を用いて，完全対数尤度の期待値

し，問題の難易度パラメータ bi は補償型の場合とは異なり，各Q(θ, θ old ) = E[ln p(y (1:T ) , z (1:T ) ; θ)]

スキル毎の値を持つベクトルとなっている．非補償型では，各(9)スキルの習熟度が掛け合わされるため，各スキルはお互いに補を，パラメータ θ について最大化する．上記の線形動的システ

完できないモデルになっている．

ムではすべてのパラメータについて解析解が求まる．

線形動的システム

2. 2.

線形動的システム (linear dynamic system, LDS) [12] は，状態が時間変化しているときに各時刻 t の状態 z (t) に応じてデータy(t)が観測される場合に用いられるモデルであり，ロボット

非補償型時系列 IRT

3.

本章では，提案手法である非補償型時系列 IRT について説明する．非補償型時系列 IRT は，線形動的システムの出力確率がガウス分布ではなく，非補償型項目応答モデルの尤度関数にの位置推定などに利用することができる．

LDS の確率モデルは，遷移確率と出力確率と初期状態確率になっている．この変更により α̂ メッセージと γ メッセージがガウス分布では表現できなくなる．そこで本研究では，局所的より変分法 [13] により，非補償型項目応答モデルの尤度関数をガウp(z (t) |z (t−1) ) = N (z (t) |Dz (t−1) , Γ)(3)p(y (t) |z (t) ) = N (y (t) |Cz (t) , Σ)(4)度関数の変分下界を指数 2 次形式で求め，変分パラメータを最) = N (z(5)

適化することでガウス近似を得る．尤度関数がガウス近似され

p(z

(1)

(1)

|µ0 , P0 )と記述される．初期状態は平均 µ0 ，共分散 P0 のガウス分布に従って発生し，各時刻における状態は線形変換 D と平均 0，共分散 Γ のガウスノイズにより遷移する．データは，状態に線形変換 C と平均 0，共分散 Σ のガウスノイズが加わり観測される．また，LDS のグラフィカルモデルは図 1 で表される．観測

ス分布で近似する．具体的には，非補償型項目応答モデルの尤ることで，線形動的システムと同様の方法で潜在状態の事後分
布とパラメータの学習が可能になる．以下ではまず，非補償型
時系列 IRT の生成モデルについて説明し，次に非補償型の尤度関数を局所ガウス近似する方法（指数 2 次形式で変分下界を求める方法）を説明し，最後に変分パラメータの最適化方法を説明する．モデルパラメータの学習は紙面の都合上省略するが，変数は黒丸，潜在変数は白丸で表す．

LDS において観測データが与えられた場合の潜在変数の事
後確率の算出方法について述べる．隠れマルコフモデルと同様

出力確率のパラメータの最適化のみモンテカルロ EM [11] を使
い，他は線形動的システムと同様に行えばよい．
生成モデルにフォワード再帰とバックワード再帰と呼ばれる式を順方向，

3. 1.

逆方向に適用することで事後確率が算出でき，前者はカルマン

非補償型時系列 IRT の生成モデルは，線形動的システムと同フィルタ (Kalman filter)，後者はカルマンスムーザ (Kalmansmoother) と呼ばれる．

(t+1)

フォワード再帰では，データ yでの潜在変数 z(t)

様に状態遷移確率，出力確率，初期状態確率により次の(1),...,y(t)p(zjが与えられた上

(6)





..
.










 (t+1)(t),ΓDi(j,t) zj + βkT x(t+1)
 (10)
i(j,t+1)
j,k




.

の事後確率

α̂(z (t) ) := p(z (t) |y (1) , . . . , y (t) )



(t)

|zj )

= N zj..

を t = 1 から順に t = T まで算出する．これを α̂ メッセージと呼ぶ．α̂ メッセージの計算には次の再帰式を利用する．

α̂(z (t) )
∝ p(y

(t)

(t)

(t)

p(yj = 1|zj ) =
|z

)

(t)

σ ai(j,t),k (zj,k − bi(j,t),k ))Qi(j,t),k

(11)

k

∫
(t)

∏ (

(1)

p(z (t) |z (t−1) )α̂(z (t−1) )dz (t−1)

(7)

(t)

(8)

(1)

(t−1)

∝ p(y (t) |z (t) )p(zj |yj , . . . , yj)

(1)p(zj ) = N (zj |µ0 , P0 )(12)で 定 義 さ れ る ．こ こ で ，学 習 者 j ∈ {1, · · · , N } の 時 刻

—3—

!"#%
%&'(
)*+

,./0

'"

#

$

!
!

図 2: 正解側の尤度関数

./1

"

&%
図 3: 不正解側の尤度関数

!"#$

&$

(t)
t ∈ {1, · · · , Tj } のスキル状態を確率変数 zj ∈ RK で表し，学習者 j が時刻 t に取り組んだ問題に正解できたかどうかを確率変
(t)

数 yj

∈ {0, 1} で表す．Qi,k ∈ {0, 1} は問題とスキルの対応表

で，問題 i ∈ {1, · · · , M } を解くためにスキル k ∈ {1, · · · , K}が必要であれば 1，不要であれば 0 である．問題とスキルの対
応表 Q はデータとして与えられているとする．i(j, t) は，学習者 i が時刻 t に取り組んだ問題インデックスを表し，以降で文(t)
脈から明らかな場合には省略する場合がある．xj,k ∈ RFk は状態遷移の共変量であり任意の情報を入れ込むことができる．例
えば，学習者の属性 (e.g., 年齢) や学習者 j がスキル k を含む問題を最後に学習した時からの経過時間から考えられる忘却率∆

2− h (∆: 経過時間,h:半減期) などである．また，問題ごとのバ(t)イアス項も xj,k に次元を用意し，1 を立てることで表現することができる．Di , Γi は共に対角行列で，Qi,k = 1 である k に対しては (Di )k,k = di,k , (Γi )k,k = γk であり，Qi,k = 0 である kに対しては，(Di )k,k = 1, (Γi )k,k = 0 とする（注 1）．

3. 2.

局所ガウス近似!"#$

図 4: 変分下界導出ステップ 1!
算出し(図 4 左)，尤度関数を勾配方向で 1 次元化する (図 4 右の尤度関数)．1 次元化された尤度関数に対して，指数 2 次下界を導出する (図 4 右の下界 2)．ここまでで，ガウス近似の中心ベクトルと勾配方向の分散までが得られたことになる．次に，ステップ 2 では勾配方向に直交する他の成分についの共分散行列を考える．ここで，尤度関数の等高線が凸関数であることに注意すると，勾配方向に直交する方向では尤度関数は上昇する．
そのため，どのような共分散行列を用いても下界となる. ここでは事前分布の勾配方向に直交する成分の共分散行列を利用する．(図 5 左) 以下，ステップ 1 とステップ 2 について順に具体的に下界を導出する．
ステップ 1 では，まずシグモイド関数の変分下界 [13] を利用できるように不正解尤度をシグモイド関数の積の和の形式に変換する．

p(y = 0|za ) = 1 −∏

本節では，不正解側，正解側の順に尤度関数の局所ガウス近k

似方法を説明する．本節と次節ではある学習者 j が時刻 t に問

2Ki −1=

題 i(j, t) を解いた場合を考え，煩雑さを避けるためにインデッ

σ (ai,k (za,k −bi,k ))∑∏l=1

クス j, t を省略する．また，スキルインデックス k は問題 i にσ (sign(l, k) (ai,k (za,k −bi,k ))) .(13)k

関係する Qi,k = 1 であるスキルのみを対象にする．学習者のス

ここで，sign(l, k) は l を二進数表記にしたときの k ビット目がキル状態には問題 i に関係する部分を切り出した変数 za ∈ RKi

0 であれば −1 で 1 であれば +1 とする．次に，z 座標系から勾を使う．Ki は，問題 i が必要とするスキル数である．za を使い，式 (11) は p(y = 1|za ) =

∏K i

k=1

σ(ai,k (za,k − bi,k )) と表現

配方向に軸を取り直した s 座標系 (図 4 左) に変換し，s1 軸方
向に 1 次元化する．z 座標系と s 座標系の変換を

(

できる．

3. 2. 1

不正解側尤度のガウス近似

不正解側尤度関数は p(y = 0|za ) = 1−za = W∏ks+

σ(ai,k (za,k −bi,k ))

[

])
0

W T ma(14)であり，図 3 の形をしている．形がガウス分布と大きく異るたで定義する．行列 W の各列は正規直交基底をなし，1 列目のベ

め，全体をガウス分布で近似しても良い α̂ メッセージの近似にクトルは勾配 ∇p(y = 0|za )|ma 方向とする．また，W T を W Tはならない．そこで，変分ロジスティック回帰 [13] で利用されから 1 行目を除いたものとする．式 (14) に s = [s1 , 0, . . . , 0]Tている局所的変分法を用いる．具体的には，尤度関数の変分下を代入することで勾配方向に 1 次元化した場合の za,k と s1 の界を指数 2 次形式で求めることで尤度関数をガウス近似する．

関係式が次の[変分下界に含まれる変分パラメータは周辺尤度を最大化することにより決定する．周辺尤度を最大化するように変分下界を最za,k = wk,:

適化するため，事前分布付近で尤度関数を局所的にガウス近似]
s1

W T ma= wk,1 s1 + wk,2: W T ma

(15)

することになる．本節では，変分パラメータの値は既知とし，で表せる．wk,: は，行列 W の k 行目の行ベクトルを表し，wk,2:

最適な変分パラメータの求め方は 3. 3 節で説明する．

を行列 W の k 行目の 2 列目以降の行ベクトルを表す．これを変分下界を求める手順は次の 2 つのステップからなる．ま式 (13) に代入すると，勾配方向に 1 次元化した尤度ず，ステップ 1 では事前分布の中心における尤度関数の勾配を=（注 1）：決定的な遷移を表現するために便宜的に分散を 0 と表現した．

∑ ∏ (2Ki −1(σ sign(l, k) ai,k (wk,1 s1 + wk,2: W T ma − bi,k )

l=1

k

—4—

))

∑ ∏

2Ki −1

=

l=1

!"#%

σ(ãl,i,k s1 − b̃l,i,k )

,-

(16)

k

が得られる．ただし，!

ãl,i,k = sign(l, k)ai,k wk,1%&'(b̃l,i,k = −sign(l, k)ai,k (wk,2: W T ma − bi,k )

&$

!"#$

とする．ここでシグモイド関数の変分下界

σ(z) ≥ σ(ξ) exp(

図 5: 変分下界導出ステップ 2)

1(z − ξ) − λ(ξ)(z 2 − ξ 2 )
2

(17)

∑

≥

と計算できる．以上により，s 座標系における尤度関数のガウス近似は N (s|ηs , Ψs ) となる．

を適用すると2Ki −1

(
∏

l=1

∑

q B
∑l l l[ (∑)
 2 l q l Al 
2
qA

 −1
l l l
ηs = 
,
Ψ
=

s
0


0

)
σ(ξl,k )exp(−Al s21 + Bl s1 + Cl )

(18)

k

..
.

(:= lb1 (s1 , ξ))となる．（図 4 右の下界 1）ただし，ξl,k = ãl,i,k ξ − b̃l,i,k ，1
λ(ξ) = 2ξ
(σ(ξ) − 21 )，

Al =

∑

Bl =

k

Cl =

∑

1
2λ(ξl,k )ãl,i,k b̃l,i,k + ãl,i,k
2

∑

[ (∑
)
q B
∑l l l
2
qA
−1
l l l
η = W  2 l ql Al  , Ψ = W
W ma

1
2
λ(ξl,k )(ξl,k− b̃2l,i,k ) − (ξl,k + b̃l,i,k )2

3. 2. 2である．ξ は変分パラメータである．さらに，対数を取り Jensenの不等式を適用すると，対数尤度の 2 次下界

k

l=1



2

= −

−1
∑

σ(ξl,k )exp(−Al s21 + Bl s1 + Cl )
ql



Ki



+

∑

∑

2Ki −1

ql Al  s21 + ][
ql Cl − ln ql +

l=1

ql ∝

∏

k

σ(ai,k (za,k − bi,k )) である (図 2)．

各 k について式（17）のシグモイド関数の変分下界を適用する≥ −λ(ξi,k )(ai,k (za,k − bi,k ))2 +ql Bl  s1

+const∑≥−]ln σ(ξl,k )(19)2(λ(ξi,k )a2i,k )2(

1(ai,k (za,k − bi,k ))2(za,k −

bi,k +1
4λ(ξi,k )ai,k))2+const．

(22)

(23)

k

以上により，正解側尤度関数のガウス近似は N (za |η, Ψ) とが得られる．（図 4 右の下界 2）ただし，[正解側尤度のガウス近似

正解側の尤度関数は

ln σ (ai,k (za,k − bi,k ))

l=1

l=1
2Ki −1WT.

Λ\s1

について平方完成することで以下が得られる．

[∏
ql ln

]
0

ことで，指数 2 次形式の下界を得る．さらに対数を取り，za,kln p(y = 0|za )

∑

0

T

k

≥

.

Λ\s1

な尤度関数のガウス近似 N (za |η, Ψ) が得られる．（図 5 右の下界）

2Ki −1]0

これを再び変換式 (14) を利用して z 座標系に戻すと，最終的λ(ξl,k )ã2l,i,k

k

∑

./

∏

なる．共分散行列は対角行列である．

]
σ(ξl,k ) exp(−Al s21 + Bl s1 + Cl )(20)

ηk = bi,k +

k

3. 3.

とする．

1
,
4λ(ξi,k )ai,k2
Ψ−1
k,k = 2λ(ξi,k )ai,k .

変分パラメータの最適化

次に，ステップ 2 ではステップ 1 で求めた勾配方向での 1

3. 2. 1 節，3. 2. 2 節で導入した変分下界には変分パラメータ次元下界を多次元での下界にする．まず，s 座標系におけるが存在する．本節では，不正解側，正解側の順に変分パラメー

多次元下界を求め，それを z 座標系に変換することで最終的タの最適化方法を述べる．

に z 座標系における多次元ガウス分布を求める．事前分布を

まず不正解側の変分パラメータは ξ である．前節で求めた s

N (za |ma , Ga,a ) とすると，変換式 (14) を用いて s 座標系での座標系における s1 成分の 1 次元下界 lb1 (s1 , ξ) と，事前分布の表現を求めると，精度行列は Λ = W T G−1
a,a W となる．さらに，

s1 についての周辺確率 p(s1 ) = N (s1 |ms1 , σs21 ) を用いると，

シューア補行列の変換を利用すると，s1 方向を周辺化した事前∫
p(y = 0|za )p(s1 )ds1 ≥分布の精度行列は−1
Λ1,2:KiΛ\s1 = Λ2:Ki ,2:Ki − Λ2:Ki ,1 Λ1,1

∫

(21)

lb1 (s1 , ξ)p(s1 )ds1(24)が成り立つ．右辺を ξ に関して最大化することにより，変分パ—5—ラメータ ξ を最適化する．右辺はガウス分布の積の形をしてい+2ξl,k ãl,i,k λ(ξl,k ) −るため指数部を s1 に関して平方完成することで，∫λ′ (ξl,k ) = −lb1 (s1 , ξ)p(s1 )ds1

∑

(
∏

l=1

k

2Ki −1

=

√

1

∑

=

(

(
∏

l=1

σ(ξl,k )

k

(

+ Bl

1
2
σs

+ 2Al

2

正解側の変分パラメータの最適化は，

∫
∫

)


+ Cl 

1
1+ ai,k (za,k − bi,k ) − ξk + λ(ξk )ξk2
2
2

(25)

2πσk2

と計算できる．次に，対数を取り Jensen の不等式を適用する

−1
∑

[
∑

2

≥

Ql

l=1

(

k
m s1
2
σs


m2s
− 21
2σs1



+ Cl  + const(26)(:= lb2 (ξ))

Ql ∝σ(ξl,k )



m s1
2
σs

(

1 + 2Al σs21

)− 12

′

m2s
− 21
2σs1



+ Cl 

が得られる．ただし，

s1{とした．次に，lb2 (ξ) を ξ に関して微分すると，

=

∑

Ql

[

∑

l

−

k

(ms1 + Bl σs21 )2 ∂Al
∂Cl
−
+
∂ξ
(1 + 2Al σs21 )2 ∂ξ

∂Bl
=
∂ξ

4.

ãl,i,k λ′ (ξl,k )k2ã2l,i,k b̃l,i,k λ′ (ξl,k )

∑ ′
∂Cl
2
=
λ (ξl,k )ãl,i,k (ξl,k− b̃2l,i,k )∂ξk2λ(ξk )a2i,k bi,k +S2=2λ(ξk )a2i,k + ϕkai,k+ ma,k ϕk
2

実験

本章では，3. 2 節で説明した局所ガウス近似により α̂ メッ

]

∑ 3

k

=

索などを用いて，微分が 0 となる ξk を求めることができる．

となる．ただし，∑S1とする．以上で微分が求まったため，不正解側と同様に 2 分探∂ξl,k
(1 − σ(ξl,k ))∂ξ

σs21
ms1 + Bl σs21 ∂Bl∂Al+
2(1 + 2Al σs1 ) ∂ξ
(1 + 2Al σs21 ) ∂ξ

∂Al
=
∂ξ

2

λ (ξk )ai,k∂lb3
1= (1 − σ(ξk )) −− λ′ (ξk )a2i,k b2i,k −
∂ξk
S2
2
(
)
S1 ′
S1
+λ′ (ξk )ξk2 + 2λ(ξk )ξk +λ (ξk )a2i,k 2bi,k −
S2
S2



)2
+ Bl

)
exp  ( 12 σ12 + 2Al∂
lb2 (ξ)
∂ξ

k

すると，

k

(

)

となる．ただし，ϕk = σ12 とした．さらに ξk に関して微分が得られる．ただし，

)

(

1
1
ln λ(ξk )a2i,k + ϕk
2
2
1
−λ(ξk )a2i,k b2i,k − ξk + λ(ξk )ξk2
2
2λ(ξk )a2i,k bi,k + 12 ai,k + ma,k ϕk++ const4λ(ξk )a2i,k + 2ϕk

lb3 (ξk ) = ln σ(ξk ) −

s1

(
∏

(29)

ることで積分が計算でき，

1
ln (1 + 2Al σs21 )
2

)2

)
+ ( 1
2 σ12 + 2Al

1(za,k − ma,k )2 dza,k2σk2

辺はガウス分布の積の形をしているため，指数部を平方完成すln σξl,k −

+ Bl

)

exp −

の右辺を ξk について最大化することによって行う．ここで右

lb1 (s1 , ξ)p(s1 )ds1
Ki}(:= lb3 (ξk ))

ことで

∫

(

1

√

1

ln

{

σ(ξk )exp −λ(ξk )a2i,k (za,k − bi,k )2

≥ ln


m2s
− 21
2σs1

σ (ai,k (za,k − bi,k )) p(za,k )dza,k

ln

)

−1
1 + 2Al σs21 2

)2

m s1
2
σs


exp  ( 1次に，正解側の変分パラメータ ξk の最適化について述べる．

1(s1 − ms1 )2 ds1
2σs21

(

(28)

なる ξ を求めることができる．

))

λ(ξl,k )1+
σ(ξl,k ) (1 − σ(ξl,k ))
ξl,k
2ξl,k

である．微分が求まったため，2 分探索などを使い微分が 0 とexp(−Al s21 + Bl s1 + Cl )σ(ξl,k )exp

2πσs212Ki −1)∫1ãl,i,k2(27)

セージがどう近似されるかを人工データにより示し，次に AS-

SISTments Data を用いた予測精度の実験結果について示す．
4. 1.

α̂ メッセージの近似

非補償型の項目応答モデルにおける尤度関数を局所ガウス近似することで得られる近似事後分布 (α̂ メッセージ) と真の事後分布を人工データにより比較する．紙面の制約上，不正解側と正解側それぞれ 2 つの場合のみを示す．
まず，図 6 に不正解側の尤度関数のガウス近似についての結果を示す．図 6a と図 6b は，ある尤度関数（左図の赤）と事前分布

—6—

(a)

(a)

(b)

(b)

図 6: 不正解側のガウス近似と近似事後分布

図 7: 正解側のガウス近似と近似事後分布
表 1: 正解予測の結果（左図の緑）が与えられた時の，尤度関数の局所ガウス近似（左図の青）と真の事後分布（右図の赤）と近似事後分布（右図の青）を示手法している．尤度関数は p(y = 0|z) = 1 − σ(z1 − 3.0)σ(z2 − 5.0)

提案手法

0.763 ± 0.005 0.827 ± 0.004 0.794 ± 0.004

である．左図を見ると，尤度関数が事前分布付近で局所ガウス

DKT

0.743 ± 0.004 0.892 ± 0.006 0.811 ± 0.003

近似され，異なるガウス近似が得られていることが分かる．ま

適合率

再現率

F値

DKVMN 0.748 ± 0.004 0.838 ± 0.024 0.790 ± 0.010

た，事前分布の中心の勾配方向と直交する方向では事前分布と表 2: 不正解予測の結果

同じ分散を尤度関数のガウス近似が持っている．右図を見ると，

手法

近似事後分布が真の事後分布をおおよそ捉えていることが見て提案手法

0.609 ± 0.013 0.510 ± 0.010 0.555 ± 0.011

DKT

0.667 ± 0.015 0.413 ± 0.026 0.510 ± 0.024

取れる．中心はほぼ等しい一方で，分散は過小評価している．
これは，尤度関数を下からガウス分布で抑えるためこのような

適合率

再現率

F値

DKVMN 0.601 ± 0.017 0.461 ± 0.035 0.521 ± 0.020

傾向となることが考えられる．
次に図 7 に正解側の尤度関数のガウス近似の結果を示す．

習ベースの DKT [7]，DKVMN [8] を用いた．評価は 5 分割交図の見方は不正解側と同様である．尤度関数は p(y = 1|z) =差検定を用いて行い，正解予測と不正解予測について適合率，σ(z1 − 3.0)σ(2.0(z2 − 5.0)) である．左図を見ると，事前分布

再現率，F 値を算出した．また DKT では n_hidden パラメーに応じて尤度関数が局所ガウス近似されていることが分かり，

タを {50, 100, 200} の中で最も良いものを選択し，DKVMN

また正解側のガウス近似は各次元ごとに行うため，ガウス分布

は memory_size と (q_embed_dim, qa_embed_dim) についの軸は z 座標系と並行になる．右図を見ると近似事後分布が真て，{20, 50, 100}, {(50, 100), (100, 200)} のグリッドで最も良いの事後分布をほぼ正確に近似できていることが分かる．正解側

ものを選択した．結果を表 1 と表 2 にそれぞれ示す．正解予測の方が不正解側の近似よりも正確であることが分かる．

では，DKT がもっと良い F 値を示し，不正解予測では提案手

4. 2.

予測精度

法が最も良い F 値を示した．応用を考えると，不正解予測の性

2009-2010 ASSISTment Data (Skill Builder) [14] を使い予能が良い方が，学習者へのメリットが大きい．また，正解予測

測精度の実験を行った．ASSISTment Data は，生徒がｅラーと不正解予測の総合 F 値を図 8 に示す．総合的には提案手法がニングシステムを用いて問題を解いた時系列の正誤ログデー

DKT, DKVMN に対して高い F 値を示していることが分かる．

タであり，問題に対して複数のスキルタグが付与されている．

Skill Builder データは，スキルを育成することを目的としたも

5.

関連研究

ので，3 問連続問題を正解しないと宿題が完了できない．Skill

モデルが解釈可能なホワイトボックス型のナレッジトレーシ

Builder データのうち，少なくとも 10 名以上が解いている問ングに関する研究について述べる．まず，1 つのスキルを扱う場

題に関するログに絞り，ユーザ数 4,106，問題数 9,151，スキ合と複数のスキルを扱う場合に大きく分かれ，複数のスキルを

ル数 106 のデータセットを作成した．比較手法として，深層学扱う場合，さらに人間が定義したスキルタグを利用する方法とそ—7—<sDE
<d
ᥦ᱌ᡭἲ
Ϭ͘ϱϬϬ

Ϭ͘ϱϱϬ

Ϭ͘ϲϬϬ

Ϭ͘ϲϱϬ

Ϭ͘ϳϬϬ

図 8: 正解と不正解の総合 F 値
うでない方法に分かれる．1 つのスキルの習熟を追跡する方法で最も有名なものが Bayesian Knowledge Tracing(BKT) [15]である．BKT は隠れマルコフモデルであり，学習者の潜在ス
キル変数をバイナリの確率変数で表す．BKT は，[16] により，学習者ごとの個々の特性を反映する拡張が行われている．
複 数 の ス キ ル を 扱 う ナ レ ッ ジ ト レ ー シ ン グ で ，人 間 が定義したスキルタグを用いる方法には [17] [18] [19] [20] がある．[17] [18] [19] は，1 つの問題に対し 1 つスキルを対応付ける一方で，[20] は 1 つの問題に複数のスキルを対応付ける．[20]は，出力確率はバイナリではなくスコア値であるが各スキルの総和でスコアを回帰するため補償型の出力確率を持つモデルである．
複数のスキルを扱うナレッジトレーシングで，人間が定義したスキルタグを必要としない方法に [6] [21] がある．[6] では，線形動的システムをベースとしており，本研究との大きな違いは
出力確率が補償型か，非補償型かの違いである．補償型の出力確率であっても事後確率を解析的に求めることはできないため，

EP 法を用いて事後確率をガウス近似する．[21] では，ナレッジ
トレーシングをテンソル分解として定式化する．これらの方法
は，人間が定義した問題とスキルの対応表を必要としないため適用範囲が広いが予測理由の説明性は低い．

6.

ま と

め

本研究では，非補償型の項目応答モデルを持つ確率モデルによる解釈可能なナレッジトレーシングを提案した．提案モデル
は，線形動的システムをベースとして出力確率が非補償型の項
目応答モデルであった．非補償型の項目応答モデルを導入する
ことで学習者の潜在スキルの事後分布が解析的に求められなくなる問題に対し，非補償型の尤度関数を局所的にガウス近似することで，近似事後分布を解析的求める方法を導出し，生成モ
デルのパラメータは LDS と同様にして学習可能である．人工データにより，提案方法が真の事後分布を近似できていることを示し，さらに ASSISTment Data を使い深層学習ベースの DKT, DKVMN に対しても高い予測性能を示すことを確認
した．
文

献

[1] Knewton - Achievement Within Reach. https://www.
knewton.com/.
[2] DreamBox Learning - Online Math Learning. https://www.
dreambox.com/.
[3] Hiroaki Ogata, Chengjiu Yin, Misato Oi, Fumiya Okubo,
Atsushi Shimada, Kentaro Kojima, and Masanori Yamada.
E-book-based learning analytics in university education. In
ICCE, pages 401–406, 2015.

[4] Hao Cen, Kenneth Koedinger, and Brian Junker. Learning
factors analysis–a general method for cognitive model evaluation and improvement. In ITS, pages 164–175. Springer,
2006.
[5] Philip I. Pavlik, Hao Cen, and Kenneth R. Koedinger. Performance factors analysis –a new alternative to knowledge
tracing. In AIED, pages 531–538. IOS Press, 2009.
[6] Andrew S Lan, Christoph Studer, and Richard G Baraniuk.
Time-varying learning and content analytics via sparse factor analysis. In SIGKDD, pages 452–461. ACM, 2014.
[7] Chris Piech, Jonathan Bassen, Jonathan Huang, Surya
Ganguli, Mehran Sahami, Leonidas J Guibas, and Jascha
Sohl-Dickstein. Deep knowledge tracing. In NeurIPS, pages
505–513, 2015.
[8] Jiani Zhang, Xingjian Shi, Irwin King, and Dit-Yan Yeung.
Dynamic key-value memory networks for knowledge tracing.
In WWW, pages 765–774, 2017.
[9] Ben Babcock. Estimating a noncompensatory irt model using metropolis within gibbs sampling. Applied Psychological
Measurement, 35(4):317–329, 2011.
[10] Run ze Wu, Qi Liu, Yuping Liu, Enhong Chen, Yu Su, Zhigang Chen, and Guoping Hu. Cognitive modelling for predicting examinee performance. In IJCAI, 2015.
[11] Greg CG Wei and Martin A Tanner. A monte carlo implementation of the em algorithm and the poor man’s data
augmentation algorithms. Journal of the American statistical Association, 85(411):699–704, 1990.
[12] Christopher M. Bishop. Pattern Recognition and Machine
Learning (Information Science and Statistics). SpringerVerlag, Berlin, Heidelberg, 2006.
[13] Tommi S Jaakkola and Michael I Jordan. Bayesian parameter estimation via variational methods. Statistics and Computing, 10(1):25–37, 2000.
[14] Mingyu Feng, Neil Heﬀernan, and Kenneth Koedinger. Addressing the assessment challenge with an online system thattutors as it assesses. User Modeling and User-Adapted Interaction, 19(3):243–266, 2009.
[15] Albert T Corbett and John R Anderson. Knowledge tracing: Modeling the acquisition of procedural knowledge. Usermodeling and user-adapted interaction, 4(4):253–278, 1994.
[16] Michael V Yudelson, Kenneth R Koedinger, and Geoﬀrey J
Gordon. Individualized bayesian knowledge tracing models.
In AIED, pages 171–180. Springer, 2013.
[17] Kevin H Wilson, Yan Karklin, Bojian Han, and Chaitanya
Ekanadham. Back to the basics: Bayesian extensions of
irt outperform neural networks for proficiency estimation.
arXiv preprint arXiv:1604.02336, 2016.
[18] Chaitanya Ekanadham and Yan Karklin. T-skirt: Online
estimation of student proficiency in an adaptive learning
system. arXiv preprint arXiv:1702.04282, 2017.
[19] R.V. Lindsey, M. Khajah, and M.C. Mozer. Automatic discovery of cognitive skills to improve the prediction of student learning. NeurIPS, 2:1386–1394, 01 2014.
[20] Yuying Chen, Qi Liu, Zhenya Huang, Le Wu, Enhong Chen,
Runze Wu, Yu Su, and Guoping Hu. Tracking knowledge
proficiency of students with educational priors. In CIKM,
pages 989–998. ACM, 2017.
[21] Shaghayegh Sahebi, Yu-Ru Lin, and Peter Brusilovsky. Tensor factorization for student modeling and performance prediction in unstructured domain. EDM, 2016.

—8—