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.35e-05   1.08e+04       -122   2.66e-15      0.167
        4       5.59   3.99e-05   1.21e+04       -119    2.6e-07      0.378
        8       6.04   5.88e-05    1.6e+04       -115    0.00129      0.329
       16       7.27   0.000104   2.16e+04       -118     0.0134      0.193
       32       6.97   0.000192    2.8e+04       -114      0.107      0.225
       64       7.03   0.000368   3.88e+04       -117     0.0531      0.219
      128       7.23   0.000741   5.79e+04       -114     0.0944      0.196
      256       7.05     0.0014   6.31e+04       -115       0.13      0.217
      512       7.14     0.0027   6.66e+04       -115      0.171      0.207
 1.02e+03       7.19    0.00537   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 and allc: the time (in second) and allocation (in bytes) used in each round.
  • log(Z₁/Z₀): the stepping_stone() estimator for the log of the normalization constant, see the documentation page on approximation of the normalization constant.
  • min(α) and mean(α): 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.124   9.06e+06       -131   3.85e-13      0.311          1          1      0.284      0.536
        4          0       6.18    0.00173    7.9e+04       -123   7.56e-05      0.314      0.827       0.37       0.41      0.535
        8          0       6.91    0.00308   8.76e+04       -115     0.0183      0.232      0.637      0.502      0.438      0.567
       16          0       7.14    0.00611   1.42e+05       -114     0.0311      0.207      0.547      0.296      0.476      0.537
       32          0       6.66      0.012   1.85e+05       -116     0.0441       0.26      0.514      0.279      0.485      0.539
       64          0       6.79     0.0239   2.74e+05       -116     0.0963      0.246      0.617      0.334      0.463      0.535
      128          0       7.27     0.0474      4e+05       -115      0.126      0.192      0.357      0.295      0.451      0.533
      256          2       7.22      0.095   6.82e+05       -114      0.136      0.198      0.444      0.307      0.475      0.533
      512          4       7.19       0.19   9.08e+05       -116       0.15      0.202      0.408      0.302      0.485      0.531
 1.02e+03         12       7.21      0.378   1.48e+06       -115      0.167      0.199      0.428       0.33      0.483      0.531
───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
  • max|ρ| and mean|ρ|: 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 to step!().
  • min(αₑ) and mean(αₑ): 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
10×8 DataFrame
Rowroundn_scansn_tempered_restartsglobal_barrierglobal_barrier_variationallast_round_max_timelast_round_max_allocationstepping_stone
Int64Int64Int64?Float64?Float64?Float64?Float64?Float64?
112missing7.49502missing2.2322e-510784.0-122.071
224missing5.5942missing3.6719e-512096.0-119.001
338missing6.03849missing6.0943e-516000.0-114.872
4416missing7.26595missing0.00010403421600.0-118.087
5532missing6.97496missing0.00019327127952.0-113.96
6664missing7.03149missing0.00036868838768.0-116.851
77128missing7.23408missing0.00074587357920.0-113.569
88256missing7.0475missing0.0013897463104.0-115.08
99512missing7.13997missing0.002744266560.0-114.586
10101024missing7.18702missing0.0053612573216.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)
20×4 DataFrame
Rowroundfirstsecondmean
Int64Int64Int64Float64
11122.66334e-15
21231.42777e-6
31348.3416e-5
41450.00380696
51560.246034
61670.233044
71780.178188
81890.396346
919100.447479
102122.59743e-7
112235.68511e-6
122340.0113721
132450.323006
142780.484048
152671.0
162560.462578
1729100.519665
182890.60512
193120.00129063
203230.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");