简介
本篇博客主要介绍如何利用EM算法进行多维正态缺失数据的参数估计，并进行R代码的实现。这里主要是使用Sweep Operator来实现。首先感谢我的队友JB大哥、xiaojj舍友以及杰哥，大家的共同努力，才完成了这份作业。
关于Sweep Operator来做Multiple Regression，说的比较好的是下面这个网站，大家可以先进行学习（因为后面算法的核心就是这个）：
- The SWEEP Algorithm for Multiple Regression
E步
首先我们写出complete data likelihood，然后写出其充分统计量。接着我们再对充分统计量求条件期望（基于观测），具体公式如下：
M步
M步就更加简单了，直接将参数化为前面充分统计量的表示形式，然后有缺失的用其期望代替。（无缺失的期望就是本身）
下面结合R代码继续说明。
定义相关的函数
首先我们可以照着公式打sweep operator。注意这里是整个EM算法的核心，需要反复用到，所以必须需要尽可能的提速，这样才可在最后将我们的算法速度进行进一步提升。所以我们在这里使用了向量化的操作。
## Sweep operator
Swp <- function(A, k){
A[-k, -k] <- A[-k, -k] - tcrossprod(A[-k, k]) / A[k, k] # tcrossprod: A[-k, k] %*% t(A[-k, k])
A[-k, k] <- A[-k, k] / A[k, k]
A[k, -k] <- A[-k, k]
A[k, k] <- -1 / A[k, k]
return(A)
}下面我们记录下每个样本中，未缺失的变量，记为向量mark，然后迭代使用前面的算子。
## Repeated calculation
SwpEstR <- function(para, mark){
for (i in mark) {
para <- Swp(para, i + 1)
}
return(para)
}需要注意的是，由于我们整个算法最核心，用的次数最多的是前面两个函数，所以为了进一步提速，我们可以“开挂”——将前面的R代码修改为C++代码。其实由于Rcpp这一神包的存在，我们将其修改为C++代码，具体如下（下面代码写在Rcpp文件里面，source进来即可）。
#include <Rcpp.h>
using namespace Rcpp;
NumericMatrix Swp (NumericMatrix G, int k) {
int n = G.rows();
int m = G.cols();
NumericMatrix H(n, m);
for (int j = 0; j < n; j++) {
for (int l = 0; l < m; l++) {
if (j == k && l == k) {
H(j, l) = - 1 / G(j, l);
}
else if (j != k && l == k) {
H(j, k) = G(j, k) / G(k, k);
H(k, j) = H(j, k);
}
else if (j != k && l != k) {
H(j, l) = G(j, l) - G(j, k) * G(k, l) / G(k, k);
}
}
}
return(H);
}
// [[Rcpp::export]]
NumericMatrix SwpEstC (NumericMatrix para, NumericVector mark) {
int len_mark = mark.size();
for (int i = 0; i < len_mark; i++) {
para = Swp(para, mark(i));
}
return(para);
}下面我们开始构造一些开始用于插补的函数。
这里首先在Impute中，我们添加了一个选择——是基于C++还是R的，前者速度会快一倍以上（这个后面会绘图进行展示） 。
接着就到了这个函数非常精华的部分了。由于需要对样本进行遍历，我们这里使用了for循环。但是为了提速，我们使用向量化的操作，只对缺失的样本进行了遍历（当p很小的时候，将缺失模式进行分类，然后再进行循环可以进一步提速，但对于p很大的情形就没有意义了，甚至加上缺失模式进行分类的判断可能还会影响速度。
这里为了偷懒只写了对缺失样本进行遍历的方法，但我们其实可以根据n与p的关系，写一个判别，然后根据不同的情况，调用不同的循环方式）。代码里面的这两行其实是最核心的部分：
Y_na[i, -mark] <- A[-Mark, Mark] %*% c(1, Y_na[i, mark])
C[-mark, -mark] <- C[-mark, -mark] + A[-Mark, -Mark]我们巧妙的运用了向量化的编程技巧，一个是为了插补缺失的Y，用于估计均值。第二个是为了估计后面的协方差阵。（当两个变量都缺失时，我们就需要估计两者乘积的期望，可以化为方差 + 期望的乘积，而C就是我们期望的乘积）
这里用的非常巧妙，不需要像传统的用法，将未缺失的排列到sweep operator处理过后矩阵的左上角，而是直接使用index进行选取，最后插补Y，并且计算C。
## Imputation of each step of EM
Impute <- function(Y, para, base_C){
p <- ncol(Y)
C <- matrix(0, p, p)
if (base_C == F) {
SwpEst <- SwpEstR
} else {
SwpEst <- SwpEstC
}
ind_na <- (Y)
Y_na_ind <- which(apply(ind_na, 1, any))
Y_na <- Y[Y_na_ind, ]
for (i in 1:nrow(Y_na)) {
mark <- which(!ind_na[Y_na_ind[i], ])
A <- SwpEst(para, mark)
Mark <- c(1, mark + 1)
Y_na[i, -mark] <- A[-Mark, Mark] %*% c(1, Y_na[i, mark])
C[-mark, -mark] <- C[-mark, -mark] + A[-Mark, -Mark]
}
Y[Y_na_ind, ] <- Y_na
return(list(Y = Y, C = C))
}有了上面的Impute，下面我们就可以很轻松的进行一步EM。核心是下面的代码，这就是我们的M步：
mu_t <- colMeans(Y)
sigma_t <- (tcrossprod(t(Y) - mu_t) + C) / n然后再整理成我们的参数矩阵，便可完成一步EM的迭代了。
## one-step EM
EM1 <- function(para, Y, base_C) {
n <- nrow(Y)
Y_imp <- Impute(Y, para, base_C)
Y <- Y_imp$Y
C <- Y_imp$C
mu_t <- colMeans(Y)
sigma_t <- (tcrossprod(t(Y) - mu_t) + C) / n
para[-1, 1] <- para[1, -1] <- mu_t
para[-1, -1] <- sigma_t
return(para)
}之后我们就可以设置收敛条件（max((para_2 - para_1) ^ 2) >= eps），进行多步EM算法的更新迭代。
这里加了一个判断：如果你给了初值，我们就是你给定的初始条件。但是如果没有给定初值，我们就使用均值向量为0，方差为单位矩阵。
重复多次EM1即可。
## EM
EM <- function(Y, eps = 1e-3, base_C = FALSE, ini_mu = NULL, ini_Sigma = NULL) {
p <- ncol(Y)
if (is.null(ini_mu) && is.null(ini_Sigma)) {
para_1 <- diag(p + 1)
para_1[1, 1] <- -1
} else {
para_1 <- matrix(-1, p + 1, p + 1)
para_1[-1, 1] <- para_1[1, -1] <- ini_mu
para_1[-1, -1] <- ini_Sigma
}
t <- Sys.time()
while(1) {
para_2 <- EM1(para_1, Y, base_C = base_C)
if (max((para_2 - para_1) ^ 2) < eps) break
para_1 <- para_2
}
delta_t <- Sys.time() - t
return(list(mu_ets = para_2[-1, 1], Sigma_ets = para_2[-1, -1], delta_t = delta_t))
}到这里EM算法就完成了，同时我们可以进行均值向量与协方差阵的估计了。但其实对于一个完整的模拟而言，到这还远远没有完成，我们还需要对估计出来的参数的方差进行估计（看看我们估计的是否准确），同时还需要生成数据进行检测我们的算法究竟估计的准不准。
