[论文翻译][2015][28]Bayesian Estimation of the DINA Model With Gibbs Sampling(基于Gibbs采样的DINA模型贝叶斯参数估计方法)

下载论文

摘要

基于已经提出的DINA模型的贝叶斯表达公式，可以使用Gibbs采样通过联合后延分布来拟合试题的猜对率、失误率、被试属性和潜在认知状态的概率。本程序扩展了Be guain等人的方法以估计三参数和四参数正则模型中的失误率和猜对率。模拟实验表明该模型具备较好的参数拟合能力。该技术已被应用与心理旋转测试，研究人员可以通过“DINA”R软件包免费获取该算法和插件。
关键词： 认知诊断; 贝叶斯统计; 蒙特卡洛；空间认知

引言

认知诊断模型(CDMs)是诊断学生技掌握状态的重要手段。研究者提出了大量的认知诊断模型，如DINA模型、NIDA模型、DINO模型、RUM模型、GDINA模型、LCDM模型和通用诊断模型等。上述模型的复杂程度不同，且先前的研究表明，诸如DINA模型、NIDA模型、DINO模型可以被应用在更为广泛的认知诊断框架中。DINA模型因其具备较好的可解释性成为当今比较流行的认知诊断模型。本篇文章的目的是运用Gibbs采样，并引入引入DINA模型的贝叶斯表达式来进行参数估计。

本文包括5个章节。第一章介绍DINA模型的基本概念并基于之前的研究讨论本文推荐的贝叶斯表达式的优势。第二章详细阐述贝叶斯表达式和如何基于Gibbs采样运用全条件分布来估计被试（技能掌握情况和属性信息）和试题（猜对率、失误率）的相关参数。第三章通过模拟数据研究，来展示模型的参数拟合效果。第四章展示模型在心里选择测试中的应用。最后一章包含讨论、未来的研究方向和总结回顾。

DINA模型

令 $i=1,\dots ,N$代表每个被试的编号，$j=1,\dots ,J$代表试题的编号，$K=1,\dots K$代表技能编号，$c=1,\dots ,C$代表潜在技能组合的编号。认知诊断模型假定知道每到试题覆盖了哪些技能点，则$J$道试题会产生一个$J\times K$结构化矩阵，矩阵中的元素可以为0或者1( q j , k ∈ { 0 , 1 } q_{j,k} \in \lbrace 0,1\rbrace qj,k​∈{0,1})代表答对试题$j$是否需要技能点$k$。也就是说$q_{j,k}=1$代表答对试题$j$需要技能点$k$,如果为0,则代表不需要技能点$k$。令$\mathbf{q}=(q_{j,1},\dots ,q_{j,K})$即Q矩阵中的一行，代表试题$j$所需要的技能向量。

