在 R 中有 运行 Monte Carlo 模拟的有效方法吗?

Is there an efficient way of running Monte Carlo simulations in R?

我有一个 table,其中包含 'events' 的列表。每个事件都有每年的平均发生率和两个形状参数,用于从 Beta 分布生成随机值。

每年的总平均事件数是 table 中 Rate 列的总和,并假设这是泊松分布。

我想生成一个 Monte Carlo 模拟,这样对于每个模拟,我都会从泊松分布中生成随机数量的事件,然后对于每个随机生成的事件,我想从原始 FreqTable 中随机抽取一行与其中的速率成比例,并合并该行的形状参数。然后我想使用这些形状参数从 Beta 分布生成随机损失。

我在下面提供了一些可重现的代码,展示了我使用 for 循环对此进行的初步尝试。我陷入了需要通过在 FreqTable$CumulativeRate 中找到与 'rand' 值最接近的匹配来合并 FreqTable 中的形状参数的点。我试过使用 purr 库中的 detect_index,但我无法让它在循环中工作。

我打算在代码运行时生成至少 100k 次模拟,我知道 for 循环可能不是最有效的方法,但(由于我的编码能力有限)这就是我的全部能想到!

FreqTable <- data.frame(Event=c(1:10),
                        Rate=seq(0.1, 1, length.out = 10),
                        shape1=c(3:12),
                        shape2=c(5:14))


annual_freq <- sum(FreqTable$Rate)
FreqTable$CumulativeRate <- cumsum(FreqTable$Rate)/annual_freq

sims = 10

LossesList <- list()

for (i in 1:sims){
  num_losses <- rpois(1, annual_freq)
  sim_number <- rep(i, num_losses)
  rand <- runif(num_losses, 0, 1)  
  losses_list <- cbind(sim_number, rand)

  # PSEUDO CODE
  #look up shape parameters from FreqTable by using the nearest match on FreqTable$CumulativeRate from the rand value generated above
  #then generate a random variable from the Beta distribution for each row in the losses_list using the looked up shape parameters for that row

  LossesList[[i]] <- losses_list
}

SimulatedResults <- unlist(losses_list)

将 Roland 的建议形式化,您不需要使用 for 循环 - 依靠矢量化函数将提高性能。

让我知道这是否适合您:

FreqTable <- data.frame(Event=c(1:10),
                        Rate=seq(0.1, 1, length.out = 10),
                        shape1=c(3:12),
                        shape2=c(5:14))


annual_freq <- sum(FreqTable$Rate)
FreqTable$CumulativeRate <- cumsum(FreqTable$Rate)/annual_freq

sims = 10


sim.num <- 1:sims
num_losses <- rpois(sims, annual_freq)
rand.vals <- sapply(num_losses, runif) #list of num_losses[i] random values
sapply(rand.vals, cat)
table.ids2 <- lapply(rand.vals, function(x) {sapply(x, function(el){which.min(abs(FreqTable$CumulativeRate-el))}) })

sim.param.df <- data.frame()
for (i in 1:length(table.ids2)){

  closest.freq.table.event <- table.ids2[[i]]
  beta.shape1 <- sapply(table.ids2[[i]], FUN = function(x){FreqTable[FreqTable$Event == x, "shape1"]})
  beta.shape2 <- sapply(table.ids2[[i]], FUN = function(x){FreqTable[FreqTable$Event == x, "shape2"]})
  sim.id.vec <- rep(i, length(beta.shape1))

  temp.df <- data.frame(
    sim.id = sim.id.vec,
    closest.event = closest.freq.table.event,
    shape1 = beta.shape1,
    shape2 = beta.shape2
  )
  sim.param.df <- rbind.data.frame(sim.param.df, temp.df, make.row.names = FALSE)

}