★ FEATURED ARTICLE 做认知诊断项目十有八九是从一张对不齐的表格开始崩的。我印象最深的一次是把一份数学测评数据丢进模型跑完之后结果里有将近一半的学生被判定成“所有属性都没掌握”。查了两三天才发现Q矩阵的行顺序和作答数据的题目顺序错位了模型本身没问题是数据喂错了。DINA、HO-DINA、P-DINA、G-DINA这几个名字听起来像四套彼此独立的方法实际上是一条线上的四个台阶DINA是地基负责用slip和guess两个参数把“会”和“不会”拆开HO-DINA在属性之间引入一个高阶连续特质解决属性强相关时的估计不稳P-DINA处理多级计分题让选择题之外的量表题、步骤分题也能进来G-DINA则是一个饱和框架前三个都可以看作它在特定约束下的特例。这套东西主要的落地场景是教育测评、自适应练习系统的知识点定位、题库的属性标签体系、职业能力认证的分项报告以及各类心理量表的分维度解读。下面按“Q矩阵→DINA→HO-DINA→P-DINA→G-DINA→工程落地→踩坑排查”的顺序把数学定义、参数含义、代码实现和我实际踩过的坑一次讲透。1. Q矩阵先钉死再谈模型1.1 Q矩阵是什么为什么它决定模型上限Q矩阵是认知诊断里唯一一个完全靠人拍板的东西也是唯一一个模型学不出来的东西。它是一个J × K的0/1矩阵J是题目数K是属性数也就是你想诊断的知识点或技能维度q_jk 1表示第j题考了第k个属性0表示没考。比如一道“解一元二次方程”的题可能同时标了“因式分解”和“求根公式”两个属性这一行就是两个1。它的作用是把抽象的属性映射到具体的作答行为上。DINA模型的核心连接函数是η_ij Π_k α_ik^{q_jk}意思是学生i只有在第j题涉及的所有属性上都掌握了α_ik1对所有q_jk1的k成立理想状态下才会答对这道题此时η_ij 1只要缺一个η_ij 0。这就是所谓的“非补偿性”——属性之间不能互相顶替。会因式分解但不会求根公式这道题在模型眼里就是不会。这个设计的直接后果是Q矩阵错了模型学得再好也没意义。它会把“其实考了属性3”的题当成只考属性1然后强行给一批学生打上错误的标签。我一般跟团队说Q矩阵是模型的天花板你的Q标错10%后面的诊断报告就有10%的噪声而且这个噪声是有系统性的不是随机误差。另一个容易被忽略的点是Q矩阵决定了属性的可识别性。如果某个属性从来没有单独出现过每道涉及它的题都同时考了别的属性那这个属性的掌握状态在数据里就没有独立信息估计出来的结果基本靠先验撑着。1.2 Q矩阵的三种来源与验收清单实际项目里Q矩阵一般来自三个地方可靠性依次下降第一种是领域专家人工标注。教研老师对着每道题逐条勾选考了哪些知识点这是最靠谱的但成本高、一致性难保证。两个人标同一道题经常一个标2个属性一个标3个。我的做法是先让两个人独立标算一下Cohen‘s kappa低于0.6的题目拉出来开会定别硬上。第二种是从题库已有的标签体系映射。很多题库本身就有知识点标签但那些标签是给检索和组卷用的颗粒度往往和诊断需求不匹配。比如标签体系里“分数运算”是一个标签但诊断上你可能需要拆成“同分母加减”“异分母通分”“约分”三个属性。直接映射的结果通常是属性太粗诊断报告没法用。第三种是数据驱动反推用残差分析或者Q矩阵修正算法从作答数据里估计哪些题考了哪些属性。这条路我只在已有Q矩阵需要局部修正时用不会从零开始全靠它。原因是数据驱动方法对样本量和题目质量很敏感样本小的时候它会给你一堆看似合理但实际上没法学解释的标注。不管哪种来源上线前我都跑一遍这张检查表检查项判定标准不满足的后果空行检查每行至少1个1该题对诊断零贡献纯浪费属性覆盖每列至少3个1属性不可识别估计结果不可信单属性题每列最好有1道只考该属性的题属性间信息纠缠解释困难重复行不存在两行完全相同允许但题目区分度会下降重复列不存在两列完全相同两个属性完全共线等于一个属性矩阵秩非零行的秩 ≥ K属性维度不足模型欠定包含关系属性k的题目集合不被属性k’完全覆盖两属性难以分离“每列至少3个1”这条是经验值不是定理但在K4到6的常见规模下低于3基本估不准。“单属性题”这条很多人不当回事实际上有它的话属性掌握概率的估计方差会明显小一截因为单属性题的作答几乎直接反映了那个属性的状态。1.3 属性数与样本量动手前必须算清楚的一笔账DINA的E步要在所有2^K种属性掌握模式上算后验概率这个数字随K指数增长。K4是16种K6是64种K8是256种K12就是4096种。计算量本身还能忍真正的问题在样本模式越多每种模式分到的样本越少参数估计越飘。下面这张表是我按“每种模式平均至少要有15到25个被试”倒推出来的经验值只作为项目立项时的参考不是统计定理属性数K掌握模式数 2^K建议最小样本量实际项目中的常见规模38200300–500416400600–12005327001000–200066412002000–3500712820003000–5000825635005000以上且要谨慎样本不够的时候有三个应对办法我按推荐顺序排一是砍属性数把相关性极高的两个属性合并二是上HO-DINA用高阶结构把属性间的相关约束起来等效减少自由参数三是加先验用贝叶斯估计把slip和guess往合理区间拉。第三种最省事但也最“作弊”因为先验太强的时候数据说什么已经不重要了。2. DINA两个参数撑起的诊断逻辑2.1 slip与guess的物理含义DINA对每道题只估两个参数。s_j叫失误率slip意思是学生把这道题涉及的属性全掌握了却还是答错了的概率g_j叫猜测率guess意思是学生至少缺一个属性却蒙对了的概率。作答概率写成P(X_ij 1 | α_i) (1 - s_j)^{η_ij} · g_j^{1-η_ij}当η_ij 1全掌握概率是1 - s_j当η_ij 0有缺失概率是g_j。就这两行没有斜率参数没有连续潜变量非常干脆。这两个参数不是纯技术摆设它们有很强的诊断价值。一道题如果估出来s_j 0.05, g_j 0.08说明它区分度很好会的人基本做对不会的人基本做错。反过来如果g_j 0.35这道题就是个“送分题”不会的人也有三成多能蒙对它对诊断的贡献很有限。我在项目里会把s_j 0.3或g_j 0.3的题单独列出来要么改题要么重新审Q矩阵标注因为这种参数往往意味着“题目实际考的属性和你标的不一样”。还有一类情况值得警惕s_j和g_j同时很大比如都超过0.3。这时候题目几乎不携带属性信息1 - s_j - g_j接近0它在似然里就是个常数对后验概率没有影响。这类题留在测验里只会增加作答负担。2.2 EM算法实现E步到M步的推导DINA的标准估计方法是EM算法因为它有闭式M步实现起来不算复杂。整个流程是先给slip、guess和属性先验一个初值然后反复做两件事。E步对每个被试算出他在2^K种属性模式上的后验概率。给定模式α似然是Π_j P(x_ij|α)^{x_ij} (1-P(x_ij|α))^{1-x_ij}乘上模式的先验概率p(α)DINA默认假设属性间独立p(α) Π_k π_k^{α_k}(1-π_k)^{1-α_k}归一化就得到后验。M步用后验概率做权重更新参数。对题目j令W_ij P(η_ij 1 | x_i)也就是被试i在“该题涉及的属性全掌握”的所有模式上的后验概率之和。那么s_j Σ_i W_ij · (1 - x_ij) / Σ_i W_ijg_j Σ_i (1 - W_ij) · x_ij / Σ_i (1 - W_ij)直觉很好懂失误率就是“应该答对的人里实际答错的比例”猜测率就是“应该答错的人里实际答对的比例”。属性先验π_k的更新就是所有被试在属性k上后验概率的均值。这里有个工程上必须处理的细节s_j或g_j一旦被估到0取对数就会炸。所以我在实现里统一加伪计数把M步改成s_j (num a0) / (den a0 b0)a0 b0 1相当于给了一个Beta(1,1)的弱先验。这个改动看起来微不足道但它能把很多“参数跑到边界导致不收敛”的情况直接消掉。另一条路是硬截断到[0.001, 0.999]但这会让梯度信息丢失不如加伪计数干净。收敛判据我一般用对数似然的增量连续两次迭代的差值小于1e-6就停同时设一个max_iter 2000兜底。实践中大部分数据集在100到300次迭代内收敛跑到2000次还没收敛的基本都是数据或Q矩阵有问题不是算法问题。2.3 从属性后验到诊断报告模型跑完之后每个学生手里拿到的是一组后验概率P(α_ik 1 | x_i)k从1到K。这一步到“掌握/未掌握”的二值判定有两种常见做法区别不小。一种是取后验概率大于0.5就判掌握。这个做法简单但有个问题它忽略了属性之间的相关性可能给出一个数据上概率很高但模式上很别扭的结果。比如三个属性里前两个后验都是0.9第三个是0.45按阈值判就是“会会不”但模型可能认为“会会会”和“会会不会”这两个模式的联合后验更高。另一种是取联合后验概率最大的模式MAP。这个做法在模式层面自洽但有个已知的偏差MAP会偏向那些先验概率高的模式在属性掌握率不均衡的时候容易系统性误判。我做过一个对比当某个属性的掌握率只有20%时MAP判定会把一部分边缘案例全部推向“未掌握”。我现在的做法是两条腿走路对外报告用0.5阈值因为教研老师看得懂、好解释对内做模型评估和分类一致性检验时用MAP因为它是一个完整的模式。另外我会多输出一个“后验熵”熵高的学生说明模型对他的判断不确定这类学生适合推送更多题目再测一轮而不是直接给结论。报告呈现上光给“掌握/未掌握”没多大用。我会把属性按后验概率排序标出“确定掌握0.85”“倾向掌握0.6–0.85”“存疑0.4–0.6”“未掌握0.4”四档然后给一条推荐路径从存疑的属性开始补因为它的边际收益最高——已经掌握的不用管完全没掌握的补起来太费劲存疑的往往是差一点点。3. HO-DINA属性之间强相关时的高阶解法3.1 独立假设什么时候会崩DINA默认属性之间是独立的先验直接写成各属性掌握概率的连乘形式。这个假设在数学上让EM变得极其简单但在现实中经常不成立。数学能力下面的“运算”“推理”“空间想象”三个属性大概率是正相关的——运算强的学生推理通常也不弱。属性相关一强独立假设就会带来两个后果一是参数估计的方差被高估二是属性掌握模式的分布被扭曲模型会给出一些现实中很少见的模式。属性数越多、属性之间越同质这个问题越严重。我见过一个K6的项目属性之间两两相关普遍在0.5以上用标准DINA跑出来有将近三成的学生被分到了“部分掌握且模式很零散”的类别里教研老师一看就说不对劲——实际教学里学生的薄弱点通常是成片出现的不会这么碎。3.2 高阶结构的数学形式与识别条件HO-DINA的思路是既然属性相关那就假设它们背后有一个共同的连续高阶特质θ可以理解成“总体学业能力”属性掌握与否由这个特质决定。常见的有两种连接形式一种是阈值形式α_ik 1当且仅当θ_i ≥ c_k。每个属性有一个掌握阈值c_kθ超过阈值就算掌握。这个形式直观但c_k的尺度需要和θ一起识别。另一种是logistic形式也是我更常用的P(α_ik 1 | θ_i) 1 / (1 exp(-(λ_k θ_i β_k)))。λ_k是属性k在高阶特质上的载荷反映这个属性和总体能力的关联强度β_k是截距控制该属性的基础掌握率。θ服从标准正态分布均值固定为0、方差固定为1这是识别所必需的——不固定的话λ和θ的尺度可以同时放大缩小模型不可识别。在这个结构下属性之间的边际相关由共享的θ自动产生不需要额外估相关矩阵。原来DINA里隐含的2^K - 1个模式概率被压缩成K个载荷加K个截距参数大幅减少小样本下稳定性明显更好。识别条件上还有几条必须满足属性数K至少为3K2时高阶结构和一个双因子模型没区别估了也没意义每个属性至少要有一道题在测这是老生常谈λ_k的符号需要固定通常约束为正否则会出现符号翻转的等价解。3.3 什么时候不该用HO-DINAHO-DINA不是万能升级包用错场景比不用更糟。第一种不该用的情况属性分属不同的能力维度。比如一个测评同时测“语言理解”和“空间推理”这两个属性背后不是一个共同特质硬套单高阶因子会强行把它们拉成相关结果是两个维度的信息互相污染。判断方法很简单先跑独立DINA看属性间的残差相关如果明显分成两簇就该考虑多维度高阶模型而不是单因子。第二种属性数太少。K3的时候高阶结构省不下多少参数反而增加了积分计算的开销和收敛难度收益不明显。第三种样本量本来就充足比如5000以上、属性相关也不强。这种情况下独立DINA估得很稳上高阶结构属于过度建模还会让结果解释变复杂——你得同时解释高阶特质分数和属性掌握状态教研老师不一定买账。实操上HO-DINA的估计需要在EM的E步里对θ做数值积分常用高斯-埃尔米特积分积分点数取21或31就够再多边际收益很低。这意味着计算量大约是标准DINA的20到30倍K6、N2000的数据在普通笔记本上要跑十几分钟。我会先用独立DINA跑一遍看结果合理性确认数据没问题再上高阶模型。4. P-DINA多级评分题目的处理方案4.1 二值化硬套会丢什么真实测评里不可能全是0/1题。量表题是4级或5级评分主观题有步骤分操作题有等级评定。最省事的做法是设个阈值二值化比如4级量表把3和4算1、1和2算0。这个做法我用过能用但丢的东西比想象中多。丢的第一样是信息量。一个学生答“完全符合”和另一个答“比较符合”二值化之后都是1模型看不出差别。而在多级模型里这两个人的属性后验概率是不一样的前者更接近确定掌握。丢的第二样是题目的区分能力。多级题的价值恰恰在于中间类别它能把“半懂不懂”的学生和“完全不懂”的学生分开。二值化把这个中间地带一刀切没了题目的诊断效率下降得很明显。我做过对比同一批量表题二值化之后题目的1 - s - g平均掉了将近三成。丢的第三样是类别之间的有序性。二值化之后你没法检验“选3的人比选2的人掌握概率更高”这个假设而这条假设如果不成立说明题目本身有问题是需要被发现的。4.2 累积类别建模与单调约束P-DINA多级计分DINA的标准做法是累积类别形式。假设第j题有m_j 1个类别0到m_j对每个切分点c定义累积概率P(X_ij ≥ c | α_i) g_jc (1 - s_jc - g_jc) · η_ijc 1, 2, …, m_jη_ij的定义和二值DINA完全一样。每个切分点有自己的失误率s_jc和猜测率g_jc。类别概率由累积概率相减得到P(X_ij c) P(X_ij ≥ c) - P(X_ij ≥ c1)边界上P(X ≥ 0) 1、P(X ≥ m_j1) 0。这个形式要成立必须满足单调约束0 ≤ g_j1 ≤ g_j2 ≤ … ≤ 1 - s_j。翻译成人话就是“掌握程度越高答到高类别的概率越大”这是有序性的数学表达。估计的时候如果不加这个约束很容易出现交叉模型会给出“掌握的学生反而更容易选低类别”这种荒谬结果。参数个数上一道4级题m3要估3个slip加3个guess共6个参数是二值题的3倍。所以多级题的样本需求也水涨船高我的经验是每道多级题至少要150到200个有效作答才能估稳。类别数越多越吃样本5级以上量表如果没有上千人的数据估计结果会很飘。4.3 步骤法工程上最省事的替代方案如果不想上完整的多级模型有个折中方案我一直推荐步骤法也叫逐级分解。把一道m1级的题拆成m个虚拟的二值“步骤题”每一步问“有没有达到这个等级或以上”所有步骤共用同一行Q矩阵然后直接跑二值DINA。举个例子一道4级题拆成三个二值题达到1级及以上、达到2级及以上、达到3级及以上。这样一道题变成三道题参数从6个变成6个每个虚拟题2个看起来没省但实现上完全复用了二值DINA的代码不需要改任何估计逻辑工程成本几乎为零。这个方案的代价是引入了“步骤间局部独立”的假设也就是在给定属性掌握状态后三个虚拟题的回答互相独立。现实中这个假设经常被违反——第一步做对了第二步做对的概率天然更高。违反的后果是标准误被低估参数估计看起来比实际更精确。我的实际经验是类别数不超过3也就是4级及以下量表、题目本身是累积性质比如“能完成基本操作”“能独立完成”“能优化改进”步骤法跑出来的结果和完整P-DINA差别很小属性掌握判定的符合率通常在90%以上。类别数更多或者题目不是严格累积的老老实实上完整模型。还有一类题目要特别处理主观题的步骤分。这种分数往往不是“等级”而是“采分点累加”两个学生都是3分但一个是在第一步拿了3分另一个是三步各拿1分掌握状态可能完全不同。这种数据我建议不要直接当多级变量喂进去而是把采分点拆成多个二值题分别建模信息损失最小。5. G-DINA一个框架统一四种模型5.1 饱和形式的参数含义G-DINA的思路是不预设属性之间是“与”还是“或”的关系而是让数据自己说话。对第j题设它涉及K*_j个属性把这些属性的所有主效应和交互项都放进线性组合P(α*_l) δ_j0 Σ_k δ_jk α_lk Σ_{kk} δ_jkk α_lk α_lk … δ_j(12…K*) Π_k α_lk其中α*_l是第l种属性模式在这道题涉及属性上的取值组合共2^{K*_j}种。δ_j0是截距δ_jk是属性k的主效应δ_jkk是两两交互最后一项是最高阶交互。这个式子本质上是把2^{K*_j}个作答概率用一个饱和的参数化写出来参数个数正好等于模式数所以它是完全饱和的拟合能力最强。但参数一多可解释性就下来了——你很难跟教研老师解释“属性1和属性2的交互项系数是0.13”是什么意思。G-DINA真正的价值不在饱和形式本身而在于它是一个母框架。通过对δ施加不同约束可以退化出一大批经典模型这样模型比较就在同一个框架、同一套估计代码、同一个拟合指标下进行避免了跨软件比较带来的口径差异。5.2 约束如何退化出DINA、DINO、RRUM、ACDM下面这张表是我自己整理对照用的把这几个模型在G-DINA里的约束关系列清楚模型联结假设在G-DINA中的约束DINA非补偿全掌握才会主效应全为0所有低阶交互为0仅δ_j0 g_j、最高阶δ_j(12…K*) 1 - s_j - g_jDINO补偿会一个就行主效应均等于1 - s_j - g_j各阶交互按阶数交替取负号RRUM非补偿的乘法形式主效应δ_jk π_j(r_jk - 1)所有交互项为0ACDM补偿的加法形式所有交互项为0主效应自由估计G-DINA无约束全部参数自由DINA的约束最狠它把一道题的所有信息压缩成两个自由度所以它最稳但也最不灵活。DINO适合那种“多知识点交叉、答对一个就能过”的题比如综合应用题。RRUM和ACDM介于两者之间ACDM的每个属性主效应可以不同相当于允许不同属性对答对的贡献权重不一样。模型选择上如果一个测验里不同题目遵循不同的联结假设完全可以混着用——G-DINA框架允许逐题指定模型。我做过一个混合测评其中三分之二是标准知识点题用DINA三分之一是综合应用题用DINO整体拟合比全用DINA好一截AIC降了将近200。5.3 参数爆炸与模型选择的实操建议G-DINA的代价很直接一道题如果涉及4个属性饱和形式就有2^4 16个参数。20道题、每题涉及3个属性参数总量就是160个再加上属性分布参数样本没有两三千根本估不动。而且参数一多边界解、局部极值的概率大幅上升同一个数据集跑两次不同初值结果可能差很远。所以我的实操原则是优先用约束模型G-DINA只在需要检验“联结假设是否成立”的时候用。具体流程是三步。第一步全部题目先按理论假设指定DINA跑一遍看整体拟合。如果绝对拟合指标SRMSR、RMSEA已经在可接受范围就没必要上G-DINA直接收工。第二步如果整体拟合不好用itemfit看是哪些题拖后腿。只对那些拟合差的题目放开约束逐题尝试DINA、DINO、ACDM、RRUM四种约束看哪个的BIC最低。不要一上来就全放开那是给自己找麻烦。第三步如果逐题放开之后还有明显 misfit再考虑对个别题目用饱和G-DINA。这时候要接受一个事实这几个参数估计出来主要是为了拟合不要过度解读它们的实际含义。模型比较的指标我一般看三个BIC或SABIC惩罚参数个数越小越好、SRMSR标准化残差均方根小于0.05算好、RMSEA小于0.05算好。AIC在样本量大时会过度倾向复杂模型我不太用。另外一定要看参数的合理性——如果放开约束之后估出s_j 0.45这种值说明这道题本身有问题而不是模型选错了这时候该去改题。6. 工程落地全流程6.1 数据前处理的六道关卡模型跑不起来九成问题出在数据。我现在固定跑六道检查一道不过就不往下走。第一道编码一致性。作答数据里除了0和1经常混着9缺失、2多选、-1未作答这些编码。DINA只认0和1其他值必须显式处理。缺失值要么删被试、要么删题目、要么用EM在似然里积分掉三条路选一条别让它们悄悄进模型。第二道题目对齐。作答数据的列名和Q矩阵的行名必须一一对应而且顺序要一致。我现在的做法是强制用题目ID做索引读进来第一件事就是stopifnot(all(colnames(resp) rownames(Q)))不通过直接报错。这一行代码帮我省过至少三次通宵排查。第三道零方差题。某道题所有人都答对或都答错它的slip或guess没法估必须剔除。全对的题在某些场景下是“必答题”但对诊断零贡献。第四道被试有效性。作答题目数少于总题数一半的、全部选同一个选项的、作答时间异常短的这些记录会严重干扰参数估计我一般先剔掉再建模剔了多少在报告里说明。第五道属性掌握率的合理性。如果某个属性估计出来的掌握率低于5%或高于95%要么是Q矩阵标错了要么是样本本身有偏这时候要回去查不要接受这个结果。第六道样本量与K的匹配。回到第1.3节那张表对一遍不匹配就先降K别硬跑。6.2 手写EM实现与R包调用如果只是要结果R生态里的成熟包已经够用。数据准备好之后标准DINA的调用大概是这样library(CDM) resp - read.csv(responses.csv, row.names 1) # 行被试列题目 Q - as.matrix(read.csv(q_matrix.csv, row.names 1)) # 行题目列属性 # 对齐检查这一步千万别省 stopifnot(all(colnames(resp) rownames(Q))) stopifnot(all(resp %in% c(0, 1))) fit_dina - din(resp, q.matrix Q, maxit 2000, conv.crit 1e-6, progress FALSE) summary(fit_dina) round(fit_dina$slip, 3) round(fit_dina$guess, 3) head(fit_dina$person) # 每个被试的属性后验概率 fit_dina$attribute.patt # 掌握模式的分布高阶模型和多级评分换个包更顺手高阶结构通过属性分布参数指定library(GDINA) # HO-DINA属性由高阶连续特质驱动 fit_ho - GDINA(dat resp, Q Q, model DINA, att.dist higher.order, nstarts 30) # 多级评分直接把有序整数矩阵传进去模型指定为DINA约束 fit_poly - GDINA(dat resp_poly, Q Q, model DINA, nstarts 20)nstarts这个参数不要设太小。DINA的似然面有局部极值尤其在高阶模型和多级模型里单次初值跑出来的结果可能离全局解很远。我一般设20到30次随机初值取似然最大的那组。如果团队用的是Python技术栈或者需要对EM过程做定制比如加自定义先验、加约束、做增量更新手写一版更可控。下面是我常用的最小实现能跑通、能复现可以直接当基线import numpy as np def build_patterns(K): 生成 2^K 种属性掌握模式 return np.array([[(i k) 1 for k in range(K)] for i in range(2**K)], dtypeint) def eta_matrix(Q, pats): eta[a, j] 1 表示模式 a 在题目 j 涉及的属性上全部掌握 # Q: J x K, pats: 2^K x K - 2^K x J return np.prod(pats[None, :, :] ** Q[:, None, :], axis2).T def em_dina(X, Q, K, max_iter2000, tol1e-6, a01.0, b01.0, seed0): rng np.random.default_rng(seed) N, J X.shape pats build_patterns(K) eta eta_matrix(Q, pats) # 2^K x J L len(pats) S rng.uniform(0.05, 0.25, J) # slip G rng.uniform(0.05, 0.25, J) # guess pi np.full(K, 0.5) # 属性先验掌握率 old_ll -np.inf for it in range(max_iter): P (1 - S)[None, :] ** eta * G[None, :] ** (1 - eta) # 2^K x J # ---- E 步 ---- logP np.log(np.clip(P, 1e-12, 1)) logQ np.log(np.clip(1 - P, 1e-12, 1)) loglik X logP.T (1 - X) logQ.T # N x 2^K logprior (pats * np.log(pi) (1 - pats) * np.log(1 - pi)).sum(axis1) post loglik logprior[None, :] # N x 2^K post - post.max(axis1, keepdimsTrue) post np.exp(post) post / post.sum(axis1, keepdimsTrue) # N x 2^K ll np.log(np.clip((post * np.exp(loglik logprior[None, :] - loglik.max(axis1, keepdimsTrue))).sum(axis1), 1e-12, 1)).sum() # ---- M 步 ---- for j in range(J): w post * eta[:, j][None, :] # N x 2^K den_s w.sum() den_g post.sum() - den_s num_s (w * (1 - X[:, j])[:, None]).sum() num_g ((post - w) * X[:, j][:, None]).sum() S[j] (num_s a0) / (den_s a0 b0) G[j] (num_g a0) / (den_g a0 b0) pi (post.T pats) / N pi np.clip(pi, 1e-3, 1 - 1e-3) if abs(ll - old_ll) tol: break old_ll ll person_prob post pats # N x K 边际后验 return {slip: S, guess: G, pi: pi, person_prob: person_prob, pattern_post: post, loglik: old_ll, iter: it 1}这段代码有三处是我特意加的不是可有可无np.clip防止log(0)post - post.max(...)是标准的log-sum-exp技巧属性数一大似然连乘会下溢成0不减最大值整个后验就废了M步里的a0/b0就是前面说的伪计数没有它参数跑到0之后下一次迭代直接nan。跑之前建议先用模拟数据验一遍确认代码能把真参数找回来K, J, N 4, 20, 1000 rng np.random.default_rng(42) Q np.zeros((J, K), dtypeint) for j in range(J): Q[j, rng.choice(K, size2, replaceFalse)] 1 s_true rng.uniform(0.05, 0.20, J) g_true rng.uniform(0.05, 0.20, J) alpha (rng.random((N, K)) 0.6).astype(int) eta np.prod(alpha[:, None, :] ** Q[None, :, :], axis2) # N x J prob (1 - s_true)[None, :] ** eta * g_true[None, :] ** (1 - eta) X (rng.random((N, J)) prob).astype(int) res em_dina(X, Q, K) print(平均|slip误差|:, np.abs(res[slip] - s_true).mean()) print(平均|guess误差|:, np.abs(res[guess] - g_true).mean())在N1000、J20、K4这个配置下平均绝对误差通常能压到0.03以内。如果误差明显超出这个量级先别怀疑算法去查Q矩阵和模拟代码里的属性生成逻辑通常是这两处对不上。6.3 拟合检验与结果报告模板模型跑完只是开始得证明它靠谱。我固定的检验清单是四层。第一层绝对拟合。SRMSR和RMSEA这两个指标最直观SRMSR小于0.05、RMSEA小于0.05算拟合良好0.05到0.08之间是可接受超过0.08就要去看是哪些题的问题。MADcor和MADQ3作为辅助主要看残差相关的整体水平。第二层相对拟合。AIC、BIC、SABIC、CAIC放在一起看重点看BIC和SABIC它们在样本量大时更保守。同一批数据下比较DINA、DINO、ACDM、G-DINA选BIC最低且参数可解释的那个。第三层题目拟合。逐题算χ²或者用参数化的itemfitp值小于0.05的题目单独拉出来审。注意多重比较问题20道题按0.05算平均会有一道假阳性所以我会用Bonferroni校正后的阈值0.05/20 0.0025来判定避免被单道题带偏。第四层分类一致性。这个最容易被忽略但最重要。做法是用bootstrap抽一批样本重跑模型看同一个学生的属性判定在两次之间是否一致。一致率低于0.85的属性说明它的测量精度不够报告里要标注“仅供参考”不能让教研老师拿去做教学决策。报告我一般输出三张表
阅读完成 · 觉得有帮助？
