Guide
Basics
Introduction
A probabilistic program is Julia code wrapped in a @model
macro. It can use arbitrary Julia code, but to ensure correctness of inference it should not have external effects or modify global state. Stack-allocated variables are safe, but mutable heap-allocated objects may lead to subtle bugs when using task copying. By default Libtask deepcopies Array
and Dict
objects when copying task to avoid bugs with data stored in mutable structure in Turing models.
To specify distributions of random variables, Turing programs should use the ~
notation:
x ~ distr
where x
is a symbol and distr
is a distribution. If x
is undefined in the model function, inside the probabilistic program, this puts a random variable named x
, distributed according to distr
, in the current scope. distr
can be a value of any type that implements rand(distr)
, which samples a value from the distribution distr
. If x
is defined, this is used for conditioning in a style similar to Anglican (another PPL). In this case, x
is an observed value, assumed to have been drawn from the distribution distr
. The likelihood is computed using logpdf(distr,y)
. The observe statements should be arranged so that every possible run traverses all of them in exactly the same order. This is equivalent to demanding that they are not placed inside stochastic control flow.
Available inference methods include Importance Sampling (IS), Sequential Monte Carlo (SMC), Particle Gibbs (PG), Hamiltonian Monte Carlo (HMC), Hamiltonian Monte Carlo with Dual Averaging (HMCDA) and The No-U-Turn Sampler (NUTS).
Simple Gaussian Demo
Below is a simple Gaussian demo illustrate the basic usage of Turing.jl.
# Import packages.
using Turing
using StatsPlots
# Define a simple Normal model with unknown mean and variance.
@model function gdemo(x, y)
s² ~ InverseGamma(2, 3)
m ~ Normal(0, sqrt(s²))
x ~ Normal(m, sqrt(s²))
return y ~ Normal(m, sqrt(s²))
end
gdemo (generic function with 2 methods)
Note: As a sanity check, the prior expectation of s²
is mean(InverseGamma(2, 3)) = 3/(2 - 1) = 3
and the prior expectation of m
is 0. This can be easily checked using Prior
:
p1 = sample(gdemo(missing, missing), Prior(), 100000)
Chains MCMC chain (100000×5×1 Array{Float64, 3}):
Iterations = 1:1:100000
Number of chains = 1
Samples per chain = 100000
Wall duration = 1.72 seconds
Compute duration = 1.72 seconds
parameters = s², m, x, y
internals = lp
Summary Statistics
parameters mean std mcse ess_bulk ess_tail
rha ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Flo
at6 ⋯
s² 3.0088 12.2842 0.0389 98431.3947 99067.4063 1.
000 ⋯
m 0.0120 1.7364 0.0055 101749.7416 97648.9695 1.
000 ⋯
x 0.0134 2.4815 0.0078 101451.9995 99792.9764 1.
000 ⋯
y 0.0018 2.4357 0.0077 100267.8453 97045.8311 1.
000 ⋯
2 columns om
itted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
s² 0.5423 1.1057 1.7757 3.0986 12.2915
m -3.3670 -0.8953 0.0075 0.9152 3.3949
x -4.7249 -1.2641 0.0059 1.2807 4.7926
y -4.8172 -1.2773 0.0080 1.2791 4.7550
We can perform inference by using the sample
function, the first argument of which is our probabilistic program and the second of which is a sampler. More information on each sampler is located in the API.
# Run sampler, collect results.
c1 = sample(gdemo(1.5, 2), SMC(), 1000)
c2 = sample(gdemo(1.5, 2), PG(10), 1000)
c3 = sample(gdemo(1.5, 2), HMC(0.1, 5), 1000)
c4 = sample(gdemo(1.5, 2), Gibbs(PG(10, :m), HMC(0.1, 5, :s²)), 1000)
c5 = sample(gdemo(1.5, 2), HMCDA(0.15, 0.65), 1000)
c6 = sample(gdemo(1.5, 2), NUTS(0.65), 1000)
Chains MCMC chain (1000×14×1 Array{Float64, 3}):
Iterations = 501:1:1500
Number of chains = 1
Samples per chain = 1000
Wall duration = 1.62 seconds
Compute duration = 1.62 seconds
parameters = s², m
internals = lp, n_steps, is_accept, acceptance_rate, log_density, h
amiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error,
tree_depth, numerical_error, step_size, nom_step_size
Summary Statistics
parameters mean std mcse ess_bulk ess_tail rhat
e ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64
⋯
s² 1.9868 1.8407 0.0793 507.1843 530.8663 1.0005
⋯
m 1.1809 0.7571 0.0307 618.5481 589.6552 1.0002
⋯
1 column om
itted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
s² 0.5822 1.0126 1.5065 2.3121 5.9207
m -0.3145 0.7107 1.1555 1.6179 2.7291
The MCMCChains
module (which is re-exported by Turing) provides plotting tools for the Chain
objects returned by a sample
function. See the MCMCChains repository for more information on the suite of tools available for diagnosing MCMC chains.
# Summarise results
describe(c3)
# Plot results
plot(c3)
savefig("gdemo-plot.png")
The arguments for each sampler are:
- SMC: number of particles.
- PG: number of particles, number of iterations.
- HMC: leapfrog step size, leapfrog step numbers.
- Gibbs: component sampler 1, component sampler 2, ...
- HMCDA: total leapfrog length, target accept ratio.
- NUTS: number of adaptation steps (optional), target accept ratio.
For detailed information on the samplers, please review Turing.jl's API documentation.
Modelling Syntax Explained
Using this syntax, a probabilistic model is defined in Turing. The model function generated by Turing can then be used to condition the model onto data. Subsequently, the sample function can be used to generate samples from the posterior distribution.
In the following example, the defined model is conditioned to the data (arg1 = 1, arg2 = 2) by passing (1, 2) to the model function.
@model function model_name(arg_1, arg_2)
return ...
end
The conditioned model can then be passed onto the sample function to run posterior inference.
model_func = model_name(1, 2)
chn = sample(model_func, HMC(..)) # Perform inference by sampling using HMC.
The returned chain contains samples of the variables in the model.
var_1 = mean(chn[:var_1]) # Taking the mean of a variable named var_1.
The key (:var_1
) can be a Symbol
or a String
. For example, to fetch x[1]
, one can use chn[Symbol("x[1]")]
or chn["x[1]"]
.
If you want to retrieve all parameters associated with a specific symbol, you can use group
. As an example, if you have the
parameters "x[1]"
, "x[2]"
, and "x[3]"
, calling group(chn, :x)
or group(chn, "x")
will return a new chain with only "x[1]"
, "x[2]"
, and "x[3]"
.
Turing does not have a declarative form. More generally, the order in which you place the lines of a @model
macro matters. For example, the following example works:
# Define a simple Normal model with unknown mean and variance.
@model function model_function(y)
s ~ Poisson(1)
y ~ Normal(s, 1)
return y
end
sample(model_function(10), SMC(), 100)
Chains MCMC chain (100×3×1 Array{Float64, 3}):
Log evidence = -18.015931779125733
Iterations = 1:1:100
Number of chains = 1
Samples per chain = 100
Wall duration = 1.58 seconds
Compute duration = 1.58 seconds
parameters = s
internals = lp, weight
Summary Statistics
parameters mean std mcse ess_bulk ess_tail rhat
e ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64
⋯
s 4.9900 0.1000 0.0100 100.0801 NaN 1.0000
⋯
1 column om
itted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
s 5.0000 5.0000 5.0000 5.0000 5.0000
But if we switch the s ~ Poisson(1)
and y ~ Normal(s, 1)
lines, the model will no longer sample correctly:
# Define a simple Normal model with unknown mean and variance.
@model function model_function(y)
y ~ Normal(s, 1)
s ~ Poisson(1)
return y
end
sample(model_function(10), SMC(), 100)
Sampling Multiple Chains
Turing supports distributed and threaded parallel sampling. To do so, call sample(model, sampler, parallel_type, n, n_chains)
, where parallel_type
can be either MCMCThreads()
or MCMCDistributed()
for thread and parallel sampling, respectively.
Having multiple chains in the same object is valuable for evaluating convergence. Some diagnostic functions like gelmandiag
require multiple chains.
If you do not want parallelism or are on an older version Julia, you can sample multiple chains with the mapreduce
function:
# Replace num_chains below with however many chains you wish to sample.
chains = mapreduce(c -> sample(model_fun, sampler, 1000), chainscat, 1:num_chains)
The chains
variable now contains a Chains
object which can be indexed by chain. To pull out the first chain from the chains
object, use chains[:,:,1]
. The method is the same if you use either of the below parallel sampling methods.
Multithreaded sampling
If you wish to perform multithreaded sampling and are running Julia 1.3 or greater, you can call sample
with the following signature:
using Turing
@model function gdemo(x)
s² ~ InverseGamma(2, 3)
m ~ Normal(0, sqrt(s²))
for i in eachindex(x)
x[i] ~ Normal(m, sqrt(s²))
end
end
model = gdemo([1.5, 2.0])
# Sample four chains using multiple threads, each with 1000 samples.
sample(model, NUTS(), MCMCThreads(), 1000, 4)
Be aware that Turing cannot add threads for you -- you must have started your Julia instance with multiple threads to experience any kind of parallelism. See the Julia documentation for details on how to achieve this.
Distributed sampling
To perform distributed sampling (using multiple processes), you must first import Distributed
.
Process parallel sampling can be done like so:
# Load Distributed to add processes and the @everywhere macro.
using Distributed
# Load Turing.
using Turing
# Add four processes to use for sampling.
addprocs(4; exeflags="--project=$(Base.active_project())")
# Initialize everything on all the processes.
# Note: Make sure to do this after you've already loaded Turing,
# so each process does not have to precompile.
# Parallel sampling may fail silently if you do not do this.
@everywhere using Turing
# Define a model on all processes.
@everywhere @model function gdemo(x)
s² ~ InverseGamma(2, 3)
m ~ Normal(0, sqrt(s²))
for i in eachindex(x)
x[i] ~ Normal(m, sqrt(s²))
end
end
# Declare the model instance everywhere.
@everywhere model = gdemo([1.5, 2.0])
# Sample four chains using multiple processes, each with 1000 samples.
sample(model, NUTS(), MCMCDistributed(), 1000, 4)
Sampling from an Unconditional Distribution (The Prior)
Turing allows you to sample from a declared model's prior. If you wish to draw a chain from the prior to inspect your prior distributions, you can simply run
chain = sample(model, Prior(), n_samples)
You can also run your model (as if it were a function) from the prior distribution, by calling the model without specifying inputs or a sampler. In the below example, we specify a gdemo
model which returns two variables, x
and y
. The model includes x
and y
as arguments, but calling the function without passing in x
or y
means that Turing's compiler will assume they are missing values to draw from the relevant distribution. The return
statement is necessary to retrieve the sampled x
and y
values.
Assign the function with missing
inputs to a variable, and Turing will produce a sample from the prior distribution.
@model function gdemo(x, y)
s² ~ InverseGamma(2, 3)
m ~ Normal(0, sqrt(s²))
x ~ Normal(m, sqrt(s²))
y ~ Normal(m, sqrt(s²))
return x, y
end
gdemo (generic function with 2 methods)
Assign the function with missing
inputs to a variable, and Turing will produce a sample from the prior distribution.
# Samples from p(x,y)
g_prior_sample = gdemo(missing, missing)
g_prior_sample()
(-4.399191238263897, -0.06599138237512658)
Sampling from a Conditional Distribution (The Posterior)
Treating observations as random variables
Inputs to the model that have a value missing
are treated as parameters, aka random variables, to be estimated/sampled. This can be useful if you want to simulate draws for that parameter, or if you are sampling from a conditional distribution. Turing supports the following syntax:
@model function gdemo(x, ::Type{T}=Float64) where {T}
if x === missing
# Initialize `x` if missing
x = Vector{T}(undef, 2)
end
s² ~ InverseGamma(2, 3)
m ~ Normal(0, sqrt(s²))
for i in eachindex(x)
x[i] ~ Normal(m, sqrt(s²))
end
end
# Construct a model with x = missing
model = gdemo(missing)
c = sample(model, HMC(0.01, 5), 500)
Chains MCMC chain (500×14×1 Array{Float64, 3}):
Iterations = 1:1:500
Number of chains = 1
Samples per chain = 500
Wall duration = 1.48 seconds
Compute duration = 1.48 seconds
parameters = s², m, x[1], x[2]
internals = lp, n_steps, is_accept, acceptance_rate, log_density, h
amiltonian_energy, hamiltonian_energy_error, numerical_error, step_size, no
m_step_size
Summary Statistics
parameters mean std mcse ess_bulk ess_tail rhat
e ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64
⋯
s² 2.9678 1.5805 1.4518 1.3576 21.2342 2.1102
⋯
m -2.1296 0.5047 0.3984 1.6426 14.2176 1.6582
⋯
x[1] -1.0668 0.4309 0.3121 2.2861 20.4932 1.3787
⋯
x[2] -1.9494 0.2403 0.1924 1.5666 20.5654 1.8044
⋯
1 column om
itted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
s² 1.0698 1.4199 3.1969 3.9949 6.2587
m -2.7673 -2.4849 -2.3049 -1.8122 -1.0943
x[1] -1.7950 -1.4818 -1.0092 -0.6820 -0.3575
x[2] -2.4318 -2.1281 -1.9049 -1.7819 -1.5198
Note the need to initialize x
when missing since we are iterating over its elements later in the model. The generated values for x
can be extracted from the Chains
object using c[:x]
.
Turing also supports mixed missing
and non-missing
values in x
, where the missing ones will be treated as random variables to be sampled while the others get treated as observations. For example:
@model function gdemo(x)
s² ~ InverseGamma(2, 3)
m ~ Normal(0, sqrt(s²))
for i in eachindex(x)
x[i] ~ Normal(m, sqrt(s²))
end
end
# x[1] is a parameter, but x[2] is an observation
model = gdemo([missing, 2.4])
c = sample(model, HMC(0.01, 5), 500)
Chains MCMC chain (500×13×1 Array{Float64, 3}):
Iterations = 1:1:500
Number of chains = 1
Samples per chain = 500
Wall duration = 1.43 seconds
Compute duration = 1.43 seconds
parameters = s², m, x[1]
internals = lp, n_steps, is_accept, acceptance_rate, log_density, h
amiltonian_energy, hamiltonian_energy_error, numerical_error, step_size, no
m_step_size
Summary Statistics
parameters mean std mcse ess_bulk ess_tail rhat
e ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64
⋯
s² 1.3781 0.4416 0.2028 6.4822 30.9557 1.2614
⋯
m 0.1970 0.2748 0.0794 11.5829 37.9162 1.0441
⋯
x[1] -0.0341 0.3759 0.2862 1.8543 20.6421 1.5484
⋯
1 column om
itted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
s² 0.7452 1.0905 1.2838 1.5927 2.3255
m -0.4449 0.1057 0.2209 0.3650 0.6837
x[1] -0.6546 -0.3279 -0.0249 0.1583 0.6807
Default Values
Arguments to Turing models can have default values much like how default values work in normal Julia functions. For instance, the following will assign missing
to x
and treat it as a random variable. If the default value is not missing
, x
will be assigned that value and will be treated as an observation instead.
using Turing
@model function generative(x=missing, ::Type{T}=Float64) where {T<:Real}
if x === missing
# Initialize x when missing
x = Vector{T}(undef, 10)
end
s² ~ InverseGamma(2, 3)
m ~ Normal(0, sqrt(s²))
for i in 1:length(x)
x[i] ~ Normal(m, sqrt(s²))
end
return s², m
end
m = generative()
chain = sample(m, HMC(0.01, 5), 1000)
Chains MCMC chain (1000×22×1 Array{Float64, 3}):
Iterations = 1:1:1000
Number of chains = 1
Samples per chain = 1000
Wall duration = 2.8 seconds
Compute duration = 2.8 seconds
parameters = s², m, x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8],
x[9], x[10]
internals = lp, n_steps, is_accept, acceptance_rate, log_density, h
amiltonian_energy, hamiltonian_energy_error, numerical_error, step_size, no
m_step_size
Summary Statistics
parameters mean std mcse ess_bulk ess_tail rhat
e ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64
⋯
s² 1.8298 0.6807 0.2239 8.1502 32.1267 1.0221
⋯
m 0.5251 0.4781 0.3069 2.6269 19.1169 1.9548
⋯
x[1] 1.6652 0.5137 0.1935 7.4635 14.6955 1.0493
⋯
x[2] 2.2582 0.4409 0.2889 2.5551 41.6556 2.1114
⋯
x[3] 1.3261 0.3186 0.1450 4.8317 23.8646 1.2789
⋯
x[4] 0.4437 0.4003 0.2279 3.5036 15.2000 1.4068
⋯
x[5] 0.9751 0.5595 0.3180 3.5974 20.4467 1.4211
⋯
x[6] -0.7887 0.9681 0.6234 2.5598 20.3470 1.9945
⋯
x[7] -0.7925 0.5634 0.2325 7.7980 20.3102 1.1742
⋯
x[8] 0.0311 0.3516 0.1177 9.3554 14.6621 1.0151
⋯
x[9] -0.4124 0.3103 0.0909 11.7480 21.5038 1.0069
⋯
x[10] -0.5153 0.5122 0.2918 3.1825 20.4875 1.4673
⋯
1 column om
itted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
s² 0.8404 1.2987 1.6768 2.2308 3.3380
m -0.2981 0.1201 0.4692 0.9802 1.3657
x[1] 0.6698 1.2801 1.7757 2.1017 2.4874
x[2] 1.4809 1.8996 2.2318 2.6190 3.0276
x[3] 0.8294 1.0448 1.3172 1.6028 1.9281
x[4] -0.3291 0.1074 0.5452 0.7726 1.0685
x[5] 0.0676 0.4128 1.0454 1.3136 2.0963
x[6] -2.2310 -1.6432 -0.9077 0.1924 1.0623
x[7] -2.2262 -1.0075 -0.6444 -0.4235 -0.0190
x[8] -0.7566 -0.2200 0.0449 0.3194 0.6116
x[9] -0.9650 -0.5982 -0.4246 -0.2361 0.2853
x[10] -1.6072 -0.8621 -0.4466 -0.0912 0.2863
Access Values inside Chain
You can access the values inside a chain several ways:
- Turn them into a
DataFrame
object - Use their raw
AxisArray
form - Create a three-dimensional
Array
object
For example, let c
be a Chain
:
DataFrame(c)
convertsc
to aDataFrame
,c.value
retrieves the values insidec
as anAxisArray
, andc.value.data
retrieves the values insidec
as a 3DArray
.
Variable Types and Type Parameters
The element type of a vector (or matrix) of random variables should match the eltype
of the its prior distribution, <: Integer
for discrete distributions and <: AbstractFloat
for continuous distributions. Moreover, if the continuous random variable is to be sampled using a Hamiltonian sampler, the vector's element type needs to either be:
Real
to enable auto-differentiation through the model which uses special number types that are sub-types ofReal
, or- Some type parameter
T
defined in the model header using the type parameter syntax, e.g.function gdemo(x, ::Type{T} = Float64) where {T}
. Similarly, when using a particle sampler, the Julia variable used should either be: - An
Array
, or - An instance of some type parameter
T
defined in the model header using the type parameter syntax, e.g.function gdemo(x, ::Type{T} = Vector{Float64}) where {T}
.
Querying Probabilities from Model or Chain
Consider first the following simplified gdemo
model:
@model function gdemo0(x)
s ~ InverseGamma(2, 3)
m ~ Normal(0, sqrt(s))
return x ~ Normal(m, sqrt(s))
end
# Instantiate three models, with different value of x
model1 = gdemo0(1)
model4 = gdemo0(4)
model10 = gdemo0(10)
DynamicPPL.Model{typeof(Main.var"##WeaveSandBox#501".gdemo0), (:x,), (), ()
, Tuple{Int64}, Tuple{}, DynamicPPL.DefaultContext}(Main.var"##WeaveSandBox
#501".gdemo0, (x = 10,), NamedTuple(), DynamicPPL.DefaultContext())
Now, query the instantiated models: compute the likelihood of x = 1.0
given the values of s = 1.0
and m = 1.0
for the parameters:
prob"x = 1.0 | model = model1, s = 1.0, m = 1.0"
0.39894228040143265
prob"x = 1.0 | model = model4, s = 1.0, m = 1.0"
0.39894228040143265
prob"x = 1.0 | model = model10, s = 1.0, m = 1.0"
0.39894228040143265
Notice that even if we use three models, instantiated with three different values of x
, we should obtain the same likelihood. We can easily verify that value in this case:
pdf(Normal(1.0, 1.0), 1.0)
0.3989422804014327
Let us now consider the following gdemo
model:
@model function gdemo(x, y)
s² ~ InverseGamma(2, 3)
m ~ Normal(0, sqrt(s²))
x ~ Normal(m, sqrt(s²))
return y ~ Normal(m, sqrt(s²))
end
# Instantiate the model.
model = gdemo(2.0, 4.0)
DynamicPPL.Model{typeof(Main.var"##WeaveSandBox#501".gdemo), (:x, :y), (),
(), Tuple{Float64, Float64}, Tuple{}, DynamicPPL.DefaultContext}(Main.var"#
#WeaveSandBox#501".gdemo, (x = 2.0, y = 4.0), NamedTuple(), DynamicPPL.Defa
ultContext())
The following are examples of valid queries of the Turing
model or chain:
-
prob"x = 1.0, y = 1.0 | model = model, s = 1.0, m = 1.0"
calculates the likelihood ofx = 1
andy = 1
givens = 1
andm = 1
. -
prob"s² = 1.0, m = 1.0 | model = model, x = nothing, y = nothing"
calculates the joint probability ofs = 1
andm = 1
ignoringx
andy
.x
andy
are ignored so they can be optionally dropped from the RHS of|
, but it is recommended to define them. -
prob"s² = 1.0, m = 1.0, x = 1.0 | model = model, y = nothing"
calculates the joint probability ofs = 1
,m = 1
andx = 1
ignoringy
. -
prob"s² = 1.0, m = 1.0, x = 1.0, y = 1.0 | model = model"
calculates the joint probability of all the variables. -
After the MCMC sampling, given a
chain
,prob"x = 1.0, y = 1.0 | chain = chain, model = model"
calculates the element-wise likelihood ofx = 1.0
andy = 1.0
for each sample inchain
. -
If
save_state=true
was used during sampling (i.e.,sample(model, sampler, N; save_state=true)
), you can simply doprob"x = 1.0, y = 1.0 | chain = chain"
.
In all the above cases, logprob
can be used instead of prob
to calculate the log probabilities instead.
Maximum likelihood and maximum a posterior estimates
Turing provides support for two mode estimation techniques, maximum likelihood estimation (MLE) and maximum a posterior (MAP) estimation. Optimization is performed by the Optim.jl package. Mode estimation is currently a optional tool, and will not be available to you unless you have manually installed Optim and loaded the package with a using
statement. To install Optim, run import Pkg; Pkg.add("Optim")
.
Mode estimation only works when all model parameters are continuous -- discrete parameters cannot be estimated with MLE/MAP as of yet.
To understand how mode estimation works, let us first load Turing and Optim to enable mode estimation, and then declare a model:
# Note that loading Optim explicitly is required for mode estimation to function,
# as Turing does not load the opimization suite unless Optim is loaded as well.
using Turing
using Optim
@model function gdemo(x)
s² ~ InverseGamma(2, 3)
m ~ Normal(0, sqrt(s²))
for i in eachindex(x)
x[i] ~ Normal(m, sqrt(s²))
end
end
gdemo (generic function with 6 methods)
Once the model is defined, we can construct a model instance as we normally would:
# Create some data to pass to the model.
data = [1.5, 2.0]
# Instantiate the gdemo model with our data.
model = gdemo(data)
DynamicPPL.Model{typeof(Main.var"##WeaveSandBox#501".gdemo), (:x,), (), (),
Tuple{Vector{Float64}}, Tuple{}, DynamicPPL.DefaultContext}(Main.var"##Wea
veSandBox#501".gdemo, (x = [1.5, 2.0],), NamedTuple(), DynamicPPL.DefaultCo
ntext())
Mode estimation is typically quick and easy at this point. Turing extends the function Optim.optimize
and accepts the structs MLE()
or MAP()
, which inform Turing whether to provide an MLE or MAP estimate, respectively. By default, the LBFGS optimizer is used, though this can be changed. Basic usage is:
# Generate a MLE estimate.
mle_estimate = optimize(model, MLE())
# Generate a MAP estimate.
map_estimate = optimize(model, MAP())
ModeResult with maximized lp of -4.62
[0.9074074068278611, 1.1666666673479593]
If you wish to change to a different optimizer, such as NelderMead
, simply place your optimizer in the third argument slot:
# Use NelderMead
mle_estimate = optimize(model, MLE(), NelderMead())
# Use SimulatedAnnealing
mle_estimate = optimize(model, MLE(), SimulatedAnnealing())
# Use ParticleSwarm
mle_estimate = optimize(model, MLE(), ParticleSwarm())
# Use Newton
mle_estimate = optimize(model, MLE(), Newton())
# Use AcceleratedGradientDescent
mle_estimate = optimize(model, MLE(), AcceleratedGradientDescent())
Some methods may have trouble calculating the mode because not enough iterations were allowed, or the target function moved upwards between function calls. Turing will warn you if Optim fails to converge by running Optim.converge
. A typical solution to this might be to add more iterations, or allow the optimizer to increase between function iterations:
# Increase the iterations and allow function eval to increase between calls.
mle_estimate = optimize(
model, MLE(), Newton(), Optim.Options(; iterations=10_000, allow_f_increases=true)
)
More options for Optim are available here.
Analyzing your mode estimate
Turing extends several methods from StatsBase
that can be used to analyze your mode estimation results. Methods implemented include vcov
, informationmatrix
, coeftable
, params
, and coef
, among others.
For example, let's examine our ML estimate from above using coeftable
:
# Import StatsBase to use it's statistical methods.
using StatsBase
# Print out the coefficient table.
coeftable(mle_estimate)
─────────────────────────────
estimate stderror tstat
─────────────────────────────
s 0.0625 0.0625 1.0
m 1.75 0.176777 9.8995
─────────────────────────────
Standard errors are calculated from the Fisher information matrix (inverse Hessian of the log likelihood or log joint). t-statistics will be familiar to frequentist statisticians. Warning -- standard errors calculated in this way may not always be appropriate for MAP estimates, so please be cautious in interpreting them.
Sampling with the MAP/MLE as initial states
You can begin sampling your chain from an MLE/MAP estimate by extracting the vector of parameter values and providing it to the sample
function with the keyword init_params
. For example, here is how to sample from the full posterior using the MAP estimate as the starting point:
# Generate an MAP estimate.
map_estimate = optimize(model, MAP())
# Sample with the MAP estimate as the starting point.
chain = sample(model, NUTS(), 1_000; init_params=map_estimate.values.array)
Beyond the Basics
Compositional Sampling Using Gibbs
Turing.jl provides a Gibbs interface to combine different samplers. For example, one can combine an HMC
sampler with a PG
sampler to run inference for different parameters in a single model as below.
@model function simple_choice(xs)
p ~ Beta(2, 2)
z ~ Bernoulli(p)
for i in 1:length(xs)
if z == 1
xs[i] ~ Normal(0, 1)
else
xs[i] ~ Normal(2, 1)
end
end
end
simple_choice_f = simple_choice([1.5, 2.0, 0.3])
chn = sample(simple_choice_f, Gibbs(HMC(0.2, 3, :p), PG(20, :z)), 1000)
Chains MCMC chain (1000×3×1 Array{Float64, 3}):
Iterations = 1:1:1000
Number of chains = 1
Samples per chain = 1000
Wall duration = 16.64 seconds
Compute duration = 16.64 seconds
parameters = p, z
internals = lp
Summary Statistics
parameters mean std mcse ess_bulk ess_tail rhat
e ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64
⋯
p 0.4096 0.2074 0.0318 54.7426 143.2946 1.0747
⋯
z 0.1520 0.3592 0.0193 345.1864 NaN 1.0091
⋯
1 column om
itted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
p 0.0633 0.2501 0.3835 0.5554 0.8294
z 0.0000 0.0000 0.0000 0.0000 1.0000
The Gibbs
sampler can be used to specify unique automatic differentiation backends for different variable spaces. Please see the Automatic Differentiation article for more.
For more details of compositional sampling in Turing.jl, please check the corresponding paper.
Working with filldist and arraydist
Turing provides filldist(dist::Distribution, n::Int)
and arraydist(dists::AbstractVector{<:Distribution})
as a simplified interface to construct product distributions, e.g., to model a set of variables that share the same structure but vary by group.
Constructing product distributions with filldist
The function filldist
provides a general interface to construct product distributions over distributions of the same type and parameterisation.
Note that, in contrast to the product distribution interface provided by Distributions.jl (Product
), filldist
supports product distributions over univariate or multivariate distributions.
Example usage:
@model function demo(x, g)
k = length(unique(g))
a ~ filldist(Exponential(), k) # = Product(fill(Exponential(), k))
mu = a[g]
return x .~ Normal.(mu)
end
demo (generic function with 2 methods)
Constructing product distributions with arraydist
The function arraydist
provides a general interface to construct product distributions over distributions of varying type and parameterisation.
Note that in contrast to the product distribution interface provided by Distributions.jl (Product
), arraydist
supports product distributions over univariate or multivariate distributions.
Example usage:
@model function demo(x, g)
k = length(unique(g))
a ~ arraydist([Exponential(i) for i in 1:k])
mu = a[g]
return x .~ Normal.(mu)
end
demo (generic function with 2 methods)
Working with MCMCChains.jl
Turing.jl wraps its samples using MCMCChains.Chain
so that all the functions working for MCMCChains.Chain
can be re-used in Turing.jl. Two typical functions are MCMCChains.describe
and MCMCChains.plot
, which can be used as follows for an obtained chain chn
. For more information on MCMCChains
, please see the GitHub repository.
describe(chn) # Lists statistics of the samples.
plot(chn) # Plots statistics of the samples.
There are numerous functions in addition to describe
and plot
in the MCMCChains
package, such as those used in convergence diagnostics. For more information on the package, please see the GitHub repository.
Changing Default Settings
Some of Turing.jl's default settings can be changed for better usage.
AD Chunk Size
ForwardDiff (Turing's default AD backend) uses forward-mode chunk-wise AD. The chunk size can be set manually by setchunksize(new_chunk_size)
.
AD Backend
Turing supports four packages of automatic differentiation (AD) in the back end during sampling. The default AD backend is ForwardDiff for forward-mode AD. Three reverse-mode AD backends are also supported, namely Tracker, Zygote and ReverseDiff. Zygote
and ReverseDiff
are supported optionally if explicitly loaded by the user with using Zygote
or using ReverseDiff
next to using Turing
.
For more information on Turing's automatic differentiation backend, please see the Automatic Differentiation article.
Progress Logging
Turing.jl
uses ProgressLogging.jl to log the progress of sampling. Progress
logging is enabled as default but might slow down inference. It can be turned on
or off by setting the keyword argument progress
of sample
to true
or false
, respectively. Moreover, you can enable or disable progress logging globally by calling setprogress!(true)
or setprogress!(false)
, respectively.
Turing uses heuristics to select an appropriate visualization backend. If you use Juno, the progress is displayed with a progress bar in the Atom window. For Jupyter notebooks the default backend is ConsoleProgressMonitor.jl. In all other cases, progress logs are displayed with TerminalLoggers.jl. Alternatively, if you provide a custom visualization backend, Turing uses it instead of the default backend.