Моделирование максимального правдоподобия в R, MaxLik

Question

Я пытаюсь оценить модель путем моделирования максимальной вероятности с помощью пакета MaxLik в R. К сожалению, с увеличением размера данных я сталкиваюсь с серьезными проблемами производительности. Может ли кто-нибудь посоветовать следующее:

Есть ли способ ускорить мой код (он уже векторизован, поэтому я не знаю, как его улучшить)? Есть ли способ реализовать процесс оптимизации через R cpp, чтобы ускорить его? Есть ли более разумный способ реализовать имитацию максимальной вероятности с помощью настраиваемой функции правдоподобия?

Я уже пробовал doParallel на экземпляре AWS, но это не сильно ускоряет процесс.

Я создал воспроизводимый пример и прокомментировал наиболее важные части:

#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 <- function (b) {...} на </p>

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)])

Однако это решает лишь частично проблема, поскольку с увеличением размера данных время оценки снова увеличивается: (

Мне, вероятно, придется иметь дело с этим, если у кого-то нет элегантного решения?