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