Interpreting pigeons' standard output
During the execution of parallel tempering, interim diagnostics can be computed and printed to standard out at the end of every iteration (this can be disabled using show_report = false
):
using Pigeons
pigeons(target = toy_mvn_target(100))
┌ Info: Neither traces, disk, nor online recorders included.
│ You may not have access to your samples (unless you are using a custom recorder, or maybe you just want log(Z)).
└ To add recorders, use e.g. pigeons(target = ..., record = [traces; record_default()])
────────────────────────────────────────────────────────────────────────────
scans Λ time(s) allc(B) log(Z₁/Z₀) min(α) mean(α)
────────── ────────── ────────── ────────── ────────── ────────── ──────────
2 7.5 2.81e-05 1.1e+04 -122 2.66e-15 0.167
4 5.56 4.58e-05 1.18e+04 -119 2.6e-07 0.382
8 5.82 6.66e-05 1.61e+04 -114 0.00137 0.353
16 7.02 0.000104 2.15e+04 -116 0.000205 0.219
32 7.14 0.000179 2.76e+04 -112 0.0379 0.206
64 7.18 0.000345 3.96e+04 -116 0.0442 0.202
128 7.2 0.000684 5.86e+04 -114 0.142 0.2
256 7.21 0.00131 5.66e+04 -115 0.151 0.199
512 7.14 0.00263 7.45e+04 -115 0.19 0.207
1.02e+03 7.12 0.00513 7.02e+04 -115 0.173 0.209
────────────────────────────────────────────────────────────────────────────
The functions called to emit each of these can be found at all_reports()
. Some key quantities:
Λ
: the global communication barrier, as described in Syed et al., 2021 and estimated using the sum of rejection estimator analyzed in the same reference. Syed et al., 2021 also developed a rule of thumb to configure the number of chains: PT should be set to roughly 2Λ.time
andallc
: the time (in second) and allocation (in bytes) used in each round.log(Z₁/Z₀)
: thestepping_stone()
estimator for the log of the normalization constant, see the documentation page on approximation of the normalization constant.min(α)
andmean(α)
: minimum and average swap acceptance rates over the PT chains.
Additional statistics can be shown when more recorders are added. For example, to accumulate other constant-memory summary statistics:
pigeons(target = toy_mvn_target(100), record = record_online(), explorer = AutoMALA())
───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
scans restarts Λ time(s) allc(B) log(Z₁/Z₀) min(α) mean(α) max|ρ| mean|ρ| min(αₑ) mean(αₑ)
────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
2 0 6.2 0.153 1.21e+07 -131 3.84e-13 0.311 1 1 0.284 0.536
4 0 6.21 0.00172 8.1e+04 -122 0.00023 0.31 0.827 0.345 0.41 0.539
8 0 7.19 0.0031 8.67e+04 -117 0.0114 0.201 0.617 0.361 0.441 0.582
16 0 7.14 0.00626 1.86e+05 -119 0.00252 0.207 0.586 0.268 0.459 0.523
32 0 7.1 0.0124 2.55e+05 -111 0.0483 0.211 0.63 0.24 0.451 0.511
64 0 7.09 0.0242 3.3e+05 -115 0.083 0.212 0.521 0.292 0.474 0.535
128 0 7.27 0.0477 3.82e+05 -115 0.112 0.192 0.462 0.303 0.468 0.533
256 1 7.34 0.0953 6.24e+05 -113 0.106 0.184 0.436 0.306 0.487 0.534
512 3 7.08 0.19 9.53e+05 -115 0.164 0.214 0.449 0.332 0.474 0.532
1.02e+03 11 7.18 0.379 1.57e+06 -115 0.161 0.202 0.472 0.332 0.48 0.534
───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
max|ρ|
andmean|ρ|
: maximum and average (across chains) correlation of the random variables $L^t_i = V(X_i)$ and $L^{t+1}_i = V(X_i)$ where $V = \log \pi_N / \pi_1$, $X_i \sim \pi_{\beta_i}$, and $t, t+1$ are indices just before and after a call tostep!()
.min(αₑ)
andmean(αₑ)
: minimum and average (across chains) of the explorer's acceptance rates.
Programmatic access
The tables described above can also be accessed as a DataFrame
via:
using Pigeons
pt = pigeons(target = toy_mvn_target(100))
pt.shared.reports.summary
Row | round | n_scans | n_tempered_restarts | global_barrier | global_barrier_variational | last_round_max_time | last_round_max_allocation | stepping_stone |
---|---|---|---|---|---|---|---|---|
Int64 | Int64 | Int64? | Float64? | Float64? | Float64? | Float64? | Float64? | |
1 | 1 | 2 | missing | 7.49502 | missing | 1.9717e-5 | 11008.0 | -122.071 |
2 | 2 | 4 | missing | 5.5631 | missing | 3.1239e-5 | 11776.0 | -118.917 |
3 | 3 | 8 | missing | 5.82286 | missing | 6.7015e-5 | 16128.0 | -114.195 |
4 | 4 | 16 | missing | 7.02455 | missing | 0.000108002 | 21504.0 | -116.005 |
5 | 5 | 32 | missing | 7.1422 | missing | 0.000190537 | 27552.0 | -111.707 |
6 | 6 | 64 | missing | 7.18489 | missing | 0.00035806 | 39632.0 | -115.937 |
7 | 7 | 128 | missing | 7.202 | missing | 0.00071067 | 58560.0 | -114.282 |
8 | 8 | 256 | missing | 7.21158 | missing | 0.00133582 | 56608.0 | -114.816 |
9 | 9 | 512 | missing | 7.13659 | missing | 0.00258838 | 74528.0 | -115.087 |
10 | 10 | 1024 | missing | 7.123 | missing | 0.00510923 | 70224.0 | -115.257 |
Detailed statistics can be accessed via these DataFrames. For example to obtain mean swap acceptance probabilities for each round and pair of communicating chains, use:
first(pt.shared.reports.swap_prs, 20)
Row | round | first | second | mean |
---|---|---|---|---|
Int64 | Int64 | Int64 | Float64 | |
1 | 1 | 1 | 2 | 2.66334e-15 |
2 | 1 | 2 | 3 | 1.42777e-6 |
3 | 1 | 3 | 4 | 8.3416e-5 |
4 | 1 | 4 | 5 | 0.00380696 |
5 | 1 | 5 | 6 | 0.246034 |
6 | 1 | 6 | 7 | 0.233044 |
7 | 1 | 7 | 8 | 0.178188 |
8 | 1 | 8 | 9 | 0.396346 |
9 | 1 | 9 | 10 | 0.447479 |
10 | 2 | 1 | 2 | 2.59743e-7 |
11 | 2 | 2 | 3 | 5.68511e-6 |
12 | 2 | 3 | 4 | 0.0113721 |
13 | 2 | 4 | 5 | 0.323006 |
14 | 2 | 7 | 8 | 0.484048 |
15 | 2 | 6 | 7 | 1.0 |
16 | 2 | 5 | 6 | 0.462578 |
17 | 2 | 9 | 10 | 0.550769 |
18 | 2 | 8 | 9 | 0.60512 |
19 | 3 | 1 | 2 | 0.00137441 |
20 | 3 | 2 | 3 | 0.254306 |
Creating a plot from this:
using StatsPlots
plotlyjs()
my_plot = @df pt.shared.reports.swap_prs StatsPlots.plot(:round, :mean, group = :first)
StatsPlots.savefig(my_plot, "swap_prs.html");