Блог

Я пытаюсь оценить совокупные неопределенности, связанные с различными входными параметрами, используя метод Марковской цепи Монте-Карло в R. Другими словами, используя параметры неопределенности, указанные в документации по входным данным, я пытаюсь генерировать распределения для каждого из наборов данных, создавая 1000 случайных значений в пределах используемого распределения (нормальное распределение или усеченное нормальное распределение). Однако я не знаю, как это сделать с помощью функций purrr::map() быстрее и без истощения оперативной памяти. Набор данных содержит 2 строки и 80 столбцов. Вот упрощенный пример:

 library(tidyverse); library(truncnorm);library(data.table); library(dtplyr)
n <- 1000 # number of simulations
n_obs <- 10000 # number of observations. Does not work if e.g. 50000
 

Создайте фрейм данных.

 dt <- data.frame(
  var1 = runif(n_obs, 0, 100),
  var2_low = runif(n_obs, 0, 1),
  var2_mean = runif(n_obs, 0, 5),
  var2_up = runif(n_obs, 0, 10)
  )
 

Преобразуйте в ленивую таблицу данных, чтобы ускорить процесс

 dt1 <- dt %>% as.data.table() %>%
  lazy_dt(., immutable = FALSE)
 

Имитировать

 dt_sim <- dt1 %>% 
  mutate(mean_val =  rep(1, nrow(.)), # just row of 1 
         var1_rnorm = map(.x = mean_val,~rnorm(n, mean = .x, sd = 0.10)), # normal distribution with given sd
         sim_var1 = map2(.x = var1, .y = var1_rnorm, ~(.x*.y))) %>%  # multiply the data with simulated distribution
  
  # add truncated normal distribution for each row (var2)
  mutate(sim_var2 = pmap(.,~ rtruncnorm(n, 
                                        a = dt$var2_low, 
                                        b = dt$var2_up,
                                        mean =dt$var2_mean))) %>% 
 
  # multiply simulated variables sim_var1 and sim_var2
  mutate(sim_vars_multiplied = 
           pmap(list(x = .$sim_var1,
                     y = .$sim_var2),
                function(x,y) (x*y))) %>% 
  
  # derive coefficient of variation
  mutate(var_mean =map(.x = sim_vars_multiplied, ~ mean(.x, na.rm = TRUE)),
         var_sd = map(.x = sim_vars_multiplied, ~ sd(.x, na.rm = TRUE)),
         var_cv = unlist(var_sd) / unlist(var_mean)) %>% 
  
  # select only the variables needed
  dplyr::select(var_cv)
  

# collect the results
sim_results <- dt_sim %>% as.data.table() 
 

Это может помочь

 library(data.table)
#n <- 1000 # number of simulations
n_obs <- 10000 
dt <- data.table(
  var1 = runif(n_obs, 0, 100),
  var2_low = runif(n_obs, 0, 1),
  var2_mean = runif(n_obs, 0, 5),
  var2_up = runif(n_obs, 0, 10)
  )

dt[, mean_val := rep(1,.N)]
dt[, var1_rnorm := rnorm(.N, mean = mean_val, sd = 0.10)]
dt[, sim_var1 := var1 * var1_rnorm]

dt[, sim_var2 := truncnorm::rtruncnorm(.N, a = var2_low, b = var2_up, mean = var2_mean)]
dt[, sim_vars_multiplied := sim_var1 * sim_var2]

dt[, var_mean := mean(sim_vars_multiplied, na.rm=TRUE)]
dt[, var_sd := sd(sim_vars_multiplied, na.rm=TRUE)]
dt[, var_cv := var_sd / var_mean]

sim_results <- dt[,var_cv]