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

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

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

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

$Q ={\int L(u, a, b, \theta)}g(\theta|u, a,b)d\theta$

其中

$u$ 是观察变量，

$a$ 是斜率，

$b$ 是截距，

$\theta$ 是潜变量，

$L$ 是似然函数，

$g$ 是条件概率密度函数。

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

$9765625\times10$ ，这个矩阵元素数量为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》)，问题简化为求下式，即

$\theta_i$ 的后验条件概率密度

$P(\theta_i|u, a, b) = \frac{L(u, a, b, \theta_i)g(\theta_i)}{\int L(u, a, b, \theta_i)g(\theta_i)d\theta}$

上式中

$g$ 是

$\theta_i$ 的先验概率密度函数，

$\theta_i$ 的先验信息默认为

$\theta_i\sim N(0, I)$ ，

$L$ 是似然函数，当我们使用高斯厄米特积分的时候，

$\theta_i$ 的取值是高斯厄米特的节点(node)，同时概率密度

$g$ 也相应的是高斯厄米特积分的权重(weight)。

注意到

$\int L(u, a, b, \theta_i)g(\theta_i)d\theta = E(L(u, a, b, \theta_i))$ ，而先验信息

$\theta\sim N(0, I)$ ，所以我们可以随机的从标准多元正态分布中采集

$K$ 个

$\theta$ 的样本，所以

$\int L(u, a, b, \theta_i)g(\theta_i)d\theta \approx \sum_{k=1}^{K}{L(u, a, b, \theta_i))\frac{1}{k}}$

所以

$g$ 的概率密度值为

$\frac{1}{k}$ ，所以

$P(\theta_k) = \frac{\frac{1}{K}L(u, a, b, \theta_k)}{\frac{1}{K}\sum_{k=1}^{K}L(u, a, b, \theta_k)}= \frac{L(u, a, b, \theta_k)}{\sum_{k=1}^{K}L(u, a, b, \theta_k)}$

M步

M步的形式和高斯厄米特积分的M步一致，可以用irls或牛顿迭代或BFGS等求解。

著名的R包mirt的MCEM方法用的就是上述的形式。这种形式的参数估计问题是很不稳定，非常容易出现类似梯度爆炸的问题，另外难以判断收敛(下图比较了LSAT数据集下的MCEM的收敛性和GQ(高斯厄米特积分)的收敛性)。优势是代码实现简单，非常容易推广到其他的IRT模型。MCEM的收敛不平稳GQ的收敛很平稳

最无脑方式的对数似然值趋势

最弱智的方法

另一种方法是利用极大似然估计值渐进多元正态分布的性质。即

$\theta_i \sim N(\hat{\theta_i}, \hat{\Sigma_i}^{-1})$

$\hat{\theta_i}$ 是

$\theta_i$ 的估计，

$\Sigma_i$ 是信息矩阵，利用这个分布抽取

$K$ 个

$\theta_i$ 样本，然后求极大

$Q = {\frac{1}{K}\sum_{k=1}^{K}\sum_{j=1}^{J}\sum_{i=1}^{N}logP(u_{ij}, a_{j}, b_{j}, \theta_{ij})}$

对于二分类模型

$P(u_{ij}, a_{j}, b_{j}, \theta_{ik}) =P(a_j, b_j, \theta_{ik})^{u_{ij}}(1 - P(a_j, b_j, \theta_{ik}))^{1 - u_{ij}}$

这种方法是最弱智的方法。通常不会有人使用这种方法进行mcem计算，已经计算了

$\hat\theta$ 和

$\hat\Sigma$ ，为什么不直接使用LA和AGQ方法。另外，这种方法其实就是数据增强的坐标下降法(也即联合极大似然法)，聊胜于无。

最有趣的方法

最有趣的采样方法是Stephen Schilling提出的，Schilling参考的是1991年的一篇讲MCEM贝叶斯二分类probit估计的文章，Schilling原文有些公式没有给出推导，可能会误导初学者，我在下面给出一些详细的推导。

对于二分类probit回归(binary probit progression)，我们有下面的公式

$P(u_{ij} =1 | \theta_i, a_j, b_j) = P(x_{ij} > 0|\theta_i, a_j, b_j)$ =

$\Phi(\theta_i^Ta_j + b_j)$

所以依据广义线性模型对方差量纲的限制(设定为1)，

$x_{ij}$ 服从截断正态分布

$x_{ij} \sim N(\theta_i^Ta_j + b_j, 1)$ ，如果

$u_{ij} = 1$ ，

$x_{ij} \geq 0$ ，如果

$u_{ij} = 0$ ，

