1

I'm doing some simulations with different probability distributions and with different parameters. Example: distribiution - normal, parameter mean and sd. Mean is fixed (eg 0) and sd changes from eg 1 to 5. I wrote some function to deal with this, simple version which will give you the idea is below:

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)
    }
}

In my function parse is taking half of the simulation time (13sec out of 25sec). I need to improve this, any idea how? You can run Rprof function and check the result:

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

Rprof results of my simulation (only part)

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

One idea to speed up this is use internal parse function:

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
1
  • PS did you really mean y*y + y*y and not x*x + y*y ?
    – Ben Bolker
    Jan 12, 2015 at 14:28

1 Answer 1

2

How about

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)
    }
}

(It's probably a bad idea to call a numeric vector seq: it will work most of the time, but when it fails it will fail in a very confusing way!)

Benchmarks:

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

I think

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

should be much faster, but it's not entirely clear what you're actually trying to do.

2
  • Added some benchmarks. Your solutions looks really promising :) Jul 3, 2014 at 22:52
  • 1
    Simulation time has decreased from 25sec to 9sec! Thank you very much! Jul 3, 2014 at 23:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.