DINA模型假定每位测试者在$K$个技能点上都有掌握和不掌握两种状态。如果被试$i$掌握技能$k$,则令${\alpha }_{i,k}=1$,否则令${\alpha }_{i,k}=0$。则被试$i$的技能状态可以用$K$维的二分向量来标识${\mathbf{\alpha }}_{\mathbf{i}}=({\alpha }_{i,1},\dots ,({\alpha }_{i,K})$，$K$个技能点，共可以产生$2^K$个技能点组合，${\mathbf{\alpha }}_{\mathbf{i}}$便是其中的一种。对于技能组合(分类)$c$来说，其对应的技能组合维${\mathbf{\alpha }}_{\mathbf{c}}=({\alpha }_{c,1},\dots ,({\alpha }_{c,K})$,若技能组合$c$包含技能点$k$，则令${\alpha }_{c,k}=1$，否则令其为0。

DINA模型是一款连接模型，这意味着，不考虑猜对的情况下，被试需要掌握向量${\mathbf{q}}_{\mathbf{j}}$中所有的知识点才能答对试题。我们令被试$i$的理想答题结果为${\mathbf{\eta }}_{\mathbf{i}}=({\eta }_{i,1},\dots ,{\eta }_{i,J})$,其中
$${\eta }_{i,j}=\prod _{k=1}^K{\alpha }_{i,k}^{q_{j,k}}=\mathcal{I}({{\mathbf{\alpha }}_{\mathbf{i}}}^T\cdot {\mathbf{q}}_{\mathbf{j}}={{\mathbf{q}}_{\mathbf{j}}}^T\cdot {\mathbf{q}}_{\mathbf{j}})$$

其中$\mathcal{I}$为指示函数。当${\eta }_{i,j}=1$表示被试$i$掌握了试题$j$所需要的所有技能点，当${\eta }_{i,j}=0$代表被试$i$未完全掌握答题试题$j$所需的技能点。在DINA模型中，被试$i$答对试题$j$的概率是一个关于被试技能组合${\mathbf{\alpha }}_{\mathbf{i}}$和试题参数的函数：
$$P(y_{i,j}=1|{\mathbf{\alpha }}_{\mathbf{i}},{\mathbf{\Omega }}_{\mathbf{j}})=\{\begin{array}{ll}{g}_j,& {\eta }_{i,j}={\prod }_{k=1}^K{\alpha }_{i,k}^{q_{j,k}}=0\\ 1-s_j,& {\eta }_{i,j}={\prod }_{k=1}^K{\alpha }_{i,k}^{q_{j,k}}=1\end{array}\text{\hspace{1em}(1)}$$
${\mathbf{\Omega }}_{\mathbf{j}}=(g_j,s_j)$是试题参数的行向量，其中$g_j$为未掌握试题$j$所需的知识点但是能猜对该题的概率$g_j=P(y_{i,j}=1|{\eta }_{i,j}=0)$； 相反，$s_j$代表掌握了试题$j$所需的知识点，但是失误答错的概率$s_j=P(y_{i,j}=0|{\eta }_{i,j}=1)$。

DINA模型的响应函数可以用于构建可观察的试题得分结果的似然函数。令${\mathbf{y}}_{\mathbf{i}}=(y_{i,1},y_{i,2},\dots ,y_{i,J})$，用该$J$维向量来表示第$i$个学生$J$题的答题结果，其中$y_{i,j}$为1代表第$i$个学生答对了$j$题，为0代表未大对该题。令$\mathbf{\Omega }=({\mathbf{\Omega }}_1,{\mathbf{\Omega }}_2,\dots ,{\mathbf{\Omega }}_{\mathbf{J}})$代表$J$题的猜对率和失误率。同时，令$P_{j,c}=P(y=1|{\mathbf{\alpha }}_{\mathbf{c}},{\mathbf{\Omega }}_{\mathbf{j}})$代表知识状态为${\mathbf{\alpha }}_{\mathbf{c}}$时，答对试题$j$的概率，注意${\alpha }_i={\alpha }_c$代表被试$i$为第$c$种(类)知识状态。假定试题的答题结果条件依赖给定的${\mathbf{\alpha }}_{\mathbf{i}}$，则在给定被试为第$c$中知识状态下，观察到答题结果$y_i$的概率为：
$$p(y_i|{\mathbf{\alpha }}_{\mathbf{i}}={\mathbf{\alpha }}_{\mathbf{c}},\mathbf{\Omega })=\prod _{j=1}^J{P}_{j,c}^{y_{i,j}}(1-P_{j,c}{)}^{(1-y_{i,j})}\text{\hspace{1em}(2)}$$

令${\pi }_c=P({\mathbf{\alpha }}_{\mathbf{c}})$代表第$c$种知识状态出现的概率, 令$\mathbf{\pi }=({\pi }_1,{\pi }_2,\dots ,{\pi }_C)$代表C种知识状态各自出现概率的向量。在给定被试的答题结果$y_i$，试题的参数$\mathbf{\Omega }$，各种知识状态出现的概率$\mathbf{\pi }$,则依据边际概率定义可知：
$$P(y_i|\mathbf{\Omega },\mathbf{\pi })=\sum _{c=1}^C{\pi }_{c}P(y_i|{\mathbf{\alpha }}_{\mathbf{i}}={\mathbf{\alpha }}_{\mathbf{c}},\mathbf{\Omega })\text{\hspace{1em}(3)}$$

则观察到的$N$个被试$J$题的答题结果发送的似然函数维：
$$p(\mathbf{y}|\mathbf{\Omega },\mathbf{\pi })=\prod _{i=1}^N\sum _{c=1}^C{\pi }_{c}P(y_i|{\mathbf{\alpha }}_{\mathbf{i}}={\mathbf{\alpha }}_{\mathbf{c}},\mathbf{\Omega })\text{\hspace{1em}(4)}$$

其中$\mathbf{y}=({\mathbf{y}}_1,{\mathbf{y}}_2,\dots ,{\mathbf{y}}_{\mathbf{N}})$是观察到的 $N\times J$ 被试答题结果矩阵。

贝叶斯公式主要关注参数的后验概率分布，对于 DINA模型来说：
$$p(\mathbf{\Omega },\mathbf{\pi }|\mathbf{y})\propto p(\mathbf{y}|\mathbf{\Omega },\mathbf{\pi })p(\mathbf{\Omega })p(\mathbf{\pi })\text{\hspace{1em}(5)}$$

其中$p(\mathbf{\Omega })$未试题参数的先验分布，$p(\mathbf{\pi })$为各种知识状态的先验概率分布。Maris¹指出了CMD参数后验分布难以处理的特性，并提出了使用蒙特卡洛采样曾向进行后验参数分布估计的方法。因此，众多后续研究采用了MCMC程序来对CDM参数进行后验分布估计²。Roussos等人对参数估计方法进行了总结(包括基于贝叶斯模型和非贝叶斯模型得参数估计)³。同时，之前得研究者还讨论了如何使用标准得统计软件对认知诊断模型进行参数估计。

贝叶斯参数估计优势

本文为DINA认知诊断模型开发得贝叶斯参数估计相对当前已有得程序具备如下优势。

第一、 模型公式可以得出学生和试题得可分析的完整条件分布，从而可以使用吉布斯才进行进行有效的计算。DINA模型和其他认知诊断模型之前采用的是基于马氏链的蒙特卡洛采样方法(MCMC-MH),这种方法需要调整$\mathbf{\Omega }$和$\mathbf{\pi }$的分布参数。因此每次使用MH算法对新数据进行分析时，需要重新对参数进行调整，具体来说，研究人员需要调整参数来确保击中率超过50%（拒绝率小于50%）。尽管诊断给定的数据仅需要数次调整参数即可，但是对于MH算法经验欠缺的人来说，会使得软件应用起来非常繁琐。相反，随后讨论的基于全条件分布的Gibbs采样，无需手工调整参数即可轻松使用。

再举一个例子，考虑后续模拟研究中，基于MH算法的额外需求。在编写算法代码并确定模拟条件后，需要确定$\mathbf{\Omega }$和$\mathbf{\pi }$的候选参数。后续参数需要反复实验进行调整，以确保接受率能达到25%到50%，因此需要花费大量的精力来进行参数调整。相反，使用本文开发的Gibbs采样程序只需提供先验参数即可。本文随后提出的程序，对$\mathbf{\Omega }$和$\mathbf{\pi }$引入了无信息的统一先验，因此可以快速使用算法进行参数估计。

第二、 本文提出的方案有利于进行Gibbs采样。比如，之前的CDM研究者们使用OpenBUGS⁴ Gibbs采样来对DINA模型进行参数估计，使用WinBUGS⁵ Gibbs采样来对HO-DINA模型进行参数估计。然而，Open BUGS和WinBUGS都存在一个缺陷，他们都没有明确规定DINA模型的单调性限制，即对所有试题而言 $0\le g_j+s_j\le 1$⁶。因为先前的Gibbs采样方法没有严格的执行可识别的条件，因此可能会导致应用程序的收敛性。随后讨论的模型扩展了三参数贝叶斯公式和四参数贝叶斯公式，来明确模型的单调性需求。因此本文的策略可以被扩展到其他包含失误率和猜对率的认知诊断模型中，如NIDA，NIDO和RUM。

第三、 与接受/拒绝采样方法相比(如MH算法)，Gibbs采样一个可以被证明的优势是改进了链的收敛。也就是说，先前的研究表明，在25%到50%的迭代之间接受后选择，则MH采样收敛效果最好。因此，假设后验采样的变化少于迭代次数的一半，MH采样可能导致后验采样的收敛速度变慢。相反，Gibbs采样有潜力利用完整的条件分布来进行采样来改善收敛效果。但是，正如一位匿名的评论家所言
，Gibbs采样不能保证所有的模型和参数收敛速度更快。随后给出的蒙特卡洛证据表明，MCMC算法在350到750次迭代后，参数的后验分布完成收敛，相反之前的使用Gibbs采样进行参数估计的方法，要进行1000到5000步迭代才能收敛。

第四、 如下文所述，学生潜在的属性一般是使用先验模型进行建模的。这么做的一个优点是，指定分类先验提供一种自然的方法来计算学生获取技能的可能性和学生知识状态潜在的分类。也就是说，可以总结每个学生技能分布的后延分布，来显示学生归属各个潜在类别的概率。

第五、 本文提供的算法已经以 ‘DINA’ R 软件包 的形式免费开放给研究者。本文提出的Gibbs采样使用C++语言编写，参数估计具有更高的效率。而且，R软件包中含有说明文档和插图，方便研究者上手和使用。

DINA R软件包地址：https://cran.r-project.org/web/packages/dina/index.html

基于贝叶斯的DINA模型参数估计方法

本章将介绍DINA模型的贝叶斯公式。第一部分将展现模型，并介绍先验概率。第二部分将为你呈现Gibbs采样的全条件分布实现细节。

贝叶斯模型表达式

依据DINA认知诊断模型的定义，可以得出如下贝叶斯公式：

$$y_{i,j}|{\mathbf{\alpha }}_{\mathbf{i}},{\mathbf{\Omega }}_{\mathbf{j}}\sim Bernoulli\left((1-s_j{)}^{{\eta }_{i,j}}{g}_j^{(1-{\eta }_{i,j})}\right),{\eta }_{i,j}=\mathcal{I}({{\mathbf{\alpha }}_{\mathbf{i}}}^T\cdot {\mathbf{q}}_{\mathbf{j}}={{\mathbf{q}}_{\mathbf{j}}}^T\cdot {\mathbf{q}}_{\mathbf{j}})\text{\hspace{1em}(6)}$$$$p({\mathbf{\alpha }}_{\mathbf{i}}|\mathbf{\pi })\propto \prod _{c=1}^C{\pi }_c^{\mathcal{I}({\alpha }_i={\alpha }_c)},0\le {\pi }_c\le 1,\sum _{c=1}^C{\pi }_c=1\text{\hspace{1em}(7)}$$$$p(\mathbf{\pi })=p({\pi }_1,{\pi }_2,\dots ,{\pi }_C)\propto Dirichlet({\mathbf{\delta }}_0),{\mathbf{\delta }}_0=({\delta }_{0,1},\dots ,{\delta }_{0,C})\text{\hspace{1em}(8)}$$ p ( Ω j ) = p ( s j , g j ) ∝ s j ( α s − 1 ) ( 1 − s j ) ( β s − 1 ) g j ( α g − 1 ) ( 1 − g j ) ( β g − 1 ) I ( ( s j , g j ) ∈ P ) P = { ( s , g ) : 0 ≤ s + g ≤ 1 , 0 ≤ s ≤ 1 , 0 ≤ g ≤ 1 } (9) p(\boldsymbol{\Omega_j}) = p(s_j, g_j) \propto s_j^{(\alpha_s - 1)}(1 - s_j)^{(\beta_s - 1)}g_j^{(\alpha_g - 1)}(1 - g_j)^{(\beta_g - 1)}\mathcal{I} \left((s_j,g_j) \in \mathcal{P} \right) \\ \mathcal{P} = \{(s,g): 0 \le s + g \le 1 , 0 \le s \le 1, 0 \le g \le 1\} \tag{9} p(Ωj​)=p(sj​,gj​)∝sj(αs​−1)​(1−sj​)(βs​−1)gj(αg​−1)​(1−gj​)(βg​−1)I((sj​,gj​)∈P)P={(s,g):0≤s+g≤1,0≤s≤1,0≤g≤1}(9)式6假定在给定学生知识状态${\mathbf{\alpha }}_{\mathbf{i}}$和试题参数${\mathbf{\Omega }}_{\mathbf{j}}$的情况下，试题的响应函数服从伯努利分布。事实上，式6只是对等式1的重新描述。式7假定被试知识向量${\mathbf{\alpha }}_{\mathbf{i}}$的分类先验，即属于类别$c$的概率： ${\pi }_c=p({\mathbf{\alpha }}_{\mathbf{i}}={\mathbf{\alpha }}_{\mathbf{c}}|\mathbf{\pi })$.

式8显示$\mathbf{\pi }$的共轭先验是以${\mathbf{\delta }}_0=({\delta }_{0,1},\dots ,{\delta }_{0,C})$为参数的Dirichlet分布。回想一下Dirchlet分布，它是Beta分布的多元版本，支持单位间隔。${\mathbf{\delta }}_0$中的元素值指明了各个分类的常见先验信息。如果没有$\mathbf{\pi }$的先验信息，可以将${\mathbf{\delta }}_0$设置为1，研究者也可以根据世界情况调整${\mathbf{\delta }}_0$的值。比如${\pi }_c$的边际分布服从$Beta({\delta }_{0,c},{\sum }_{c^{\prime }\ne c}{\delta }_{0,c^{\prime }})$，词${\mathbf{\delta }}_0$的值可以选用能体现观察到的${\pi }_c$的值。 提示一下，下个章节讨论的蒙特卡洛方法的结果表明，对于任意的$c$,当${\delta }_{0,c}=1$时，模型具有很强的鲁棒性。

式9包含${\mathbf{\Omega }}_{\mathbf{j}}$的联合先验。Junker等人构建的DINA模型⁶需要$0\le s_j+g_j\le 1$ ,以此来提供潜在分类的区分信息。式9为${\mathbf{\Omega }}_{\mathbf{j}}$的先验在试题参数上明确施加了单调性限制。特别的，$p({\mathbf{\Omega }}_{\mathbf{j}})$的先验通过以${\alpha }_s、{\beta }_s、{\alpha }_g、{\beta }_g$为参数的$s_j$和$g_j$的Beta密度函数的积进行了重构，并通过$\mathcal{P}$来对参数$s_j$和$g_j$的取值范围上增加线性截断限制。通过选择${\alpha }_s、{\beta }_s、{\alpha }_g、{\beta }_g$可以将失误率和猜对率的先验合并到不集中。请注意,选择无信息先验${\alpha }_s={\beta }_s={\alpha }_g={\beta }_g=1$与对$s_j$和$g_j$的先验采用线性截断是一致的。随后讨论的模拟实验将正式使用对失误率和猜对率使用无信息先验是有效的。

全条件分布

对于${\mathbf{\alpha }}_{\mathbf{i}},\mathbf{\pi },{\mathbf{\Omega }}_{\mathbf{j}})$的全条件分布公式如下：
$$p({\mathbf{\alpha }}_{\mathbf{i}}|{\mathbf{y}}_{\mathbf{i}},\mathbf{\pi },\mathbf{\Omega })=\prod _{c=1}^C{\tilde{\pi }}_{i,c}^{\mathcal{I}({\mathbf{\alpha }}_{\mathbf{i}}^{\mathbf{T}}\cdot {\mathbf{q}}_{\mathbf{j}}={\mathbf{q}}_{\mathbf{j}}^{\mathbf{T}}\cdot {\mathbf{q}}_{\mathbf{j}})},{\tilde{\pi }}_{i,c}=p({\mathbf{\alpha }}_{\mathbf{i}}={\mathbf{\alpha }}_{\mathbf{c}}|{\mathbf{y}}_{\mathbf{i}},\mathbf{\pi },\mathbf{\Omega })\text{\hspace{1em}(10)}$$$$\begin{gathered} \mathbf{\pi }|{\mathbf{\alpha }}_{\mathbf{i}},\dots ,{\mathbf{\alpha }}_{\mathbf{N}}\propto Dirichlet(\tilde{\mathbf{N}}+{\mathbf{\delta }}_0) \\ \tilde{\mathbf{N}}=({\tilde{N}}_1,\dots ,{\tilde{N}}_C)=(\sum _{i=1}^N\mathcal{I}({\mathbf{\alpha }}_{\mathbf{i}}={\mathbf{\alpha }}_1),\dots ,\sum _{i=1}^N\mathcal{I}({\mathbf{\alpha }}_{\mathbf{i}}={\mathbf{\alpha }}_{\mathbf{C}}))\text{\hspace{1em}(11)} \end{gathered}$$$$\begin{gathered} p({\mathbf{\Omega }}_{\mathbf{j}}|{\mathbf{y}}_{\mathbf{j}},{\eta }_{1,j},\dots ,{\eta }_{N,j})\propto s_j^{({\tilde{\alpha }}_s-1)}(1-s_j{)}^{({\tilde{\beta }}_s-1)}{g}_j^{({\tilde{\alpha }}_g-1)}(1-g_j{)}^{({\tilde{\beta }}_g-1)}\mathcal{I}\left((s_j,g_j)\in \mathcal{P}\right) \\ {\tilde{\alpha }}_s={\mathcal{S}}_j+{\alpha }_s,{\tilde{\beta }}_s={\mathcal{T}}_j-{\mathcal{S}}_j+{\beta }_s \\ {\tilde{\alpha }}_g={\mathcal{G}}_j+{\alpha }_g,{\tilde{\beta }}_g=N-{\mathcal{T}}_j-{\mathcal{G}}_j+{\beta }_g \\ {\mathcal{T}}_j=\sum _{i=1}^N{\eta }_{i,j},{\mathcal{S}}_j=\sum _{i|y_{i,j}=0}^N{\eta }_{i,j},{\mathcal{G}}_j=\sum _{i|y_{i,j}=1}^N(1-{\eta }_{i,j})\text{\hspace{1em}(12)} \end{gathered}$$ 式10显示${\mathbf{\alpha }}_{\mathbf{i}}$的全条件分布式一个多分类分布，对于其中的分类$c$的概率用${\tilde{\pi }}_{i,c}$来表示。特别的，贝叶斯理论表明，给定答题结果${\mathbf{y}}_{\mathbf{i}}$, 试题参数$\mathbf{\Omega }$, 各个类别出现的概率$\mathbf{\pi }$的情况下学生$i$的认知状态属于类别$c$的概率如下：
$${\tilde{\pi }}_{i,c}=p({\mathbf{\alpha }}_{\mathbf{i}}={\mathbf{\alpha }}_{\mathbf{c}}|{\mathbf{y}}_{\mathbf{i}},\mathbf{\pi },\mathbf{\Omega })=\frac{p({\mathbf{y}}_{\mathbf{i}}|{\mathbf{\alpha }}_{\mathbf{i}}={\mathbf{\alpha }}_{\mathbf{c}},\mathbf{\Omega }){\pi }_c}{{\sum }_{c=1}^C{\pi }_{c}p({\mathbf{y}}_{\mathbf{i}}|{\mathbf{\alpha }}_{\mathbf{i}}={\mathbf{\alpha }}_{\mathbf{c}},\mathbf{\Omega })}\text{\hspace{1em}(13)}$$ ${\mathbf{\alpha }}_{\mathbf{i}}$的全条件分布是一个以${\tilde{\pi }}_{i,1},\dots ,{\tilde{\pi }}_{i,C}$为概率的多项式分布。

式11阐述了$\mathbf{\pi }$的全条件分布。特别的，$\mathbf{\pi }$的全条件分布可以被转化成如下形式：
$$\begin{gathered} p(\mathbf{\pi }|{\mathbf{\alpha }}_1,\dots ,{\mathbf{\alpha }}_{\mathbf{N}})\propto p({\mathbf{\alpha }}_1,\dots ,{\mathbf{\alpha }}_{\mathbf{N}}|\mathbf{\pi })p(\mathbf{\pi })=\left[\prod _{i=1}^{N}p({\mathbf{\alpha }}_{\mathbf{i}}|\mathbf{\pi })\right]p(\mathbf{\pi }) \\ =\left[\prod _{i=1}^N\prod _{c=1}^C{\tilde{\pi }}_c^{\mathcal{I}({\mathbf{\alpha }}_{\mathbf{i}}^{\mathbf{T}}\cdot {\mathbf{q}}_{\mathbf{j}}={\mathbf{q}}_{\mathbf{j}}^{\mathbf{T}}\cdot {\mathbf{q}}_{\mathbf{j}})}\right]\prod _{c=1}^C{\pi }_c^{{\delta }_{0,c}}=\prod _{c=1}^C{\pi }_c^{{\tilde{N}}_c+{\delta }_{0,c}}\text{\hspace{1em}(14)} \end{gathered}$$ 式14显示，$\mathbf{\pi }$的全条件分布是${\mathbf{\alpha }}_{\mathbf{i}}$多项式先验的乘积和$\mathbf{\pi }$的Dirichlet的先验分布。在MCMC算法迭代过程中，用向量$\tilde{\mathbf{N}}=({\tilde{N}}_1,\dots ,{\tilde{N}}_C)$来表示$C$种类别的学生人数，其中${\tilde{N}}_c={\sum }_{i=1}^N\mathcal{I}({\mathbf{\alpha }}_{\mathbf{i}}={\mathbf{\alpha }}_{\mathbf{c}})$代表认知状态为$c$类的学生人数。相应的，式14中潜在类别概率的全条件分布服从$Dirichlet(\tilde{\mathbf{N}}+{\mathbf{\delta }}_0)$。

式12展示了试题参数$s_j$、$g_j$的全条件分布被简化成二项Beta分布。更精确的说，$s_j$和$g_j$的全条件分布可以转化成：
$$\begin{gathered} p(s_j,g_j|y_{1,j},\dots ,y_{N,j},{\mathbf{\alpha }}_1,\dots ,{\mathbf{\alpha }}_{\mathbf{N}})\propto \left[\prod _{i-1}^{N}p(y_{i,j}|{\mathbf{\alpha }}_{\mathbf{i}},s_j,g_j)\right]p(s_j,g_j) \\ \propto \left\{\prod _{i=1}^N{\left[(1-s_j{)}^{{\eta }_{i,j}}{g}_j^{(1-{\eta }_{i,j})}\right]}^{y_{i,j}}{\left[s_j^{{\eta }_{i,j}}(1-g_j{)}^{1-{\eta }_{i,j}}\right]}^{(1-y_{i,j})}\right\}\cdot \\ {s}_j^{({\alpha }_s-1)}(1-s_j{)}^{({\beta }_s-1)}{g}_j^{({\alpha }_g-1)}(1-g_j{)}^{({\beta }_g-1)}\mathcal{I}\left((s_j,g_j)\in \mathcal{P}\right) \\ =s_j^{({\tilde{\alpha }}_s-1)}(1-s_j{)}^{({\tilde{\beta }}_s-1)}{g}_j^{({\tilde{\alpha }}_g-1)}(1-g_j{)}^{({\tilde{\beta }}_g-1)}\mathcal{I}\left((s_j,g_j)\in \mathcal{P}\right)\text{\hspace{1em}(15)} \end{gathered}$$式15的最终结果通过简化个体的乘积和采集$s_j,1-s_j,g_j,1-g_j$的指数项进行构建。定义完式12后，可以回顾一下${\tilde{\alpha }}_s,{\tilde{\beta }}_s,{\tilde{\alpha }}_s{\tilde{\alpha }}_s$。在MCMC迭代过程中，令${\mathcal{T}}_j={\sum }_{i=1}^N{\eta }_{i,j}$代表掌握答对试题$j$所需知识点的学生人数。${\mathcal{S}}_j={\sum }_{i|y_{i,j}=0}^N{\eta }_{i,j}$代表掌握试题$j$所需知识点，但是没有答对的学生人数,$g_j={\sum }_{i|y_{i,j}=1}^N(1-{\eta }_{i,j})$代表没有掌握试题$j$所需知识点，但是猜对的学生人数。$s_j,g_j$的全条件分布简化成了以 ${\tilde{\alpha }}_s,{\tilde{\beta }}_s,{\tilde{\alpha }}_s{\tilde{\alpha }}_s$为参数的Beta分布。附录中展示了在$0\le s_j+g_j\le 1$单调性约束条件下，$s_j$和$g_j$的联合采样方法。

蒙特卡洛仿真研究

本节将阐述蒙特卡洛仿真的结果以证实MCMC算法的准确性，并提供有关收敛的其他信息。仿真研究类似于Torre等人设计的方案。仿真研究中，$K=5,J=30,N=1000$,并采用Torre等人之前使用的$Q$矩阵。仿真研究关注三个话题：

第一、研究猜对率和失误率的两个条件，以此来理解潜在技能和试题如何影响参数的恢复(估计)。高诊断场景(潜在技能与试题之间具有强关系，猜对率和失误率都很低)将使用$g_j=s_j=0.1$,低诊断场景将使用$g_j=s_j=0.2$

第二、需要考虑两种场景下潜在类别的概率分布。回顾一下，$\mathbf{\pi }$表示学生的技能分布，${\pi }_c$代表认知状态为潜在类别${\mathbf{\alpha }}_{\mathbf{c}}$的学生占比。在第一种场景下，潜在类别的概率使用常数${\pi }_c=2^{-K}$(被试的技能分布是平均的，即各种认知状态发生的概率相同)。在第二种场景下，包含3个或以上知识点的认知状态类别发生的可能性是小于3个知识点的认知状态类别的两倍。也就是说，如果${\mathbf{\alpha }}_{\mathbf{c}}$中包含的知识点小于3个，则该知识状态发生的概率为${\pi }_c=0.021$, 如果包含的知识点大于等于3个，则该知识状态发生的概率为${\pi }_c=0.042$。

第三、仿真研究需对多项式后验分布收敛的迭代次数进行估计。对于每种场景，运行了5个链，并使用’coda’ R程序包计算了多项式潜在规模的缩小系数$\hat{R}$，来评估$\mathbf{\Omega }$和$\mathbf{\pi }$的收敛情况。$\hat{R}$的值小于1.1或1.2时就证明了该条例已经收敛了。注意$\hat{R}$是针对50到2000的链长计算得来的，需要通过burn-in(去除前面数个不稳定的样本点)来达到收敛状态。

MCMC算法实现需要5000次迭代，其中前面1000次迭代会被丢掉(burn-in)。通过重复计算50次的$s_j,g_j,\pi$的平均值和均方根误差(RMSE)来评估准确性。Beta分布和Dirichlet分布使用无信息先验${\alpha }_s={\beta }_s={\alpha }_g={\beta }_g=1,{\mathbf{\delta }}_0=1_{32}$,其中$1_{32}$代表这一个32维全是1的向量。

图1提供了仿真后验分布收敛的正确。特别的，图1包含8个绘图面板，绘制了50个副本的$\hat{R}$。图1整理了诊断条件会影响收敛，但是总体技能分布的本质不会影响$\hat{R}$。比如，在高诊断场景中($s_j=g_j=0.1$),经过300次迭代$\hat{R}<1.1$，相反在第诊断场景中($s_j=g_j=0.2$)大约经过750次迭代$\hat{R}<1.1$。因此，蒙特卡洛证据表明，对于仿真研究，采用1000的预烧应该足以满足后验分布的收敛。
$$图1：蒙特卡洛仿真条件下多元潜在规模减小因子\hat{R}(\hat{R}<1.1代表已经收敛)$$

图2和图3绘制了$s_j、g_j、\mathbf{\pi }$蒙特卡洛估计的期望值，以及失误率和猜对率的根均方误差RMSE。图2表明算法在8中仿真条件下取得了较小的偏差。图3绘制了8中仿真条件下，失误率和猜对率的根均方误差。图2表明算法可以对试题参数进行无偏估计，这意味着图3中的根均方误差可以量化抽样变异性。
$$图2：试题失误率和猜对率估计结果$$
如之前的研究类似，图3表明失误率参数的抽样变异性与试题需要的技能数量相关(试题所需知识点越多，失误率越大？)。有些包含三个知识点的试题，其失误率的抽样变异性是包含一个知识点试题的两倍。而且，失误率抽样变异性在平坦的$\mathbf{\pi }$分布场景（相对于高技能分布场景）更高。因此图3表明试题${\mathbf{q}}_{\mathbf{j}}$中的1的数量和分布影响着试题失误率的抽样变异性。相反，图3表明，需要更多技能的试题，其猜对率$g_j$的抽样变异性更小。
$$图3：失误率和猜对率的根均方误差$$
包含三个知识点的试题的失误率根均方误差RMSE超过了0.04，在偏差接近于0请情况下，误差范围约为$\pm 0.08$。存在两种可能的解释，一种是算法虽然是无偏的，单位效果很差；另外一种是图3观测到采样变异性可能是DINA模型的特征。请注意，Torre在项目所需的技能和项目参数的RMSE之间发现了相似的关系，但是失误率和猜对率的最大标准误差为0.031和0.017（N=2000，$\pi$为平坦分布）。Torre发现猜对率的根均方误差RMSE与$\mathbf{q}\prime \mathbf{q}$呈现负相关，猜对率的根均方误差RMSE与$\mathbf{q}\prime \mathbf{q}$呈现正相关。另一项模仿Torre的仿真研究使用了100个样本也证明了这两种解释的合理性。在后续的仿真研究中，失误率和猜对率最大的RMSE为0.032和0.017，这为后一种解释提供了证据，即失误率(猜对率)对于包含更多(更少)技能的试题具有更大的变异性。

图4和图5给出了使用蒙特卡洛方法，潜在技能分布$\mathbf{\pi }$的偏差和RMSE。与图2类似，算法在整个防止条件得出了$\mathbf{\pi }$的无偏估计。
$$图4:潜在类别发生概率\pi 的估计结果$$
图5展示了$\mathbf{\pi }$的RMSE，并展示了其采样变异性。与图3不同的是，图5表明$\mathbf{\pi }$的RMSE与诊断条件以及潜在类别的技能分布无关。然而在图5的$h$面板中,${\pi }_c$在技能较少的潜在类别中的采样变异性却增加了。
$$图5:潜在类别发生概率\pi 的RMSE$$

简而言之，蒙特卡洛防止研究可以得出以下几个结论：首先，新的Gibbs采样方法可以准确估计模型参数，其次，当采用非先验信息时，算法时准确的，第三，结果表明模型可以在合理的迭代次数后收敛。

应用：空时可视化测试

作者使用Gibbs采样对空间可视化测试（在三维空间中通过旋转和平移从一种状态转化成另外一种状态的测试）的进行参数估计(估计的结果中有的试题猜对率和失误率竟然超过了5,0%　(゜ロ゜))。该部分不进行翻译。

讨论

本研究为广泛使用的DINA模型提供了贝叶斯表达式。并开发了一款高效的贝叶斯采样算法来进行模型的参数估计。本节将在现有的研究背景下讨论本文的研究优点，同时对未来改进方向进行展望。

首先，本文提出的公式可以进行扩展来进行$\mathbf{Q}$的估计。认知诊断模型假定研究者能精确的试题中包含的知识点进行标记。最近的研究者们开始对$\mathbf{Q}$的估计程序及$\mathbf{Q}$矩阵的错误检查进行研究。比如，Chiu提供了优化$\mathbf{Q}$的概述方法，此外，Decarlo和Templin等人不是采用固有的方式来确定$\mathbf{Q}$矩阵的值，而是使用贝叶斯方法来随机选择$\mathbf{Q}$中元素的值。注意Templin等人才用了MH采样方法，而Decarlo使用了Gibbs采样方法：OpenBUGS 。本文提出的Gibbs采样方法可以结合先前的策略来对$\mathbf{Q}$矩阵进行建模，这么做的好处有：1、优化混合$、pmb\Omega$和$\mathbf{\pi }$ 2、强制加入单调性需求，3、无需微调方便使用。

第二，DINA模型时众多通用认知诊断模型的一个特例。之后的研究可以在更多的认知诊断模型上使用贝叶斯公式。

第三，蒙特卡洛仿真研究表明，算法在估计参数的后验分布时，可以在合理的迭代次数下收敛。比如图1表明Gibbs采样可以在最多经过750次迭代后就可以收敛。而且实验表明，算法在使用$\mathbf{\Omega }$和$\mathbf{\pi }$无信息先验时，可以取得较高的精度。

计算$2^K$的潜在类别全条件分布并非没有缺点,“当$\mathbf{\alpha }$的维度超过$K=3或K=4$时，需要进行简化(Torre 2004)²”。MCMC程序基于Rcpp包使用C++语言开发，当$J=30,K=5,N=1000$时，在2.4GHz处理器、6GB内存的平台上，完成5000次迭代需要128秒。然而，Gibbs采样可能还需要更多的计算资源。例如Torre运用非线性因素分析对潜在技能背后的高阶因素进行建模，提出了HO-DINA模型，该模型中各个技能是相互独立的，通过给定的$K^{\prime }$($K^{\prime }\le K$)维连续值${\theta }_i$来进行计算，比如$p({\mathbf{\alpha }}_{\mathbf{i}}|{\mathbf{\theta }}_{\mathbf{i}})={\prod }_{k=1}^{K}p({\alpha }_{i,j}|{\mathbf{\theta }}_{\mathbf{i}})$。HO-DINA使用$K^{\prime }(K^{\prime }-1)/2$个参数来近似模拟$2^K$个潜在类别。如果$K=20$，会产生1048576个类别,本文提出的算法，在每轮迭代都会逐一遍历各个潜在类别。相反，在HO-DINA模型中，如果$K^{\prime }=20$,只需估计230个参数即可，如果近似适当，则计算上有非常大的优势。因此高维认知诊断需要跟进一步进行研究。

总的来说，本文为DINA模型提供了贝叶斯公式，为轻松估计模型参数提供了一种选择。

参考文献

[1]. Maris, E. (1999). Estimating multiple classification latent class models. Psychometrika,
64, 187–212.
作者研究了潜在多分类认知诊断模型得固有性质，首次提出使用蒙特卡洛采样(MCMC)方法对认知诊断模型进行参数估计
[2]. de la Torre, J., & Douglas, J. A. (2004). Higher-order latent trait models for cognitive diagnosis. Psychometrika, 69, 333–353.
作者提出了高阶DINA模型：HO-DINA，认为学生得知识点掌握情况受学生学习能力等高阶因素得影响，同时作者使用了MCMC方法对HO-DINA模型进行了参数估计
[3]. Roussos, L. A., DiBello, L. V., Stout, W., Hartz, S. M., Henson, R. A., & Templin, J. L. (2007). The fusion model skills diagnosis system. In J. P. Leighton & M. J. Gierl (Eds.), Cognitive diagnostic assessment for education: Theory and applications (pp.275–318). Cambridge, England: Cambridge University Press.
作者回顾了融合认知诊断模型，并对认知诊断模型得参数估计方法进行了总结该篇论文暂时不提供下载
[4]. DeCarlo, L. T. (2012). Recognizing uncertainty in the Q-matrix via a Bayesian extension of the DINA model. Applied Psychological Measurement, 36, 447–468.
作者使用贝叶斯方法对DINA模型中的Q矩阵进行合理性识别并对分数减法数据集的Q矩阵进行了矫正，该篇论文暂时不提供下载
[5]. Li F. A modified higher-order DINA model for detecting differential item functioning and differential attribute functioning |D|. uga, 2008.
文章基于HO-DINA提出了一个新的模型，并使用Gibbs采样对模型进行了参数估计。实验部分，作者于MH采样方法进行了对比。
[6]. Junker, B. W., & Sijtsma, K. (2001). Cognitive assessment models with few assumptions, and connections with nonparametric item response theory. Applied Psychological Measurement, 25, 258–272.
作者对认知诊断模型的单调性进行着论证，提出了一些假设（猜对率和失误率之和小于1），是的模型具有更好的可用性和可解释性。