[This article was first published on R snippets, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
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.Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
The model was described in my earlier posts so let us go down to a comparison:
library(e1071)< o:p>
cont.run <- function(burn.in, reps, n, d, l, s) {< o:p>
tr <- rep(0, n)< o:p>
r <- rep(0, reps)< o:p>
for (i in 1:reps) {< o:p>
sig <- rnorm(1, 0, d)< o:p>
mul <- 1< o:p>
if (sig < 0) {< o:p>
sig <- –sig< o:p>
mul <- –1< o:p>
}< o:p>
r[i] <- mul * sum(sig > tr) / (l * n)< o:p>
tr[runif(n) < s] <- abs(r[i])< o:p>
}< o:p>
return(kurtosis(r[burn.in:reps]))< o:p>
}< o:p>
system.time(replicate(10,< o:p>
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< o:p>
function cont_run(burn_in, reps, n, d, l, s)< o:p>
tr = Array(Float64, n)< o:p>
r = Array(Float64, reps)< o:p>
for i in 1:reps< o:p>
aris = 0< o:p>
sig = randn() * d< o:p>
mul = 1< o:p>
if sig < 0< o:p>
sig = –sig< o:p>
mul = –1< o:p>
end< o:p>
for k in 1:n< o:p>
if sig > tr[k]< o:p>
aris += 1< o:p>
end< o:p>
end< o:p>
ari = aris / (l * n)< o:p>
r[i] = mul * ari< o:p>
for j in 1:n< o:p>
if rand() < s< o:p>
tr[j] = ari< o:p>
end< o:p>
end< o:p>
end< o:p>
kurtosis(r[burn_in:reps])< o:p>
end< o:p>
n = 10< o:p>
t_start = time()< o:p>
k = Array(Float64, n)< o:p>
for i in 1:n< o:p>
k[i] = cont_run(1000, 10000, 1000, 0.005, 10.0, 0.01)< o:p>
end< o:p>
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 their blog: R snippets.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.