Simulatin speed: GNU R vs Julia

November 22, 2013
By

(This article was first published on R snippets, and kindly contributed to R-bloggers)

Recently there is a lot of noise about Julia. I have decided to test its speed in simulation tasks on my toy Cont model. I thought I had vectorized my GNU R code pretty well, but Julia is much faster.

The model was described in my earlier posts so let us go down to a comparison:
Here is my GNU R code:

library(e1071)

cont.run <- function(burn.in, reps, n, d, l, s) {
  tr <- rep(0, n)
  r <- rep(0, reps)
  for (i in 1:reps) {
    sig <- rnorm(1, 0, d)
    mul <- 1
    if (sig < 0) {
        sig <- -sig
        mul <- -1
    }
    r[i] <- mul * sum(sig > tr) / (l * n)
    tr[runif(n) < s] <- abs(r[i])
  }
  return(kurtosis(r[burn.in:reps]))
}

system.time(replicate(10,
    cont.run(1000, 10000, 1000, 0.005, 10.0, 0.01)))

It's execution time is a bit below 10 seconds on my laptop.

An equivalent Julia code is the following:

using Distributions

function cont_run(burn_in, reps, n, d, l, s)
    tr = Array(Float64, n)
    r = Array(Float64, reps)
    for i in 1:reps
        aris = 0
        sig = randn() * d
        mul = 1
        if sig < 0
            sig = -sig
            mul = -1
        end
        for k in 1:n
            if sig > tr[k]
                aris += 1
            end
        end
        ari = aris / (l * n)
        r[i] = mul * ari
        for j in 1:n
            if rand() < s
                tr[j] = ari
            end
        end
    end
    kurtosis(r[burn_in:reps])
end

n = 10
t_start = time()
k = Array(Float64, n)
for i in 1:n
    k[i] = cont_run(1000, 10000, 1000, 0.005, 10.0, 0.01)
end
println(time() -t_start)

And on my machine it takes a bit less than 0.7 seconds to run.

So we get over tenfold speedup. This is a significant difference for simulation experiments.
I will have to dig more into Julia in the future.

To leave a comment for the author, please follow the link and comment on his blog: R snippets.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.