Нужна помощь в построении симуляции Монте-Карло и нахождении процентилей результата с помощью R

Question

У меня есть файл CSV, содержащий набор событий (около 40 элементов), все из которых могут происходить или не происходить, в зависимости от заданной вероятности. Столбцы: название события, размер урожая, вероятность.

Меня интересуют эти данные: общий размер доходности набора (сумма всех выходов набора) и, возможно, также общая сумма доходности за событие. Таким образом, поскольку событие может произойти или нет, и, следовательно, общий размер доходности набора может отличаться, мне нужно провести симуляцию по методу Монте-Карло для набора, имея пробу Бернулли над столбцом "Вероятность".

И, наконец, мне нужно рассчитать процентили по сумме доходности всего набора или конкретного события по всем итерациям моделирования Монте-Карло (сценариям).

У меня проблемы с записью...(я все еще изучаю R, я более привык к Java/C# и т. Д.)

Код, который я сделал в настоящее время:

#Generate sample data for a set of events that I want to simulate
eventcol <- c('Event1', 'Event2', 'Event3', 'Event4', 'Event5')
yieldcol <- c(350, 200, 100, 120, 540)
problcol <- c(0.5, 0.2, 0.9, 0.4, 0.7)
events <- data.frame(Name=eventcol, Yield=yieldcol, Probability=problcol)

#Forecast function
forecast <- function(events){
  count <- nrow(events)
  data <- data.frame(Id=seq(1, count))
  data$Name <- events$Name
  data$Yield <- events$Yield
  data$Exists <- rbinom(count,1,events$Probability)
  return(data)
}

#Create Monte Carlo simulation scenarios/realizations
scenarios <- replicate(4, forecast(events))
scenarios

Вывод следующий:

> scenarios
       [,1]      [,2]      [,3]      [,4]     
Id     Integer,5 Integer,5 Integer,5 Integer,5
Name   factor,5  factor,5  factor,5  factor,5 
Yield  Numeric,5 Numeric,5 Numeric,5 Numeric,5
Exists Numeric,5 Numeric,5 Numeric,5 Numeric,5

Но я понятия не имею, как суммировать доходность по событиям, которые действительно существуют (существует == 1) по сценарию, не говоря уже о том, чтобы найти процентиль (с функцией квантиля) по суммам. Как бы вы поступили с этим?

Что касается структуры данных, у меня есть несколько идей, но я не уверен..

1. Может быть, я должен транспонировать прогноз, а затем как-то перебирать сценарии MC один за другим и суммировать данные?

2. Может быть, я должен отфильтровать события из результатов, которые не существуют (Exists == 0). Но как и где мне это сделать?

Вероятно, было бы также более разумно, если бы результаты выглядели так (но у меня также нет идей, как этого добиться):

Scenario     Name     Yield
1            Event1   350
1            Event2   200
2            Event1   350
...

Пожалуйста, поделитесь своими мыслями!

Спасибо!

Answer 1

Да, теперь вопрос гораздо яснее!

scenarios вывод представляет собой набор списков. scenarios[3,] содержит "потенциальную доходность", scenarios[4,] содержит "существует".

Фактическая доходность для каждого сценария может быть определена следующим образом:

potential_yields = simplify2array(scenarios[3,])
exists           = simplify2array(scenarios[4,])
actual_yields    = colSums(yields*exists)

Определить и построить квантили:

yield_q  = quantile(actual_yields,probs=0:100/100)
plot(x = 0:100, y = yield_q)

Возможно, это то, что вы ищете.