Plotting

Use sample_array() to convert target chain samples into a format that can then be consumed by third party plotting packages such as MCMCChains.jl and PairPlots.jl.

See below for examples of posterior densities and trace plots.

Posterior densities and trace plots

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

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

Then use the following:

using DynamicPPL
using Pigeons
using MCMCChains
using StatsPlots
plotlyjs()

# 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,
                n_rounds = 12,
                # make sure to record the trace:
                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(sample_array(pt), sample_names(pt))

# since the above line is frequently needed, Pigeons includes
# an MCMCChains extension allowinging you to use the shorter form:
samples = Chains(pt)

# create the trace plots
my_plot = StatsPlots.plot(samples)
StatsPlots.savefig(my_plot, "posterior_densities_and_traces.html");

─────────────────────────────────────────────────────────────────────────────────────────────────────────────
  scans     restarts      Λ        time(s)    allc(B)  log(Z₁/Z₀)   min(α)     mean(α)    min(αₑ)   mean(αₑ)
────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
        2          0       3.24    0.00118   7.03e+05      -8.14    0.00178       0.64          1          1
        4          0       1.64    0.00219   1.38e+06      -5.04     0.0352      0.818          1          1
        8          0       1.17    0.00419   2.69e+06      -4.42      0.708      0.871          1          1
       16          1        1.2    0.00866   5.65e+06      -4.03      0.549      0.867          1          1
       32          6       1.11     0.0167    1.1e+07      -4.77      0.754      0.877          1          1
       64         11       1.35     0.0333   2.21e+07      -4.79      0.698       0.85          1          1
      128         25        1.6     0.0651   4.37e+07      -4.97      0.725      0.823          1          1
      256         43       1.51      0.129   8.64e+07      -4.92      0.758      0.832          1          1
      512         92       1.46      0.324   1.74e+08         -5      0.806      0.838          1          1
 1.02e+03        188       1.49      0.582   3.49e+08      -4.92      0.798      0.834          1          1
 2.05e+03        384        1.5       1.16   6.98e+08      -4.96      0.811      0.834          1          1
  4.1e+03        748       1.48       2.34   1.39e+09      -4.94      0.826      0.835          1          1
─────────────────────────────────────────────────────────────────────────────────────────────────────────────

Monitoring the log density

The value of the log density is appended to each sample. Continuing the above example, this can be seen from the variable names indexing the flattened vector created by sample_array():

sample_names(pt)

3-element Vector{Symbol}:
 :p1
 :p2
 :log_density

When using the Chains(pt) constructor as shown above, the un-normalized log density is stored inside MCMCChains' "internal" storage so will not appear in plots by default. To show it, use the following:

params, internals = MCMCChains.get_sections(samples)

my_plot = StatsPlots.plot(internals)
StatsPlots.savefig(my_plot, "logdensity.html");

Posterior pair plots

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

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

using DynamicPPL
using Pigeons
using MCMCChains
using CairoMakie
using PairPlots

# same examples as last section
an_unidentifiable_model = Pigeons.toy_turing_unid_target(100, 50)

pt = pigeons(target = an_unidentifiable_model, 
                n_rounds = 12,
                # make sure to record the trace:
                record = [traces; round_trip; record_default()])

samples = Chains(pt)
# Warning: the line below only works for Julia 1.9
#          see https://sefffal.github.io/PairPlots.jl/dev/chains/ for a workaround
my_plot = PairPlots.pairplot(samples) 

CairoMakie.save("pair_plot.svg", my_plot)
nothing # hide