在 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)
}
我有一个 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)
}