R中的模拟最大似然，MaxLik

Question

我试图通过 R 中的 MaxLik 包通过模拟最大似然来估计模型。不幸的是，随着数据量的增加，我遇到了严重的性能问题。任何人都可以就以下内容提出建议：

有没有办法加快我的代码（它已经矢量化了，所以我对如何进一步改进它一无所知）？ 有没有办法通过 Rcpp 来实现优化过程以加快速度？ 有没有更聪明的方法来使用自定义似然函数来实现模拟最大似然？

我已经在 AWS 实例上尝试过 doParallel，但这并没有显着加快进程。

我创建了一个可重现的示例并对最重要的部分进行了注释：

#create data:
#Binary DV (y), 10 IDV (V3 - V12), 50 groups (g), with 100 sequential observations each (id)
set.seed(123)
n <- 5000
p <- 10
x <- matrix(rnorm(n * p), n)
g <- rep(seq(1:(n/100)),each=100)
id <- rep(seq(1:(n/max(g))),max(g))
beta <- runif(p)
xb <- c(x %*% beta)
p <- exp(xb) / (1 + exp(xb))
y <- rbinom(n, 1, p)
data <- as.data.table(cbind(id,y,x,g))

#Find starting values for MaxLik via regular glm
standard <-
  glm(
    y  ~ 
      V3 +
      V4 +
      V5 +
      V6 +
      V7 +
      V8 +
      V9 +
      V10 +
      V11 +
      V12,
    data = data,
    family = binomial(link = "logit")
  )
summary(standard)

#set starting values for MaxLik
b <- c(standard$coefficients,sd_V3=0.5,sd_V4=0.5)

#draw 50 x # of groups random values from a normal distribution
draws <- 50
#for each g in the data, 50 randomvalues are drawn
rands <- as.data.table(cbind(g=rep(g,each=draws),matrix(rnorm(length(g)*draws,0,1),length(g)*draws,2)))
colnames(rands) <- c("g","SD_V3","SD_V4")
#merge random draws to each group, so every observation is repeated x # of draws
data <- merge(data,rands,by="g",all=T,allow.cartesian=T)

#the likelihood function (for variables V3 and V4, a mean [b3] & b[4] and a SD b[12] & b[14] is estimated
loglik1 <- function(b){

#I want the standard deviations to vary only across groups (g), but all other parameters to vary across all observations, which is why I am taking the mean across g and id (remember, every observation is a cartesian product with the random draws per group)

  ll <- data[,.(gll=mean(((1/(1+exp(-(b[1]+
                                  (b[2]+b[12]*SD_V3)*V3 + 
                                  (b[3]+b[13]*SD_V4)*V4 + 
                                  (b[4])*V5 + 
                                  (b[5])*V6 + 
                                  (b[6])*V7 + 
                                  (b[7])*V8 + 
                                  (b[8])*V9 + 
                                  (b[9])*V10 + 
                                  (b[10])*V11 + 
                                  (b[11])*V12))))^y)*
                     (1-(1/(1+exp(-(b[1]+
                                    (b[2])*V3 + 
                                    (b[3])*V4 + 
                                    (b[4])*V5 + 
                                    (b[5])*V6 + 
                                    (b[6])*V7 + 
                                    (b[7])*V8 + 
                                    (b[8])*V9 + 
                                    (b[9])*V10 + 
                                    (b[10])*V11 + 
                                    (b[11])*V12)))))^(1-y))),by=.(g,id)]
  return(log(ll[,gll]))
}

co <- maxLik::maxControl(gradtol=1e-04,printLevel=2)
maxlik <- maxLik::maxLik(loglik1,start=b,method="bfgs",control=co)
summary(maxlik)

非常感谢您的建议

Answer 1

通过将 loglik1 ，我能够显着减少优化时间（几小时到几分钟）

return(data[,.(g,id,y,logit=1/(1+exp(-(b[1]+
                                      (b[2]+b[12]*SD_V3)*V3 + 
                                      (b[3]+b[13]*SD_V4)*V4 + 
                                      (b[4])*V5 + 
                                      (b[5])*V6 + 
                                      (b[6])*V7 + 
                                      (b[7])*V8 + 
                                      (b[8])*V9 + 
                                      (b[9])*V10 + 
                                      (b[10])*V11 + 
                                      (b[11])*V12))))][,mean(y*log(logit)+(1-y)*log(1)-logit),by=.(g,id)][,sum(V1)])

然而，这只是部分解决了问题，因为随着数据量的增加，估计时间再次增加:(

除非有人有优雅的解决方案，否则我可能不得不处理这个问题？

编辑：过一段时间再了解这一点，以防将来有人遇到问题......原因，脚本需要很长时间，在于包 MaxLik 和推导 Hessian 矩阵的计算时间。如果你不需要，你可以告诉 MaxLik 不要计算它。因为我确实需要它，所以我决定通过 Rcpp 来计算它。