$x_{ij} < 0$

又带有均值结构的因子分析模型为

$X_i = B + \theta_i^TA + E$

其中

$X_i$ 和

$\theta_i$ 都是潜变量，

$A$ 是因子载荷，

$E$ 是误差，由于probit回归设定了误差方差为单位矩阵

$I$ ，所以

$E \sim N(0, I)$ ，在因子分析中，通常采用tetrachoric相关矩阵计算因子载荷

$A$ 和均值

$B$ ，由于残差方差固定为了单位矩阵

$I$ ，这种计算方法在著名的潜变量计算软件mplus中又叫delta法。

注意到先验信息

$\theta_i \sim N(0, I)$ ，所以该因子模型为正交因子模型，因子方差协方差矩阵除对角线外其他元素为0，对角线元素为1，所以

$E(X_i) = B + AE(\theta_i) + E(E) = B$

$VAR(X) =E((X - B)(X - B)^{T}) = E((\theta_i^TA + E)(\theta_i^TA + E)^{T})=E(A^T\theta_i\theta^TA + A^T\theta_i E^T + E\theta_i^TA + EE^T)$

注意到

$E(\theta \theta^T) = VAR(\theta)=I$ ，

$E(E) = E(E^T) = 0$ ，

$E(EE^T) = VAR(E) = I$ ，所以

$VAR(X) = A^TA + I$

由于probit回归中潜变量

$X_i$ 服从正态分布，所以

$X\sim N(B, A^TA+I)$

若把

$\theta_i$ 当做变量

则

$A^TE(\theta_i) = X_i - B$ ，

$AA^TE(\theta_i)=A(X_i - b)$ ，所以

$E(\theta_i) = (AA^T)^{-1}A(X_i - B)$

上面公式严格来说应该通过最小二乘法或极大似然法导出，注意到

$\theta_i$ 拥有先验信息

$\theta_i \sim N(0, I)$ ，这相当于为我们之前因子公式的最小二乘形式加了一个

$L_2$ 正则项，所以我们需要把上面的公式改写成岭回归的形式(或者MAP估计的形式)，即

$E(\theta_i) = (AA^T + I)^{-1}A(X_i - B)$

同样由岭回归最小二乘法参数估计的性质，我们可以得到

$VAR(\theta_i) = (AA^T + I)^{-1}$

所以

$\theta_i \sim N((AA^T + I)^{-1}A(X_i - B), (AA^T + I)^{-1})$

于是我们找到了

$X_i$ (注意到

$X_{ij}$ 服从截断分布，所以

$X_i$ 也服从截断分布) 和

$\theta_i$ 的分布

$X_i\sim N(B, AA^T+I)$ ，如果

$u_{ij} = 1$ ，

$x_{ij} \geq 0$ ，如果

$u_{ij} = 0$ ，

$x_{ij} < 0$

$\theta_i \sim N((AA^T + I)^{-1}A(X_i - B), (AA^T + I)^{-1})$

由于

$X_i$ 的分布不依赖于

$\theta_i$ ，我们不需要通过像吉布斯采样那样的算法即可简单方便的采集

$(X_i, \theta_i)$ 联合分布的样本。但是

$X_i$ 的分布是一个截断多元正态分布，技术实现较为困难，自论文发表20多年以来，工业界和学术界均没有一个快速的解决截断多元正态分布采样的方法。所以我们需要借助IRT局部独立性这个假设，把截断多元正态分布变成多个截断正态分布，即

$x_{ij} \sim N(\theta_i^Ta_j + b_j, 1)$ ，如果

$u_{ij} = 1$ ，

$x_{ij} \geq 0$ ，如果

$u_{ij} = 0$ ，

$x_{ij} < 0$

由于

$x_{ij}$ 的分布依赖于

$\theta_i$ ，所以我们可以采用吉布斯算法采集

$(X_i, \theta_i)$ 联合分布的样本，如果对吉布斯算法不了解，可以参考《Monte Carlo Statistical Methods》这本入门书，最简单的理解就是把吉布斯算法想象成坐标下降法(貌似万事万物都可以理解成坐标下降法)。

所以我们采样算法最终如下

依据先验信息

$\theta_i \sim N(0, I)$ 采集样本

$\theta_i$

步骤一：依据多个截断正态分布采集

$X_i$

步骤二：依据岭回归计算

$\theta_i$ 的均值和方差，在多元正态分布下采集

$\theta_i$

重复步骤一和步骤二直到合适时候的结束。Schilling论文表示处理LSAT数据集，需要燃烧前10次迭代，然后每隔5次迭代采样一个样本，总共采集25个样本即可。

