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.24e-05 1.08e+04 -122 2.66e-15 0.167
4 5.59 3.9e-05 1.21e+04 -119 2.6e-07 0.378
8 6.04 5.98e-05 1.6e+04 -115 0.00129 0.329
16 7.27 0.000105 2.16e+04 -118 0.0134 0.193
32 6.97 0.000195 2.8e+04 -114 0.107 0.225
64 7.03 0.000375 3.88e+04 -117 0.0531 0.219
128 7.23 0.00073 5.79e+04 -114 0.0944 0.196
256 7.05 0.00141 6.31e+04 -115 0.13 0.217
512 7.14 0.00278 6.66e+04 -115 0.171 0.207
1.02e+03 7.19 0.00548 7.32e+04 -115 0.172 0.201
────────────────────────────────────────────────────────────────────────────
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.126 9.06e+06 -131 3.85e-13 0.311 1 1 0.284 0.536
4 0 6.18 0.00175 7.9e+04 -123 7.56e-05 0.314 0.827 0.37 0.41 0.535
8 0 6.91 0.00314 8.76e+04 -115 0.0183 0.232 0.637 0.502 0.438 0.567
16 0 7.14 0.00623 1.42e+05 -114 0.0311 0.207 0.547 0.296 0.476 0.537
32 0 6.66 0.0122 1.85e+05 -116 0.0441 0.26 0.514 0.279 0.485 0.539
64 0 6.79 0.0242 2.74e+05 -116 0.0963 0.246 0.617 0.334 0.463 0.535
128 0 7.27 0.0482 4e+05 -115 0.126 0.192 0.357 0.295 0.451 0.533
256 2 7.22 0.096 6.82e+05 -114 0.136 0.198 0.444 0.307 0.475 0.533
512 4 7.19 0.191 9.08e+05 -116 0.15 0.202 0.408 0.302 0.485 0.531
1.02e+03 12 7.21 0.381 1.48e+06 -115 0.167 0.199 0.428 0.33 0.483 0.531
───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
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 | 2.2793e-5 | 10784.0 | -122.071 |
2 | 2 | 4 | missing | 5.5942 | missing | 3.4745e-5 | 12096.0 | -119.001 |
3 | 3 | 8 | missing | 6.03849 | missing | 5.6666e-5 | 16000.0 | -114.872 |
4 | 4 | 16 | missing | 7.26595 | missing | 0.00010157 | 21600.0 | -118.087 |
5 | 5 | 32 | missing | 6.97496 | missing | 0.000190035 | 27952.0 | -113.96 |
6 | 6 | 64 | missing | 7.03149 | missing | 0.000361094 | 38768.0 | -116.851 |
7 | 7 | 128 | missing | 7.23408 | missing | 0.000700948 | 57920.0 | -113.569 |
8 | 8 | 256 | missing | 7.0475 | missing | 0.00136713 | 63104.0 | -115.08 |
9 | 9 | 512 | missing | 7.13997 | missing | 0.0026589 | 66560.0 | -114.586 |
10 | 10 | 1024 | missing | 7.18702 | missing | 0.00527699 | 73216.0 | -115.411 |
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.519665 |
18 | 2 | 8 | 9 | 0.60512 |
19 | 3 | 1 | 2 | 0.00129063 |
20 | 3 | 2 | 3 | 0.254206 |
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");