Numerical outputs and diagnostics

Use sample_array() to convert target chain samples into a format that can then be consumed by the third party package MCMCChains.jl. We outline some useful features here, read the MCMCChains.jl documentation for more information.

Quick summary of ESS, moments, etc

Make sure to have the third party packages DynamicPPL and MCMCChains installed via

using Pkg; Pkg.add("DynamicPPL", "MCMCChains")

Also make sure to record the trace, with record = [traces]:

using DynamicPPL
using Pigeons
using MCMCChains

# example target: Binomial likelihood with parameter p = p1 * p2
an_unidentifiable_model = Pigeons.toy_turing_unid_target(100, 50)

pt = pigeons(
        target = an_unidentifiable_model,
        record = [traces; round_trip; record_default()]
    )

# collect the statistics and convert to MCMCChains' Chains
# to have axes labels matching variable names in Turing and Stan
samples = Chains(pt)

samples

Chains MCMC chain (1024×3×1 Array{Float64, 3}):

Iterations        = 1:1:1024
Number of chains  = 1
Samples per chain = 1024
parameters        = p1, p2
internals         = log_density

Summary Statistics
  parameters      mean       std      mcse   ess_bulk   ess_tail      rhat   e ⋯
      Symbol   Float64   Float64   Float64    Float64    Float64   Float64     ⋯

          p1    0.7152    0.1518    0.0068   491.5895   649.6713    1.0026     ⋯
          p2    0.7185    0.1507    0.0067   525.7289   783.7176    1.0040     ⋯
                                                                1 column omitted

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5% 
      Symbol   Float64   Float64   Float64   Float64   Float64 

          p1    0.4533    0.5884    0.7068    0.8417    0.9804
          p2    0.4729    0.5932    0.7008    0.8467    0.9856

Accessing individual diagnostics and summaries

Computing a mean (but see online statistics for a constant memory alternative):

using Statistics
m = mean(samples)

Mean
  parameters      mean 
      Symbol   Float64 

          p1    0.7152
          p2    0.7185

to access an individual entry in this example and the following ones:

m[:p1, :mean]

0.715241363315093

Highest posterior density interval:

hpd(samples, alpha = 0.05)

HPD
  parameters     lower     upper 
      Symbol   Float64   Float64 

          p1    0.4729    0.9897
          p2    0.4945    0.9989

For ESS estimates:

ess(samples)

ESS
  parameters        ess   ess_per_sec 
      Symbol    Float64       Missing 

          p1   491.5895       missing
          p2   525.7289       missing