吉布斯抽样的Python代码

from numpy.dual import svd

from scipy.stats import truncnorm

import numpy as np

def multivariate_normal(mean, cov):

x = np.random.normal(0, 1, mean.shape)

u, s, v = svd(cov)

x = np.dot(x, np.sqrt(s)[:, None] * v)

x += mean

return x

def gibbs(scores, slop, threshold, point_size, burn, thin, collect_only_theta=True):

person_size = scores.shape[0]

dim = slop.shape[0]

if collect_only_theta:

theta_ar = np.zeros((point_size, dim, person_size))

else:

# 线性模型求解的充分统计量

t_z = t_x = z_t_z = z_t_x = 0

theta = np.random.normal(size=(dim, person_size))

# theta方差

v = np.linalg.inv(np.dot(slop, slop.T) + np.eye(dim))

for i in range(point_size * thin + burn):

z_val = np.dot(theta.T, slop) + threshold

lower = np.zeros_like(scores, dtype=np.float)

lower[scores == 0] = -np.inf

lower[scores == 1] = -z_val[scores == 1]

upper = np.zeros_like(scores, dtype=np.float)

upper[scores == 1] = np.inf

upper[scores == 0] = -z_val[scores == 0]

# 采样x

x = truncnorm.rvs(lower, upper, loc=z_val)

# theta均值

loc = np.dot(np.dot(v, slop), (x - threshold).T)

# 采样theta

theta = multivariate_normal(loc.T, v).T

if i >= burn and (i - burn) % thin == 0:

if collect_only_theta:

theta_ar[int((i - burn) / thin)] = theta

else:

t_z += np.sum(theta, keepdims=True)

t_x += np.sum(x, axis=0, keepdims=True)

z_t_z += theta.dot(theta.T)

z_t_x += theta.dot(x)

if collect_only_theta:

return theta_ar

else:

return t_z / point_size, t_x / point_size, z_t_z / point_size, z_t_x / point_size

为了能向量化的进行多元正态分布采样，我自己写了一个进行多元正态分布随机数采样的方法，numpy原生方法不支持均值向量化操作。

采样代码

# 用logit生成模拟值

slop = np.random.uniform(0.7, 2, (2, 10))

threshold = np.random.normal(0, 1, size=(1, 10))

theta = np.random.normal(0, 1, size=(1000, 2))

z = np.dot(theta, slop) + threshold

p = 1 / (1 + np.exp(-z))

scores = np.random.binomial(1, p)

# 初值

slop_start_val = np.ones_like(slop)

threshold_start_val = np.zeros_like(threshold)

# 采样

gibbs(scores, slop_start_val, threshold_start_val, point_size=2000, burn=100, thin=1, collect_only_theta=True)

我们把采样的结果画成图采样呈现平稳分布

M步

M步的两种方法，一是只使用采集的

$\theta_i$ 样本集进行Q函数的极大计算，二是使用

$(X_i, \theta_i)$ 两个变量的联合样本集或只使用

$X_i$ 样本集，由于

$x_{ij} \sim N(\theta_i^Ta_j + b_j, 1)$ ，我们可以使用非常简单的最小二乘法求解

$a_j$ 和

$b_j$ 。

第一种方法，还是对完全数据似然函数求期望

$Q = {\frac{1}{K}\sum_{k=1}^{K}logL(a, b, \theta)}={\frac{1}{K}\sum_{k=1}^{K}\sum_{j=1}^{J}\sum_{i=1}^{N}(u_{ij}logP_{ijk} + (1-u_{ij})log(1- P_{ijk})})$

其中

$P_{ijk} = \Phi(a_j, b_j, \theta_{ik})$ ，

$\Phi$ 是probit函数，

$i$ 指代person，

$j$ 指代item，

$k$ 指代采样的样本，为了节省内存同时加快速度，我们一般不对Q函数直接求极大，而对每一个

$Q_{j}$ 求极大，即对

$Q_j={\frac{1}{K}\sum_{k=1}^{K}\sum_{i=1}^{N}(u_{ij}logP_{ijk} + (1-u_{ij})log(1- P_{ijk})})$ 求极大，这种方法一看就知道E步起的作用是数据增强，并且这种方法类似于最弱智的MCEM方法。

第二种方法非常的easy，直接套用最小二乘公式即可求解

$A$ ,

$B$ 。对于第二种方法，在我们采样的时候，没必要保存每一个采样的样本值，只需要计算参数估计的充分统计量期望即可，即

$E(\theta^T\theta)$ ,

$E(\theta^TX)$ ,

