Блог

Я делаю несколько симуляций с разными распределениями вероятностей и с разными параметрами. Пример: распределение — нормальное, среднее значение параметра и sd. Среднее значение исправлено (например, 0), а sd изменяется, например, с 1 на 5. Я написал некоторую функцию для решения этой проблемы, простая версия, которая даст вам представление, приведена ниже:

 test <- function(x, fixed, value, changed, seq){
    expr <- eval(parse(text=stri_flatten(c(fixed,"=",value))))
    for(i in seq_along(seq)){
        expr2 <- eval(parse(text=stri_flatten(c(changed,"=",seq[i]))))
        y <- dnorm(x, expr, expr2)
        sqrt(y*y   y*y)
    }
}
  

В моей функции анализ занимает половину времени моделирования (13 секунд из 25 секунд). Мне нужно это улучшить, есть идеи, как? Вы можете запустить функцию Rprof и проверить результат:

 {Rprof();test(-10:10,"mean",0,"sd",seq(1,5,length.out = 1001));Rprof(NULL);}
summaryRprof()
  

Rprof результаты моего моделирования (только часть)

 summaryRprof()
$by.self
                       self.time self.pct total.time total.pct
"parse"                    13.48    53.71      15.42     61.43
"ifelse"                    2.26     9.00      23.64     94.18
"structure"                 1.56     6.22       2.66     10.60
".External"                 0.80     3.19      17.90     71.31
"$<-"                       0.54     2.15       0.54      2.15
"options"                   0.52     2.07       0.54      2.15
"makeRestartList"           0.50     1.99       1.98      7.89
"c"                         0.46     1.83       0.46      1.83
"match"                     0.40     1.59       0.40      1.59
"doWithOneRestart"          0.38     1.51       1.60      6.37
"eval"                      0.32     1.27      15.84     63.11
"stopifnot"                 0.26     1.04       0.34      1.35
".Call"                     0.26     1.04       0.28      1.12
"phirsch"                   0.24     0.96      24.46     97.45
"t"                         0.24     0.96       0.26      1.04
"any"                       0.24     0.96       0.24      0.96
"pbeta"                     0.22     0.88      22.00     87.65
"list"                      0.20     0.80       0.20      0.80
"%in%"                      0.16     0.64       0.38      1.51
"^"                         0.14     0.56       0.14      0.56
"floor"                     0.14     0.56       0.14      0.56
  

Одна из идей ускорить это — использовать внутреннюю функцию разбора:

 text <- "a <- 5"
internalParse <- function(x) .Internal(parse(stdin(), NULL, x, NULL, NULL, "unknown"))
microbenchmark(parse(text=text), internalParse(text))
## Unit: microseconds
##                 expr     min       lq  median       uq       max neval
##   parse(text = text) 146.388 196.4995 254.677 333.5975 11487.912   100
##  internalParse(text)  75.639  96.5910 115.239 135.3885  5031.508   100
  

Как насчет

 test2 <- function(x, fixed, value, changed, seqvec){
    argList0 <- setNames(list(value),fixed)
    for(i in seq_along(seqvec)){
        argList <- c(list(x),argList0,setNames(list(seqvec[i]),changed))
        y <- do.call(dnorm,argList)
        sqrt(y*y   y*y)
    }
}
  

(Вероятно, вызывать числовой вектор — плохая идея seq : большую часть времени он будет работать, но при сбое произойдет сбой очень запутанным способом!)

Тесты:

 microbenchmark(test(-10:10,"mean",0,"sd",seq(1,5,length.out = 1001)),test2(-10:10,"mean",0,"sd",seq(1,5,length.out = 1001)),times = 10)
Unit: milliseconds
       expr       min        lq    median        uq       max neval
  test(...) 638.06613 643.23257 663.01249 682.08094 740.33703    10
 test2(...)  21.71302  22.07522  23.23993  28.96877  30.46197    10
  

Я думаю

 outer(-10:10,1:5,dnorm,mean=0)
  

должно быть намного быстрее, но не совсем ясно, что вы на самом деле пытаетесь сделать.