We’re happy to announce Turing 0.13.0, which brings a suite of much-needed features. I am particularly happy about this one because it is a feature-rich release, which is always fun.
Maximum likelihood and maximum a posteriori estimation
Turing’s model flexibility doesn’t necessarily restrict its use to only Bayesian methods – traditional maximum likelihood (MLE) or maximum a posteriori (MAP) estimation are supportable as well. If you work in a big-data style field, it’s now push-button easy to generate MAP or MLE estimates for your model, as long as there are no discrete parameters.
We recently added an interface to the world-class Optim.jl to handle this. Here is a brief example using the classic gdemo model:
using Turing
using Optim
using StatsBase
@model function gdemo()
s ~ InverseGamma(2,3)
m ~ Normal(0, sqrt(s))
1.5 ~ Normal(m, sqrt(s))
2.0 ~ Normal(m, sqrt(s))
return m, s
end
model = gdemo()
# Use the default optimizer, LBFGS.
mle1 = optimize(model, MLE())
mpa1 = optimize(model, MAP())
# Use a specific optimizer.
mle2 = optimize(model, MLE(), Optim.Newton())
mpa2 = optimize(model, MAP(), Optim.Newton())
We extended some of the StatsBase.jl methods, like coeftable, informationmatrix, stderror, vcov, etc. For example, we can examine whether our MAP estimates are statistically significant using standard errors calculated from the information matrix with coeftable:
julia> coeftable(mle1)
──────────────────────────────
estimate stderror tstat
──────────────────────────────
s 0.0625 0.0625 1.0
m 1.75 0.176777 9.89949
──────────────────────────────
or for the MAP estimate:
julia> coeftable(mpa1)
──────────────────────────────
estimate stderror tstat
──────────────────────────────
s 0.907407 0.427756 2.12132
m 1.16667 0.549972 2.12132
──────────────────────────────
You can also initialize sampling with the MAP or MLE estimates:
chain = sample(model, NUTS(), 1000, init_theta = mpa1.values.array)
We’re looking forward to seeing what people do with these new features. As always, this is new, so please report any issues you find.
Prediction
Historically, Turing has not natively supported posterior prediction, and many users had to roll their own code. No longer! Turing now implements the StatsBase.predict function that allows you to predict your data conditional on a previous Markov chain. Make sure you load StatsBase or call StatsBase.predict.
Here’s an example:
using Turing
using StatsBase
@model function linear_reg(x, y, σ = 0.1)
β ~ Normal(0, 1)
for i ∈ eachindex(y)
y[i] ~ Normal(β * x[i], σ)
end
end
f(x) = 2 * x + 0.1 * randn()
Δ = 0.1
xs_train = 0:Δ:10; ys_train = f.(xs_train);
xs_test = [10 + Δ, 10 + 2 * Δ]; ys_test = f.(xs_test);
# Infer
m_lin_reg = linear_reg(xs_train, ys_train);
chain_lin_reg = sample(m_lin_reg, NUTS(100, 0.65), 200);
# Predict on two last indices by making an uninitialized vector of length(ys_test)
# Make a new version of the model -- note that we put in our training data instead.
m_lin_reg_test = linear_reg(xs_test, Vector{Union{Missing, Float64}}(undef, length(ys_test)));
# Use the new predict function!
predictions = predict(m_lin_reg_test, chain_lin_reg)
# Get the mean predicted values.
ys_pred = collect(vec(mean(predictions[:y].value; dims = 1)))
# Get the prediction error:
ys_test - ys_pred
Forecast error:
2-element Array{Float64,1}:
0.057563331644999494
-0.07099319663700498
Prior sampling
Prior predictive checks are now much easier. You can sample directly from your prior with
chain = sample(model, Prior(), n_samples)
Not too exciting, but very necessary for good Bayesian analysis.