Есть ли эффективный способ запуска симуляции Монте-Карло в R? - PullRequest
Новая
       71

Есть ли эффективный способ запуска симуляции Монте-Карло в R?

0 голосов
/ 28 января 2020

У меня есть таблица, которая содержит список «событий». Каждое событие имеет среднюю годовую ставку вместе с двумя параметрами формы для генерации случайного значения из распределения бета-версии.

Общее среднее число событий в год представляет собой сумму столбца Rate в таблице и составляет Предполагается, что это распределение Пуассона.

Я хочу сгенерировать симуляцию Монте-Карло таким образом, чтобы для каждого моделирования я генерировал случайное число событий из распределения Пуассона, а затем для каждого случайно сгенерированного события я хочу нарисовать случайное Строка из исходного FreqTable пропорционально ставкам в нем и объединить параметры формы из этой строки. Затем я хочу сгенерировать случайную потерю из бета-распределения, используя эти параметры формы.

Ниже приведен некоторый воспроизводимый код, демонстрирующий мою первоначальную попытку сделать это с использованием a для l oop. Я застрял в точке, где мне нужно объединить параметры формы из FreqTable, найдя наиболее близкое совпадение со значением 'rand' в FreqTable $ CumulativeRate. Я пытался использовать dete_index из библиотеки purr, но я не могу заставить это работать в al oop.

. Я намереваюсь сгенерировать как минимум 100k симуляций, когда код работает, и я Я знаю, что для l oop, возможно, не самый эффективный способ сделать это, но (с моей ограниченной способностью кодирования) это все, что я мог придумать! 

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)

1 Ответ

1 голос
/ 28 января 2020

Формализуя предложение Роланда, вам не нужно использовать петли 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)

}
...