arrays - 从 R 中的另一个 3D 数组填充 3D 数组的最快方法

Question

我正在使用下面的代码从另一个 3D 数组填充 3D 数组。我用过 sapply函数在每个人(第三维)应用代码行，如 Efficient way to fill a 3D array .
这是我的代码。

ind <- 1000
    individuals <- as.character(seq(1, ind, by = 1))
    maxCol <- 7
    col <- 4
    line <- 0
    a <- 0
    b <- 0
    c <- 0

    col_array <- c("year","time", "ID", "age", as.vector(outer(c(paste(seq(0, 1, by = 1), "year", sep="_"), paste(seq(2, maxCol, by = 1), "years", sep="_")), c("S_F", "I_F", "R_F"), paste, sep="_")))
    array1 <- array(sample(1:100, length(col_array), replace = T), dim=c(2, length(col_array), ind), dimnames=list(NULL, col_array, individuals)) ## 3rd dimension = individual ID
    ## print(array1)

    col_array <- c("year","time", "ID", "age", as.vector(outer(c(paste(seq(0, 1, by = 1), "year", sep="_"), paste(seq(2, maxCol, by = 1), "years", sep="_")), c("S_M", "I_M", "R_M"), paste, sep="_")))
    array2 <- array(NA, dim=c(2, length(col_array), ind), dimnames=list(NULL, col_array, individuals)) ## 3rd dimension = individual ID
    ## print(array2)

    tic("array2")
    array2 <- sapply(individuals, function(i){

      ## Fill the first columns
      array2[line + 1, c("year", "time", "ID", "age"), i] <- c(a, b, i, c)

      ## Define column indexes for individuals S
      col_start_S_F <- which(colnames(array1[,,i])=="0_year_S_F")
      col_end_S_F <- which(colnames(array1[,,i])==paste(maxCol,"years_S_F", sep="_"))
      col_start_S_M <- which(colnames(array2[,,i])=="0_year_S_M")
      col_end_S_M <- which(colnames(array2[,,i])==paste(maxCol,"years_S_M", sep="_"))

      ## Fill the columns for individuals S
      p_S_M <- sapply(0:maxCol, function(x){pnorm(x, 4, 1)})
      array2[line + 1, col_start_S_M:col_end_S_M, i] <- round(as.numeric(as.vector(array1[line + 1, col_start_S_F:col_end_S_F, i]))*p_S_M)

      ## Define column indexes for individuals I
      col_start_I_F <- which(colnames(array1[,,i])=="0_year_I_F")
      col_end_I_F <- which(colnames(array1[,,i])==paste(maxCol,"years_I_F", sep="_"))
      col_start_I_M <- which(colnames(array2[,,i])=="0_year_I_M")
      col_end_I_M <- which(colnames(array2[,,i])==paste(maxCol,"years_I_M", sep="_"))

      ## Fill the columns for individuals I
      p_I_M <- sapply(0:maxCol, function(x){pnorm(x, 2, 1)})
      array2[line + 1, col_start_I_M:col_end_I_M, i] <- round(as.numeric(as.vector(array1[line + 1, col_start_I_F:col_end_I_F, i]))*p_I_M)

      ## Define column indexes for individuals R
      col_start_R_M <- which(colnames(array2[,,i])=="0_year_R_M")
      col_end_R_M <- which(colnames(array2[,,i])==paste(maxCol,"years_R_M", sep="_"))

      ## Fill the columns for individuals R
      array2[line + 1, col_start_R_M:col_end_R_M, i] <- as.numeric(as.vector(array2[line + 1, col_start_S_M:col_end_S_M, i])) + 
        as.numeric(as.vector(array2[line + 1, col_start_I_M:col_end_I_M, i]))

      return(array2[,,i])
      ## print(array2[,,i])

    }, simplify = "array") 
    ## print(array2)
    toc()

有没有办法提高我的代码的性能/速度(即 < 1 秒)？第 3 维有 500000 个观测值。有什么建议？

Answer 1

