A Functional Implementation of Slice Sampling in Julia

I’ve been using Clojure and Scala most recently for algorithmic work but I needed to return to Julia momentarily in order to finish a project.

Julia is not the most natural language in which to perform functional programming; it does not support tail call elimination, no native support for immutable data structures, currying interferes with its multiple dispatch model and it seems most of the design choices favoured imperativeness. Most functional functionality is provided by the third-party Lazy, FunctionalCollections and Match packages which provide lazy lists, functional data structures and pattern matching respectively.

“Stepping Out” Slice Sampling

No variables were mutated in the making of this code

function SliceSampling(g,w,x)
    function lowerbound(L)
        g(L)<y ? L : lowerbound(L-w)
    end
    function upperbound(R)
        g(R)<y ? R : upperbound(R+w)
    end
    function shrinksample(L,R)
        z=rand(Uniform(L,R))
        if g(z)>y
            z
        elseif z>x
            shrinksample(L,z)
        else
            shrinksample(z,R)
        end
    end
    y=-1*rand(Exponential(1))+g(x)
    U=rand(Uniform(0,1))
    shrinksample(lowerbound(x-U*w),upperbound(x+(1-U)*w))
end

g is the logdensity from which one is trying to sample, w is an initial width and x is the current parameter value.

Preamble

Loading required libraries and defining two function. Partial is my own implementation of partial application, a function present in Clojure but lacking in Julia. The toArray function just converts a lazy list into an array so that it can be plotted with Gadfly (this is done imperatively because the corresponding functional implementation results in a stack overflow due to Julia’s lack of support for tail call elimination).

using Distributions
using Lazy
using DataFrames
using Gadfly

function partial(f,a...)
    (b...) -> f(a...,b...)
end

function toArray(l)
  a=Array{Float64}(length(l))
  for i=1:length(l)
    a[i]=first(l)
    l=Lazy.tail(l)
  end
  a
end

Lazy Evaluation

Lets consider a test pdf from which we can sample directly to form a unit test for the above code. Here is a mixture density of three Normals…

mixture_distribution=MixtureModel([Normal(-1,0.2),Normal(0,0.2),Normal(1,0.2)],[0.2,0.6,0.2])
f = x->pdf(mixture_distribution,x)
g = x->logpdf(mixture_distribution,x)
plot(f, -4, 4)

Given a current parameter value, 0.0, one can generate a new value given the density f and an initial width 1.0 by

SliceSampling(g,1.0,0.0)

0.2485248928906587

It would be nice to create a new function, specific to f, to generate random variates which masks the distracting implementation details such as the initial width. This is the reason for the partial application function I created earlier. partial(SliceSampling,f,1.0) returns a new function which maps the current parameter value to a new value.

sample_f = partial(SliceSampling,g,1.0)
sample_f(0.5)

0.6647567078804801

This function can be used to define an infinite list of values that is lazily evaluated. The iterate function for a given function g generates an infinite lazy list as (x, g(x), g(g(x)),…). One can then take values from the beginning of the list using the take function i.e.

stream_f = iterate(sample_f,0.0);
take(stream_f,10)

(0.0 -0.21844923549096884 1.0302399899983261 1.0129526649079534 -0.5690927045181375 -1.1228947590022138 -1.1224701288130174 -0.11733086188135866 0.012546600939988967 1.1702840100996514)

take(stream_f,number_to_take) is not the most visually pleasing implementation of generating random variates from f, so once again we can use partial application to get a new function which is a lot closer to functions like rnorm in R.

rand_f = partial(take,stream_f)
rand_f(10)

(0.0 -0.21844923549096884 -0.6311260691954774 1.1600742443380705 -0.1729859710235011 -0.1801167762470986 0.15963462381533539 -0.31375918237598716 0.22438826503951037 0.02349035294676849)

Visualization

Finally, here are the density estimates from 10000 draws using slice sampling and direct methods.

plot(layer(x=(@> 10000 rand_f toArray), Geom.density,Theme(default_color=color("yellow")) ),
layer(x=rand(mixture_distribution,10000), Geom.density,Theme(default_color=color("magenta")) ),
Guide.manual_color_key("Density", [ "Slice Sampling", "Direct"], ["yellow",  "magenta"])
)