Turing2MonteCarloMeasurements

This package serves as an interface between Turing.jl and MonteCarloMeasurements.jl. Turing, as a probabilistic programming language and MCMC inference engine, produces results in the form of a Chain, a type that internally contains all the samples produced during inference. This chain is a bit awkward to work with in its natural form, why this package exists and allows for the conversion of a chain to a named tuple of Particles from MonteCarloMeasurements.jl.

Visualization

In this example, we simulate a review process where a number of reviewers are assigning scores to a number of articles. The generation of the data and the model specification are hidden under the collapsed section below.

We now focus on how to analyze the inference result. The chain is easily converted using the function Particles ```julia julia> chain = sample(reviewscore(R,nr,na), HMC(0.05, 10), 1500);

julia> cp = Particles(chain, crop=500); # crop discards the first 500 samples

julia> cp.reviewerpopbias Part1000(0.2605 ± 0.72)

julia> cp.reviewerpopgain Part1000(0.1831 ± 0.62) Particles can be plotted julia plot(cp.reviewerpopbias, title="Reviewer population bias") ![window](figs/rev_bias.svg) julia f1 = bar(articlescore, lab="Data", xlabel="Article number", ylabel="Article score", xticks=1:na) errorbarplot!(1:na, cp.truearticlescore, 0.8, seriestype=:scatter) f2 = bar(reviewerbias, lab="Data", xlabel="Reviewer number", ylabel="Reviewer bias") errorbarplot!(1:nr, cp.reviewer_bias, seriestype=:scatter, xticks=1:nr) plot(f1,f2) ```

Prediction

The linear-regression tutorial for Turing contains instructions on how to do prediction using the inference result. In the tutorial, the posterior mean of the parameters is used to form the prediction. Using Particles, we can instead form the prediction using the entire posterior distribution.

Like above, we hide the data generation under a collapsable section.

```julia @model linearregression(x, y, nobs, n_vars) = begin # Set variance prior. σ₂ ~ TruncatedNormal(0,100, 0, Inf)

# Set intercept prior.
intercept ~ Normal(0, 3)

# Set the priors on our coefficients.
coefficients = Array{Real}(undef, n_vars)
for i in 1:n_vars
    coefficients[i] ~ Normal(0, 10)
end

# Calculate all the mu terms.
mu = intercept .+ x * coefficients
for i = 1:n_obs
    y[i] ~ Normal(mu[i], σ₂)
end

end; nobs, nvars = size(x) model = linearregression(x, y, nobs, n_vars) chain = sample(model, NUTS(0.65), 2500); ```

In order to form the prediction, the original tutorial did julia function prediction(chain, x) p = get_params(chain[200:end, :, :]) α = mean(p.intercept) β = collect(mean.(p.coefficients)) return α .+ x * β end we will instead do julia cp = Particles(chain, crop=500) ŷ = x*cp.coefficients .+ cp.intercept plot(y, lab="data"); plot!(ŷ)

julia bar(coefficients, lab="True coeffs", title="Coefficients") errorbarplot!(1:n_vars, cp.coefficients, seriestype=:bar, alpha=0.5)

julia plot(plot.(cp.coefficients)..., legend=false) vline!(coefficients', l=(3,), lab="True value")

Further documentation

MonteCarloMeasurements

Turing