1

我需要对包含 31000 个观测值的数据集进行 MLE 估计，并使用对数似然函数自编码。数据集有 21 个参数需要估计。为了进行估计，我在 R. Hessian 中使用了带有 newton-raphson 方法的 maxLik 包，并且未将对数似然函数的梯度赋予函数本身。以下是如何读取数据以及如何计算对数似然函数。

[编辑@Patric，使用随机数据]

nObs <- 30000
nVeh   <- runif(nObs, 1, 2)
adult1 <- runif(nObs, 0, 1)
adult2 <- runif(nObs, 0, 1)
adult3 <- runif(nObs, 0, 1)
adult4 <- runif(nObs, 0, 1)

params <- runif(21, 1, 4)

## Matrix to keep loglikelihood records for each observation
LL <- matrix (0, ncol=1, nrow=nObs)
logLik <- function(params){
for (i in 1:nObs){
            V0 <- params[1]  
            V1 <- params[2] + params[6] * adult1[i] + params[10] * adult2[i] + params[14] * adult3[i] + params[18] * adult4[i]
            V2 <- params[3] + params[7] * adult1[i] + params[11] * adult2[i] + params[15] * adult3[i] + params[19] * adult4[i]
            V3 <- params[4] + params[8] * adult1[i] + params[12] * adult2[i] + params[16] * adult3[i] + params[20] * adult4[i]
            V4 <- params[5] + params[9] * adult1[i] + params[13] * adult2[i] + params[17] * adult3[i] + params[21] * adult4[i]

            vU <- matrix (c(V0, V1, V2, V3, V4), nrow = 1, ncol = 5)

            pChosen <- exp(vU[min(nVeh[i] + 1,5)]) / sum(exp(vU))
            LL[i] <- log(pChosen)
    }
    return(sum(LL))
}

system.time(logLik(params))
# user  system elapsed 
# 0.7     0.0     0.7

上面的代码运行很慢。MLE 估计需要几天时间。因此，我根据分层随机抽样对数据进行了抽样，并将其从 31000 个观测值减少到 5000 个观测值。

然而，它仍然运行缓慢。因此，我决定使用 doParallel 和 foreach 库的后端函数在对数似然函数中并行化 for 循环。这是我添加和更改的 for 循环：

cl <- makeCluster(8)
registerDoParallel(cl)
..... some code here to read data ......
logLik <- function(params){
    foreach (i=1:(nObs), .combine = "cbind") %do% {
       ...some more code content of log likelihood written above ...
    }
return(sum(LL))
}

我不熟悉并行计算，我阅读了很多关于并尝试在较小的函数示例中实现我学到的东西。假设 V0、V1、V2、V3 和 V4 都用常量定义，上述结构在相同的代码中工作，只有五个参数。5分钟左右比较快。然而，当我用 21 个参数运行上述模型以估计 5000 个观测值的采样数据时，它再次变得非常缓慢。

有没有其他方法可以让 maxLik 函数运行得更快？任何线索或提示，以及任何关于更快运行时间的阅读材料都值得赞赏。我使用搜索键“maxLik”、“parallel in R”等检查了网站，但找不到建议或问题。

PS：由于没有样本数据，该问题可能会被标记为负面，但由于机密性，我无法提供样本数据。

提前致谢

4

0 回答 0