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.7239    0.1520    0.0065   574.1139   703.2334    0.9999     ⋯
          p2    0.7086    0.1507    0.0062   614.1593   920.7039    0.9998     ⋯
                                                                1 column omitted

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

          p1    0.4735    0.5934    0.7108    0.8487    0.9874
          p2    0.4647    0.5850    0.6959    0.8270    0.9815

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.7239
          p2    0.7086

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

m[:p1, :mean]

0.723864629673392

Highest posterior density interval:

hpd(samples, alpha = 0.05)

HPD
  parameters     lower     upper 
      Symbol   Float64   Float64 

          p1    0.5003    0.9980
          p2    0.4803    0.9919

For ESS estimates:

ess(samples)

ESS
  parameters        ess   ess_per_sec 
      Symbol    Float64       Missing 

          p1   574.1139       missing
          p2   614.1593       missing