TL;DR:这是一个 tidyverse 解决方案，它将样本数组转换为数据帧并应用请求的更改。 编辑:我添加了步骤 1+2 将原始帖子的示例数据转换为我在步骤 3 中使用的格式。步骤 3 中的实际计算非常快(<0.1 秒)，但瓶颈是步骤 2，它需要50 万行需要 10 秒。

步骤 0:为 50 万人创建样本数据

ind <- 500000
individuals <- as.character(seq(1, ind, by = 1))
maxCol <- 7
col <- 4
line <- 0
a <- 0
b <- 0
c <- 0

col_array <- c("year","time", "ID", "age", as.vector(outer(c(paste(seq(0, 1, by = 1), "year", sep="_"), paste(seq(2, maxCol, by = 1), "years", sep="_")), c("S_F", "I_F", "R_F"), paste, sep="_")))
array1 <- array(sample(1:100, length(col_array), replace = T), dim=c(2, length(col_array), ind), dimnames=list(NULL, col_array, individuals)) ## 3rd dimension = individual ID

dim(array1)
# [1]      2     28 500000    # Two rows x 28 measures x 500k individuals

步骤 1:子集数组并转换为数据帧。

library(tidyverse)
# OP only uses first line of array1. If other rows needed, replace with "array1 %>%"
#   and adjust renaming below to account for different Var1.
array1_dt <- array1[1,,] %>% 
  as.data.frame.table(stringsAsFactors = FALSE)

第 2 步:将统计数据分成不同的列，每个年份占一行。这是最慢的一步(尤其是 spread 行)，1000 个人需要 0.05 秒，500k 需要 10 秒。如果需要，我希望 data.table 解决方案可以使它更快。

array1_dt_reshape <- array1_dt %>%
  rename(stat = Var1, ID = Var2) %>%
  filter(!stat %in% c("year", "time", "ID", "age")) %>%
  mutate(year = stat %>% str_sub(end = 1),
         col  = stat %>% str_sub(start = -3)) %>%
  select(-stat) %>%
  spread(col, Freq) %>%
  arrange(ID)

第 3 步:应用请求的转换。此函数使用两组参数计算分布，并使用这些参数来缩放输入表的列。 500k 个人需要 0.03 秒。

array_transform <- function(input_data = array1_dt_reshape, 
                           max_yr = 7, S_M_mean = 4, I_M_mean = 2) {
  tictoc::tic()
  # First calculate the distribution function values to apply to all individuals, 
  #   depending on year.
  p_S_M_vals <- sapply(0:max_yr, function(x){pnorm(x, S_M_mean, 1)})
  p_I_M_vals <- sapply(0:max_yr, function(x){pnorm(x, I_M_mean, 1)})

  # For each year, scale S_M + I_M by the respective distribution functions.
  #   This solution relies on the fact that each ID has 8 rows every time,
  #   so we can recycle the 8 values in the distribution functions.
  output <- input_data %>% 
    # group_by(ID) %>%  <-- Not needed
    mutate(S_M = S_F * p_S_M_vals,
           I_M = I_F * p_I_M_vals,
           R_M = S_M + I_M)  # %>% ungroup  <-- Not needed
  tictoc::toc()
  return(output)
}


array1_output <- array_transform(array1_dt_reshape)

结果

head(array1_output)
   ID year I_F R_F S_F          S_M        I_M         R_M
1   1    0  16  76  23 7.284386e-04  0.3640021   0.3647305
2   1    1  46  96  80 1.079918e-01  7.2981417   7.4061335
3   1    2  27  57  76 1.729010e+00 13.5000000  15.2290100
4   1    3  42  64  96 1.523090e+01 35.3364793  50.5673837
5   1    4  74  44  57 2.850000e+01 72.3164902 100.8164902
6   1    5  89  90  64 5.384606e+01 88.8798591 142.7259228
7   1    6  23  16  44 4.299899e+01 22.9992716  65.9982658
8   1    7  80  46  90 8.987851e+01 79.9999771 169.8784862
9   2    0  16  76  23 7.284386e-04  0.3640021   0.3647305
10  2    1  46  96  80 1.079918e-01  7.2981417   7.406133