创建一个函数来表征重复模拟

Question

我想创建一个函数来帮助表征某些模拟的结果。出于此 post 的目的，让模拟函数为：

example_sim <- function(time=100, npops=5){
  result <- data.frame(matrix(NA, nrow = time, ncol = npops))
  colnames(result) <- LETTERS[1:npops]
  for(i in 1:npops){
    sim <- sample.int(time, time)
    result[,i] <- sim
    result[,i] <- result[,i]*i
  }
  return(result)
}

这将根据人口数量 (npops) 和模拟时间创建一个长度和宽度不同的数据框。

我想创建一个函数，它使用此类模拟的输出并表征 n 次模拟 (nsim) 中每个总体的均值和方差。

到目前为止，我已经设法使用以下代码使其适用于两个人群：

library("matrixStats")
library("reshape2")

ensembles <- function(nsims=10, time = 100, npops = 2){
  result_N.A <- data.frame(matrix(NA, nrow = time, ncol = nsims))
  result_N.B <- data.frame(matrix(NA, nrow = time, ncol = nsims))
  for( i in 1:(nsims)){
    simulation_with_2pops <- example_sim(time=100,npops=2)
    result_N.A[,i] <- simulation_with_2pops[,1]
    result_N.B[,i] <- simulation_with_2pops[,2]
  }
  output <- simulation_with_2pops
  for( j in 1:params$ntime){
    output$meanA[j] <- rowMeans(result_N.A[j,])
  }
  for( j in 1:params$ntime){
    output$meanB[j] <- rowMeans(result_N.B[j,])
  }
  for( j in 1:params$ntime){
    output$varA[j] <- rowVars(as.matrix(result_N.A[j,]))
  }
  for( j in 1:params$ntime){
    output$varB[j] <- rowVars(as.matrix(result_N.B[j,]))
  }
  return(output)
} 
ensembles_output<- ensembles(nsims = 10)
ensembles_output

要为任意数量的人群完全实现该功能，我需要创建另一个 for 循环，我在其中创建和更新 result_N.A 对象。（大概叫做 result[i] 之类的东西。）我也考虑过创建一个 3 维对象（时间、npops、nsims）并取其一部分来计算均值和方差，但我还没有取得太大的成功。我没有为这条路线结婚，对其他建议持开放态度。

最终我想创建一个代码，其中还通过突出显示参数中的两个总体来计算协方差和相关性。（例如人口 A 和人口 E）。如果您对实施有任何想法，我将不胜感激。

感谢您考虑这个问题。

Answer 1

我认为在这种情况下使用多维数组是一个非常好的主意。

首先，您可以使用 mapply() 更便宜地获得 example_sim() 的模拟。这里有一个 time=10 和 npops=3 的例子。使用相同的 set.seed(42) 和参数并自行检查。

我在这里使用的参数要小得多，这样您就可以轻松地在脑海中检查结果。

set.seed(42)
sim <- replicate(nsims, mapply(\(time, i) sample.int(time, time)*i, 10, 1:3))

sim
# , , 1
# 
#       [,1] [,2] [,3]
#  [1,]    1   16   27
#  [2,]    5   14   30
#  [3,]   10    8    9
#  [4,]    8    2   12
#  [5,]    2   10   15
#  [6,]    4   20   18
#  [7,]    6    4    3
#  [8,]    9   12    6
#  [9,]    7   18   24
# [10,]    3    6   21
# 
# , , 2
# 
#       [,1] [,2] [,3]
#  [1,]    3   10   18
#  [2,]    1    8    6
#  [3,]    2    4   12
#  [4,]    6   16    9
#  [5,]   10    6   30
#  [6,]    8    2   15
#  [7,]    4   20   27
#  [8,]    5   14   21
#  [9,]    7   12   24
# [10,]    9   18    3
# 
# , , 3
# 
#       [,1] [,2] [,3]
#  [1,]   10    8   18
#  [2,]    8   18    6
#  [3,]    5    6   27
#  [4,]    1   16    3
#  [5,]    7   10   24
#  [6,]    4   12   15
#  [7,]    6   20   30
#  [8,]    2    4    9
#  [9,]    9    2   12
# [10,]    3   14   21

接下来，我相信您想要收集 row-wise 每个人口列 A、B、C... 的统计数据。这里你基本上想要apply(., MARGINS=1:2, FUN)。就平均而言，存在 rowMeans(., dims=2L)，速度更快。

rowMeans(sim, dims=2L)
#           [,1]      [,2] [,3]
#  [1,] 4.666667 11.333333   21
#  [2,] 4.666667 13.333333   14
#  [3,] 5.666667  6.000000   16
#  [4,] 5.000000 11.333333    8
#  [5,] 6.333333  8.666667   23
#  [6,] 5.333333 11.333333   16
#  [7,] 5.333333 14.666667   20
#  [8,] 5.333333 10.000000   12
#  [9,] 7.666667 10.666667   20
# [10,] 5.000000 12.666667   15

apply(sim, 1:2, var)
#            [,1]      [,2] [,3]
#  [1,] 22.333333 17.333333   27
#  [2,] 12.333333 25.333333  192
#  [3,] 16.333333  4.000000   93
#  [4,] 13.000000 65.333333   21
#  [5,] 16.333333  5.333333   57
#  [6,]  5.333333 81.333333    3
#  [7,]  1.333333 85.333333  219
#  [8,] 12.333333 28.000000   63
#  [9,]  1.333333 65.333333   48
# [10,] 12.000000 37.333333  108

不过我不确定为什么要使用 simulation_with_2pops 作为最终输出，因为它是 for (i in 1:nsims) 循环的最后一次迭代的结果。无论如何，希望这对你有进一步的帮助。

注意：R >= 4.1 使用。

创建一个函数来表征重复模拟

Creating a function to characterise repeated simulations

simulation

statistics

r

matrix

multidimensional-array