R vs. Julia
Lijun Wang
29 August, 2018
Let me use the Example 7.1.1 of Robert and Casella (2013) to compare the speed of R and Julia in MCMC.
We can implement this example with the following code. It is worth noting that I try to keep the same form between this two programming language as much as possible. For example, Julia does not support arbitrary mean and variance in its Gaussian sampling function randn(), while R can directly realize in rnorm(), but we both adopt the linear transformation of Gaussian distribution.

R

1
bigibbs <- function(T, rho)
2
{
3
x = numeric(T+1)
4
y = numeric(T+1)
5
for (t in 1:T){
6
x[t+1] = rnorm(1) * sqrt(1-rho^2) + rho*y[t]
7
y[t+1] = rnorm(1) * sqrt(1-rho^2) + rho*x[t+1]
8
}
9
return(list(x=x, y=y))
10
}
11
## example
12
res = bigibbs(2e6, 0.5)
Copied!

Julia

1
function bigibbs(T::Int64, rho::Float64)
2
x = ones(T+1)
3
y = ones(T+1)
4
for t = 1:T
5
x[t+1] = randn() * sqrt(1-rho^2) + rho*y[t]
6
y[t+1] = randn() * sqrt(1-rho^2) + rho*x[t+1]
7
end
8
return x, y
9
end
10
11
## example
12
x, y = bigibbs(Int64(2e6), 0.5)
Copied!

Results

The running environment is as follows:
    System: Ubuntu 18.04 (Windows subsystem for Linux)
    Processor: Intel(R) Core(TM) i7-6700 CPU @ 3.40GHz x 8
    Memory: 16 GiB
    R version: 3.4.4 (2018-03-15)
    Julia version: 1.0.0 (2018-08-08)
In terminal, use time command to get their running time:
Obviously, in our toy example, Julia outperforms much than R, nearly 16 times. Try another number of iterations, the results are similar.
Moreover, we can use PyPlot to plot in Julia v1.0, such as the histogram in this toy example:
Last modified 3yr ago
Copy link