R中网格上有效的Montecarlo模拟

Question

我正在运行Montecarlo模拟多项逻辑。因此，我有一个函数来生成数据并估计模型。此外，我希望在一个值网格上生成不同的数据集。特别是，改变个体的数量(n.indiv)和每个个体的答案数(n.choices)。

到目前为止，我已经成功地解决了这个问题，但是在某个时候，我在网格搜索个人数量(n.indiv_list)和每个人的答案数(n.choices_list)的可能值时，产生了一个嵌套的for -循环结构。最后，我非常担心在可能的值组合上运行双循环结构的最后一段代码的使用效率。可能有一种矢量化的方法可以让我错过(也许不是吗？)

最后，这主要是一种风格问题，我设法找到了一个包含来自网格搜索组合的模型的倍数对象，这些模型包含了信息丰富的名称，但如果我能够将它们全部折叠到一个列表中，但在当前的结构中，我不知道该如何做，这也是很好的。提前谢谢你！

1)函数，该函数生成数据并估计模型.

library(dplyr)
library(VGAM)
library(mlogit)

#function that generates the data and estimates the model.
mlogit_sim_data <- function(...){
  
  # generating number of (n.alter) X (n.choices)
  df <- data.frame(id= rep(seq(1,n.choices ),n.alter ))
  
  # id per individual
  df <- df %>%
    group_by(id) %>%
    mutate(altern = sequence(n()))%>%
    arrange(id)
  
  #Repeated scheme for each individual + id_ind
  df <- cbind(df[rep(1:nrow(df), n.indiv), ], id_ind = rep(1:n.indiv, each = nrow(df)))
  
  ## creating attributes
  df<- df %>%
    mutate(
      x1=rlnorm(n.indiv*n.alter),
      x2=rlnorm(n.indiv*n.alter),
    )%>%
    group_by(altern) %>%
    mutate(
      id_choice = sequence(n()))%>%
    group_by(id_ind) %>%
    mutate(
      z1 = rpois(1,lambda = 25),
      z2 = rlnorm(1,meanlog = 5, sdlog  = 0.5),
      z3 = ifelse(runif(1, min = 0 , max = 1) > 0.5 , 1 , 0)
    )
  
  # Observed utility
  df$V1 <- with(df,  b1  * x1 +   b2 * x2 )
  
  #### Generate Response Variable ####
  fn_choice_generator <- function(V){
    U <- V + rgumbel(length(V), 0, 1)
    1L * (U == max(U))
  }
  
  # Using fn_choice_generator to generate 'choice' columns 
  df <-  df %>%
    group_by(id_choice) %>%
    mutate(across(starts_with("V"), 
                  fn_choice_generator, .names = "choice_{.col}")) %>% # generating choice(s)
    select(-starts_with("V")) %>% ##drop V variables.
    select(-c(id,id_ind))
  
  
  tryCatch(
    {
      model_result <- mlogit(choice_V1 ~ 0 +  x1 + x2 |1  ,
                                                  data = df,
                                                  idx = c("id_choice", "altern"))
      return(model_result)
    },
    error = function(e){
      return(NA)
    }
  )
  
}

2)网格搜索数据的可能组合

#List with the values that varies in the simulation
  #number of individuals
  n.indiv_list <- c(1, 15, 100, 500 ) 
  #number of choice situations
  n.choices_list <- c(1, 2, 4, 8, 10)  

# Values that remains constant across simulations 
  #set number of alternatives
  n.alter   <- 3    

## Real parameters
b1 <- 1
b2 <- 2

#Number of reps
nreps <- 10 
#Set seed
set.seed(777)

#iteration over different values in the simulation 
for(i in n.indiv_list) {
  for(j in n.choices_list) {
    n.indiv <- i
    n.choices <- j
    assign(paste0("m_ind_", i, "_choices_", j), lapply(X   = 1:nreps, FUN = mlogit_sim_data))
  }
}

Answer 1

您可以使用map2包的purrr函数矢量化：

library(tidyverse)

n.indiv_list <- c(1, 15, 100, 500 ) 
#number of choice situations
n.choices_list <- c(1, 2, 4, 8, 10)
l1 <- length(n.indiv_list)
l2 <- length(n.choices_list)
v1 <- rep(n.indiv_list, each = l2)
v2 <- rep(n.choices_list, l1)  #v1, v2 generate all pairs
> v1
 [1]   1   1   1   1   1  15  15  15  15  15 100 100 100 100 100 500 500 500 500 500
> v2
 [1]  1  2  4  8 10  1  2  4  8 10  1  2  4  8 10  1  2  4  8 10
    
result <- map2(v1, v2, function(v1, v2) assign(paste0("m_ind_", v1, "_choices_", v2), lapply(X   = 1:nreps, FUN = mlogit_sim_data)))

结果将是您的函数输出列表。