$E(1^T\theta)$ ,

$E(1^TX)$ ，我们可以使用

$(X_i, \theta_i)$ 或只使用

$X_i$ 计算充分统计量，然后依据这些充分统计量的期望值，计算

$A$ ,

$B$

$A=\frac{(E(\theta^T\theta)−\frac{E(1^T\theta)^TE(1^T\theta)}{n})^{-1}}{E(\theta^TX) − \frac{E(1^T\theta)^TE(1^TX)} { n}}$

$B = \frac{E(1^TX) − E(1^T\theta)A}{n}$

Schilling方法的优势是很稳定，劣势是吉布斯采样的速度问题，虽然可以用并行采样规避这个问题，收敛准则难以判断(Schilling原文给出了一种判断mcem收敛的方法，有兴趣的可以看看，下图是估计LSAT数据集的斜率变化过程，用的是第一种方法，即求Q函数极大)。LSAT数据集斜率参数估计

Schilling方法推广到ordinal数据模型的时候，需要进行一些简单的参数化的工作，Schilling方法还能非常容易的推广到动态IRT模型和层级IRT模型等等复杂的IRT的参数估计上，总之Schilling方法是非常棒的一种参数估计方法。

三种方法的本质

三种方法的本质都是求潜变量的条件概率密度。其中最弱智方法和最有趣方法的本质也是数据增强，所以这两种方法的条件概率密度都是

$\frac{1}{K}$ (

$K$ 是每个样本集上的采样数)，这两种方法的参数估计也较为稳定，而最无脑方法采用的是计算后验概率密度，由于每次采样的潜变量值不一样，导致这种方法容易梯度爆炸，所以我也很疑惑为什么R包mirt会采用这种方法。

随机效应示例

稍后补上

代码示例

稍后补上

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

mcem r语言代码_生态学数据处理常用R语言代码
使用R来处理生态学数据越来越受到科研工作者的青睐,语义编程风格.漂亮的出图效果,能直接俘获众多用户.本文将生态学数据处理中经常会使用到的功能做个搜集整理. 本文假设读者有一些R的基础知识,对于R的编程 ...
mcem r语言代码_一个简单文本分类任务-EM算法-R语言
一.问题介绍概率分布模型中,有时只含有可观测变量,如单硬币投掷模型,对于每个测试样例,硬币最终是正面还是反面是可以观测的.而有时还含有不可观测变量,如三硬币投掷模型.问题这样描述,首先投掷硬币A,如 ...
mcem r语言代码_R语言阈值自回归模型（TAR）代码示例
原文链接: R语言时间序列TAR阈值模型分析tecdat.cn 阈值模型用于统计的几个不同区域,而不仅仅是时间序列.一般的想法是,当变量的值超过某个阈值时,过程可能表现不同.也就是说,当值大于阈值时 ...
mcem r语言代码_R语言面向对象编程:S3和R6
R语言面向对象编程:S3和R6 2017-06-10 0 R语言面向对象编程:S3和R6 一.基于S3的面向对象编程基于S3的面向对象编程是一种基于泛型函数(generic function)的实现 ...
mcem r语言代码_R语言关于if语句missing value where TRUEFALSE needed的错误
apriori df #找出频繁候选1项集 for(i in 1:length(data)) for(j in 1:length(data[[i]])){ if(any(df$candidate == ...
python 梯度提升树_机器学习：梯度提升算法｜python与r语言代码实现
梯度提升算法 10.1 GBM GBM(梯度提升机)是一种在处理大量数据以进行高预测的预测时使用的Boosting算法.Boosting实际上是一个学习算法的集合,它结合了几个基本估计量的预测,以便比 ...
excel计算二元线性回归_多重共线性及解决方法（附R语言代码）
一.多重共线性及其危害我有一组自变量: 它们满足下列关系: 那么我们称这四个变量间存在多重共线性. 这就意味着,一个变量可以被另外几个变量来解释,这就会带来两个后果 1.尽管系数估计的无偏性没有收到 ...
R语言通过loess去除某个变量对数据的影响
当我们想研究不同sample的某个变量A之间的差异时,往往会因为其它一些变量B对该变量的固有影响,而影响不同sample变量A的比较,这个时候需要对sample变量A进行标准化之后才能进行比较.标 ...
linux r的数据是存在,R语言通过loess去除某个变量对数据的影响
当我们想研究不同sample的某个变量A之间的差异时,往往会因为其它一些变量B对该变量的固有影响,而影响不同sample变量A的比较,这个时候需要对sample变量A进行标准化之后才能进行比较.标准化 ...

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

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

最新文章

热门文章