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: