mcem r语言代码_处理潜变量或随机效应的高维积分时，我们做些什么之MCEM

本文面向对象是心理和教育等社科统计的初入门者。

广义线性潜变量模型（Generalized Linear Latent Variable Models，在心理学和教育学常用的是连接函数为probit或logit的项目反应模型或因子分析模型）和广义线性混合模型（Generalized linear mixed model, GLMM，在心理和教育领域又称之为多层线性模型）是心理和教育统计的基础入门模型。

初学的人应该明白，这两个模型在参数估计的过程中会遇到难以处理的积分（intractable integral）。例如广义线性潜变量模型，EM算法的Q函数为

其中

是观察变量，

是斜率，

是截距，

是潜变量，

是似然函数，

是条件概率密度函数。

潜变量或维度较低的时候，传统的高斯厄米特积分（Gauss Hermite quadrature）足以胜任，在《Item Response Theory: Parameter Estimation Techniques》和《Generalized, Linear, and Mixed Models》这两本书中，也是推荐的高斯厄米特积分（在全息项目因子分析（多维项目反应理论）：参数估计算法这篇文章中，用的也是高斯厄米特积分），但在高维情况下，高斯厄米特积分的设计矩阵会疯狂增长，以至于我们单机内存会不够用（如果我们用5个积分点估计10个维度的高维积分，我们的设计矩阵大小为

，这个矩阵元素数量为9700万，10个积分点估计10个维度，设计矩阵的元素数量将会是1000亿），如果在生产环境部署高斯厄米特积分代码处理高维积分，无异于是拿系统稳定性开玩笑（就算是分布式系统也不要尝试），所以在生产环境中，不仅需要我们用C++、cuda、并行处理或分布式计算加速我们的计算，还需要我们寻找一种高效的处理高维积分的方法。

而处理GLLVM和GLMM模型的高维积分方法，主要有下面几种：

马尔科夫链蒙特卡洛算法（MCMC）。MCMC算法是实现起来最简单的算法，堪称无脑式方法，给初中生培训两三天估计就能上手，优点是实现简单并且对初值不敏感参数估计结果稳定，缺点是速度慢且占内存；

蒙特卡洛期望极大算法（MCEM）。MCEM方法是把EM算法中的难以处理的E步积分，用蒙特卡洛方法求解，优点与mcmc一样实现简单，运行速度比mcmc快很多倍，缺点是某些情况下容易梯度爆炸并且难以判断是否收敛；

拉普拉斯近似（Laplace Approximation, LA）和自适应高斯厄米特积分（Adaptive Gauss Hermite quadrature, AGQ）。拉普拉近似和自适应高斯厄米特积分这两个方法本质上一个方法，拉普拉斯近似其实是一个积分点的自适应高斯厄米特积分，这两个方法也约可以看做是加了重要性抽样并且固定了样本值的MCEM方法，这两个方法其实都用了极大似然估计中点估计服从协方差矩阵为信息矩阵逆矩阵的渐进正态分布的性质或正交因子模型中因子协方差矩阵的性质，优点均是计算简单；

变分近似（variational approximation, VA）算法。变分近似基于最大化ELOB，原理是找一个容易计算的普通函数近似难以处理的积分函数；

期望传播（expectation propagation）算法。期望传播算法基于最小化KL距离，期望传播和变分近似在数学形式上很像，区别是期望传播走的是直接最小化KL距离的路线，并且期望传播算法看上去更像是矩估计；

其他方法。还有一些非主流的方法，例如EM算法的变体MHRM算法，这是一种基于类似于随机梯度下降和MH采样的算法。

这些方法已经广泛应用各个统计软件，例如R语言的package中，mirt采用了MCEM和MHRM，lme4之前的版本采用了LA和AGQ，最近版本好像只有AGQ，mcmcpack包含有处理高维irt的模块，还有一些小众的包采用了MCEM和VA。

上面这张方法，无论是MCMC还是MCEM，还是AGQ，LA，VA，抑或是MHRM，你都可以看成是一种坐标下降法(coordinate descent)，第一个步骤总是求隐变量，第二个步骤求实体参数。

这一章主要聊MCEM。关于MCEM，推荐阅读《The EM Algorithm and Extensions》后面关于MCEM的章节。下面先以项目因子分析模型为例（连接函数为probit或logit的广义潜变量模型，另一种名称是多维项目反应理论）。

最无脑的方法

最无脑的做法是直接套用Bock和Aitkin的解法，Bock和Aitkin的解法是对带有积分的Q函数求梯度，然后将问题简化为求潜变量的后验条件概率密度。

E步

对带有积分形式的Q函数直接求score函数（偏导数），并经过一系列的计算和转换（详细请看《Item Response Theory: Parameter Estimation Techniques》），问题简化为求下式，即

的后验条件概率密度

上式中

是

的先验概率密度函数，

的先验信息默认为