Advanced Usage

How to Define a Customized Distribution

Turing.jl supports the use of distributions from the Distributions.jl package. By extension it also supports the use of customized distributions, by defining them as subtypes of Distribution type of the Distributions.jl package, as well as corresponding functions.

Below shows a workflow of how to define a customized distribution, using a flat prior as a simple example.

1. Define the Distribution Type

First, define a type of the distribution, as a subtype of a corresponding distribution type in the Distributions.jl package.

struct Flat <: ContinuousUnivariateDistribution
  …
end

2. Implement Sampling and Evaluation of the log-pdf

Second, define rand() and logpdf(), which will be used to run the model.

Distributions.rand(d::Flat) = rand()
Distributions.logpdf{T<:Real}(d::Flat, x::T) = zero(x)

3. Define Helper Functions

In most cases, it may be required to define helper functions, such as the minimum, maximum, rand, and logpdf functions, among others.

3.1 Domain Transformation

Some helper functions are necessary for domain transformation. For univariate distributions, the necessary ones to implement are minimum() and maximum().

Distributions.minimum(d::Flat) = -Inf
Distributions.maximum(d::Flat) = +Inf

Functions for domain transformation which may be required by multivariate or matrix-variate distributions are size(d), link(d, x) and invlink(d, x). Please see Turing’s transform.jl for examples.

3.2 Vectorization Support

The vectorization syntax follows rv ~ [distribution], which requires rand() and logpdf() to be called on multiple data points at once. An appropriate implementation for Flat are shown below.

Distributions.rand(d::Flat, n::Int) = Vector([rand() for _ = 1:n])
Distributions.logpdf{T<:Real}(d::Flat, x::Vector{T}) = zero(x)

Model Internals

The @model macro accepts a function definition and generates a Turing.Model struct for use by the sampler. Models can be constructed by hand without the use of a macro. Taking the gdemo model as an example, the two code sections below (macro and macro-free) are equivalent.

using Turing

@model gdemo(x) = begin
  # Set priors.
  s ~ InverseGamma(2,3)
  m ~ Normal(0,sqrt(s))

  # Observe each value of x.
  [x ~ Normal(m, sqrt(s))]
end

sample(gdemo([1.5, 2.0]), HMC(1000, 0.1, 5))

using Turing

# Initialize a NamedTuple containing our data variables.
data = (x = [1.5, 2.0],)

# Create the model function.
mf(vi, sampler, model) = begin
    # Set the accumulated logp to zero.
    vi.logp = 0

    # If x is provided, use the provided values.
    # Otherwise, treat x as an empty vector with
    # two entries.
    if isdefined(model.data, :x)
        x = model.data.x
    else # x is a parameter
        x = model.defaults.x
    end

    # Assume s has an InverseGamma distribution.
    s, lp = Turing.assume(sampler,
                InverseGamma(2, 3),
                Turing.VarName(vi, [:c_s, :s], ""), vi)

    # Add the lp to the accumulated logp.
    vi.logp += lp

    # Assume m has a Normal distribution.
    m, lp = Turing.assume(sampler,
                Normal(0,sqrt(s)),
                Turing.VarName(vi, [:c_m, :m], ""), vi)

    # Add the lp to the accumulated logp.
    vi.logp += lp

    # Observe each value of x[i], according to a
    # Normal distribution.
    for i = 1:length(x)
        vi.logp += Turing.observe(sampler,
                   Normal(m, sqrt(s)),
                   x[i], vi)
    end
end

# Define the default value for x when missing
defaults = (x = Vector{Real}(undef, 2),)

# Instantiate a Model object.
model = Turing.Model{Tuple{:s, :m}, Tuple{:x}}(mf, data, defaults)

# Sample the model.
chain = sample(model, HMC(1000, 0.1, 5))

Note that the Turing.Model{Tuple{:s, :m}, Tuple{:x, :y}} accepts two parameter tuples. The first set, Tuple{:s, :m}, represents parameter variables that will be generated by the model, while the second (Tuple{:x}) contains the variables to be observed.

Task Copying

Turing copies Julia tasks to deliver efficient inference algorithms, but it also provides alternative slower implementation as a fallback. Task copying is enabled by default. Task copying requires we use the CTask facility which is provided by Libtask to create tasks.

Maximum a Posteriori Estimation

Turing does not currently have built-in methods for calculating the maximum a posteriori (MAP) for a model. This is a goal for Turing’s implementation (see this issue), but for the moment, we present here a method for estimating the MAP using Optim.jl.

using Turing

# Define the simple gdemo model.
@model gdemo(x, y) = begin
 s ~ InverseGamma(2,3)
 m ~ Normal(0,sqrt(s))
 x ~ Normal(m, sqrt(s))
 y ~ Normal(m, sqrt(s))
 return s, m
end

# Define our data points.
x = 1.5
y = 2.0

# Set up the model call, sample from the prior.
model = gdemo(x, y)
vi = Turing.VarInfo()
model(vi, Turing.SampleFromPrior())

# Define a function to optimize.
function nlogp(sm)
 vi.vals .= sm
 model(vi, Turing.SampleFromPrior())
 -vi.logp
end

# Import Optim.jl.
using Optim

# Create a starting point, call the optimizer.
sm_0 = Float64.(vi.vals)
lb = [0.0, -Inf]
ub = [Inf, Inf]
result = optimize(nlogp, lb, ub, sm_0, Fminbox())