利用EM算法进行多维正态缺失数据的参数估计（使用Sweep Operator）——2. 进行估计参数的方差估计

本篇文章接着上一篇：利用EM算法进行多维正态缺失数据的参数估计（使用Sweep Operator）——1. 进行均值向量与协方差阵的估计

在利用EM算法进行多维正态缺失数据的参数估计之后，我们如何来看我们估计参数的准确性呢（也就是为这些参数再估计一个方差）？

这里就不介绍boostrap与jackknife方法来进行估计。我们介绍另一种：使用Fisher信息阵来进行估计的方法：其实核心就一个，对估计出来的Fisher信息阵求一个inverse就可以了，具体的做法如下：

我们直接按照上面文中所示的公式进行输入即可，下面是对应的R代码。

---

计算估计参数的方差估计

## the expected information matrix
EstJ <- function(Y, p, result_EM) {
  
  n <- nrow(Y)
  phi <- array(0, c(p, p, n)) # 构造三维数组
  C_l <- result_EM$Sigma_ets
  
  for(i in 1:n) {
    mark <- which(!is.na(Y[i, ]))
    phi[mark, mark, i] <- solve(C_l[mark, mark])
  }
  
  J_mu <- apply(phi, c(1, 2), sum)
  
  # J_Sigma
  n_para <- p * (p + 1) / 2
  J_Sigma <- matrix(0, n_para, n_para)
  
  trans_array <- expand.grid(1:p, 1:p)
  ind <- trans_array[, 1] <= trans_array[, 2]
  trans_array <- trans_array[ind, ]
  
  CalJSigma <- function(x, y) {
  ## 这个对应关系构造的比较巧妙
    l <- trans_array[x, 1]
    m <- trans_array[x, 2]
    r <- trans_array[y, 1]
    s <- trans_array[y, 2]
    return((2 - (l == m)) * (2 - (r == s)) / 4 * sum(phi[l, r, ] * phi[m, s, ] + phi[l, s, ] * phi[m, r, ]))
  }
  
  for(x in 1:n_para) {
    for(y in 1:x) {
      J_Sigma[x, y] <- CalJSigma(x, y)
    }
  }
  J_Sigma[upper.tri(J_Sigma)] <- t(J_Sigma)[upper.tri(J_Sigma)]
  
  # J_Sigma <- outer(1:n_para, 1:n_para, CalJSigma)
  
  J <- matrix(0, p + n_para, p + n_para)
  J[1:p, 1:p] <- J_mu
  J[-(1:p), -(1:p)] <- J_Sigma
  
  return(J)
}

---

生成模拟数据相关函数（模拟数据、进行模拟、置信度）

## Generating data
Gendata <- function(mu, Sigma_sqrt, n, p, rate){
  
  Y <- t(mu + t(Sigma_sqrt) %*% matrix(rnorm(p * n), p, n))
  
  miss_mark <- (matrix(runif(p * n), n, p) < rate) ## MCAR，控制缺失率
  Y[miss_mark] <- NA
  
  ## Exclude the rows that are all missing
  Y <- Y[!apply(miss_mark, 1, all), ]
  
  return(Y)
}

## Simulation function
MySim <- function(n, mu, Sigma_sqrt, base_C = FALSE, cal_J = FALSE, 
                  ini_mu = NULL, ini_Sigma = NULL, rate = 0.1) {
  p <- length(mu)
  Y <- Gendata(mu, Sigma_sqrt, n, p, rate = rate)
  
  result_EM <- EM(Y, base_C = base_C)
  
  if (cal_J == TRUE) {
    J <- EstJ(Y, p, result_EM)
    return(list(result_EM = result_EM, J = J))
  
  } else {
    return(result_EM)
  }
}

## Confidence interval
CIFun <- function(result, mu = mu, Sigma = Sigma, conf = 0.95) {
  
  para_vec <- c(mu, Sigma[upper.tri(Sigma, diag = T)])
  quan <- qnorm((1 - conf) / 2, lower.tail = F)
  
  mu_vec <- result$result_EM$mu_ets
  sigma_vec <- result$result_EM$Sigma_ets
  est_para_vec <- c(mu_vec, sigma_vec[upper.tri(sigma_vec, diag = T)])
  
  sd_vec <- sqrt(diag(solve(result$J)))
  
  ci_l <- est_para_vec - quan * sd_vec
  ci_r <- est_para_vec + quan * sd_vec
  
  return(para_vec < ci_r & para_vec > ci_l)
}

由于对称正定阵没有那么好生成，这里比较技巧性的是我们通过生成Sigma_sqrt，再自己乘以自己，继而生成协方差阵。大部分的代码都是比较自然的想法，关键就是尽可能使用向量化操作进行提速。

结合上一篇文章的函数，下面就可以开始我们的模拟过程。

---

模拟

首先比较基于R的代码与基于C++的时间比较。

## Simulation
n <- 1000
p <- 3
rate <- 0.1
nrep <- 100

mu <- runif(p, -5, 5)  ## Mean vector
Sigma_sqrt <- matrix(runif(p * p, -2, 2), p, p)
Sigma <- crossprod(Sigma_sqrt)  ## Covariance matrix

sim1 <- MySim(n, mu, Sigma_sqrt, base_C = F)
sim2 <- MySim(n, mu, Sigma_sqrt, base_C = T)

sim1$delta_t
sim2$delta_t

参数估计的大致比较：

mu
sim1$mu_ets
sim2$mu_ets

重复多次模拟，进行比较估计参数的方差估计与真实估计参数的方差之间的比较。

set.seed(1)
final_R <- lapply(1:nrep, MySim, n = n, mu = mu, 
                  Sigma_sqrt = Sigma_sqrt, base_C = F, cal_J = FALSE)
final_C <- lapply(1:nrep, MySim, n = n, mu = mu, 
                  Sigma_sqrt = Sigma_sqrt, base_C = T, cal_J = FALSE)

## Compare J
sim3 <- MySim(n, mu, Sigma_sqrt, base_C = T, cal_J = T)
J <- sim3$J
est_cov <- solve(J)

mu_all <- sapply(final_C, function(x) x$mu_ets)
sigma_all <- sapply(final_C, function(x) x$Sigma_ets[upper.tri(x$Sigma_ets, diag = T)])
  
diag(est_cov[1:p, 1:p])
diag(cov(t(mu_all)))

diag(est_cov[-(1:p), -(1:p)])
diag(cov(t(sigma_all)))

---

到这里，基本整个模拟就告一段落了。但若想要得到一个更好的展示，还需要进行一个非常好的可视化呈现，这个我们在下篇博客中介绍一些绚丽的展现方式，拭目以待！