Más sobre Julia (II): mi primer programa

A las entradas que he hecho sobre Julia estos últimos días, quiero añadir esta en la que publico mi primer programa en dicho lenguaje.

Me ha dado por reimplementar el programa para realizar un muestreo de Gibbs que aparece en Gibbs sampler in various languages.

Lo primero ha sido instalar Julia, para lo que basta con seguir las instrucciones que aparecen en su página de github. Y aviso: tarda bastante en descargar y compilar todas sus dependencias.

El código de mi programa, para que quede constancia pública de que vengo a ser un pionero en España usándolo, es:

function gibbs(n, thin)
        x = 0.0
        y = 0.0
        for i=1:n
                for j=1:thin
                        x = randg(3) / (y * y + 4)
                        y = randn() / sqrt(2*x + 2) + 1 / (x + 1)
                end
                println(i, " ", x, " ", y)
        end
end

gibbs(50000, 1000)

El programa en cuestión, tal cual viene en la página que cito arriba, en R, ha tardado en ejecutarse

$ time Rscript mcmc00.R > data.tab

real 8m31.259s
user 8m29.668s
sys 0m0.948s

Y en Julia,

$ time julia mcmc00.j > data.tab.j

real 0m9.268s
user 0m9.113s s
ys 0m0.144s

Es decir, del orden de 60 veces menos. Es de reseñar que este tipo de algoritmos son de los menos indicados para ser ejecutados con R: son imposibles de vectorizar. No en vano, MCMC (tal como aparece en los nombres de los programas) significa Markov Chain Monte Carlo y en una cadena de Markov, cada valor depende del estado anterior. Hay que morir a hierro y construir un bucle.

El resultado de ambas simulaciones puede apreciarse en el siguiente gráfico:

Es cierto que en este ejemplo concreto, los números aleatorios podrían precalcularse en una única llamada a rnorm o rgamma, pero no he observado una ganancia sustancial de tiempos, que se pierden, fundamentalmene, en el código interpretado dentro del bucle.