[avehtari]

This PR was assisted by Claude

Submission Checklist

• Run tests: ./runCmdStanTests.py src/test/interface/stansummary_test.cpp
• Declare copyright holder and open-source license: see below

Summary

This draft PR waits for the Stan submodule update that includes stan-dev/stan#3454. The first commit works with the current Stan submodule. The second commit updates test values that stan-dev/stan#3454 changes, so stansummary_test fails with the current submodule until that update.

Closes the remaining part of #903. MCSE, ESS, and R-hat are not meaningful for the sampler diagnostic columns accept_stat__, stepsize__, treedepth__, n_leapfrog__, divergent__, and energy__. #907 removed them in 2020; the current get_stats() computes them again for all columns. Users who check that all R-hat values are close to 1 see values for these columns. For discrete columns such as treedepth__ and divergent__, the output is a number or nan depending on whether one split chain is constant.

get_stats() has a new argument convergence_diagnostics (default true). stansummary sets it to false for the sampler diagnostics. These rows then show nan for MCSE, ESS_bulk, ESS_tail, ESS_bulk/s, and R_hat, in the console output and in the CSV file. Mean, StdDev, MAD, and the quantiles do not change. lp__ keeps all columns.

stan-dev/stan#3454 (fixes stan-dev/stan#3453) changes the ESS, R-hat, and MCSE functions that stansummary calls. Without this PR, the sampler rows would change with it, for example treedepth__ would show an ESS where it now shows nan.

Intended Effect

stansummary src/test/interface/example_output/bernoulli_chain_1.csv, sampler rows:

                Mean     MCSE   StdDev    MAD     5%   50%   95%  ESS_bulk  ESS_tail  ESS_bulk/s  R_hat

lp__            -7.3  3.7e-02     0.77   0.30   -9.0  -7.0  -6.8       519       503       22578    1.0
accept_stat__   0.90      nan  1.5e-01  0.064   0.57  0.96   1.0       nan       nan         nan    nan
stepsize__       1.0      nan  2.4e-15   0.00    1.0   1.0   1.0       nan       nan         nan    nan
treedepth__      1.4      nan  4.9e-01   0.00    1.0   1.0   2.0       nan       nan         nan    nan
n_leapfrog__     2.3      nan  9.6e-01   0.00    1.0   3.0   3.0       nan       nan         nan    nan
divergent__     0.00      nan  0.0e+00   0.00   0.00  0.00  0.00       nan       nan         nan    nan
energy__         7.8      nan  1.0e+00   0.75    6.8   7.5   9.9       nan       nan         nan    nan

How to Verify

./runCmdStanTests.py src/test/interface/stansummary_test.cpp. The expected accept_stat__ and energy__ rows in check_console_output, check_csv_output, check_csv_output_sig_figs, and check_csv_output_include_param now have nan in these columns.

The second commit changes the expected MCSE of lp__, theta, theta[6], and theta[7] in check_csv_output, check_csv_output_no_percentiles, and check_csv_output_include_param to the values of stan-dev/stan#NNNN, which computes mcse_mean on split chains as the R package posterior does.

Results (Stan Math 579253050, g++ 14.3.0, Linux x86_64):

• commit 1 with Stan develop c3db9b995: 21 of 21 tests pass.
• commit 2 with Stan fix-ess-constant-chain f324a5203: 21 of 21 tests pass.
• develop sources with the test file of commit 1: 4 tests (5 expectations) fail, check_console_output, check_csv_output, check_csv_output_sig_figs, and check_csv_output_include_param.
• develop sources with Stan fix-ess-constant-chain: 3 tests fail (check_csv_output, check_csv_output_no_percentiles, check_csv_output_include_param), because of the MCSE values that commit 2 updates.

Side Effects

The CSV file from --csv_filename has nan in these columns of the sampler rows. In the #907 discussion, CmdStanPy read this file for summary(); an interface that reads it gets nan there.

Documentation

If the stansummary examples in the CmdStan guide (stan-dev/docs) show values in these columns for the sampler rows, they need an update.

Copyright and Licensing

Please list the copyright holder for the work you are submitting (this will be you or your assignee, such as a university or company):

Aalto University

By submitting this pull request, the copyright holder is agreeing to license the submitted work under the following licenses:

• Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)
• Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)

[WardBrian]

#907 made them blank in the console output rather than nan, which might be clearer to users (I worry an nan value would also trigger concerns)

[bob-carpenter]

Why remove them? I've thought that before, but then realized reporting them makes sense because they're random variables. For example, reporting the mean and standard deviation of a random variable seems like MCSE on the mean would also make sense.

Technically, these stats can be coupled into the Markov chain by conditional generation based on position and momentum each iteration.

[avehtari]

stepsize__ is adapted separately for each chain, and during sampling fixed. Thus all chains have fixed but possibly different stepsize__, but not random. accept_stat__, stepsize__, treedepth__, n_leapfrog__, and divergent__ are random, but depend on that fixed stepsize__ per chain and if stepsize__ are different in chains, then they have different stationary distribution (e.g .small steps leads to higher acceptance, treedepth, and n_leapfrog) and reporting multihain Rhat, ESS_bulk, and ESS_tail doesn't make sense. Although stepsize__ would be different in different chains, the chains still have the same stationary target distribution with respect to parameters (and lp__).

[bob-carpenter]

Thanks, @avehtari.

If there's a natural place to add that last comment to the reference manual, that'd be very helpful to someone like me who will think through this just enough to be dangerous.

[avehtari]

I checked #907. Its description showed a sampler table without MCSE, ESS (then N_Eff) and R_hat, and the merged code kept these columns. Here is the output of CmdStan 2.34.1 and 2.40.0 for src/test/interface/example_output/bernoulli_chain_1.csv, up to energy__:

CmdStan 2.34.1, console:

                Mean     MCSE   StdDev     5%   50%   95%  N_Eff  N_Eff/s  R_hat

lp__            -7.3  3.7e-02     0.77   -9.1  -7.0  -6.8    443    19275    1.0
accept_stat__   0.90  4.6e-03  1.5e-01   0.57  0.96   1.0   1026    44597   1.00
stepsize__       1.0      nan  8.9e-16    1.0   1.0   1.0    nan      nan    nan
treedepth__      1.4  1.6e-02  4.9e-01    1.0   1.0   2.0    895    38919   1.00
n_leapfrog__     2.3  3.4e-02  9.6e-01    1.0   3.0   3.0    799    34758   1.00
divergent__     0.00      nan  0.0e+00   0.00  0.00  0.00    nan      nan    nan
energy__         7.8  5.1e-02  1.0e+00    6.8   7.5   9.9    411    17865    1.0

CmdStan 2.34.1, CSV:

name,Mean,MCSE,StdDev,5%,50%,95%,N_Eff,N_Eff/s,R_hat
"lp__",-7.2719,0.0365168,0.768874,-9.05757,-6.96978,-6.75008,443.328,19275.1,1.00037
"accept_stat__",0.896992,0.00460736,0.14756,0.568365,0.956737,1,1025.73,44597.1,0.999281
"stepsize__",1.02703,nan,8.88623e-16,1.02703,1.02703,1.02703,nan,nan,nan
"treedepth__",1.384,0.0162641,0.486601,1,1,2,895.135,38918.9,0.999067
"n_leapfrog__",2.31,0.0340858,0.963758,1,3,3,799.445,34758.5,0.999969
"divergent__",0,nan,0,0,0,0,nan,nan,nan
"energy__",7.78428,0.0508815,1.0314,6.80383,7.46839,9.88601,410.898,17865.1,1.00075

CmdStan 2.40.0, console:

                Mean     MCSE   StdDev    MAD     5%   50%   95%  ESS_bulk  ESS_tail  ESS_bulk/s  R_hat

lp__            -7.3  3.7e-02     0.77   0.30   -9.0  -7.0  -6.8       519       503       22578    1.0
accept_stat__   0.90  4.6e-03  1.5e-01  0.064   0.57  0.96   1.0      1284       941       55805   1.00
stepsize__       1.0      nan  2.4e-15   0.00    1.0   1.0   1.0       nan       nan         nan    nan
treedepth__      1.4  1.6e-02  4.9e-01   0.00    1.0   1.0   2.0       nan       nan         nan    nan
n_leapfrog__     2.3  3.4e-02  9.6e-01   0.00    1.0   3.0   3.0       791       787       34374   1.00
divergent__     0.00      nan  0.0e+00   0.00   0.00  0.00  0.00       nan       nan         nan    nan
energy__         7.8  5.1e-02  1.0e+00   0.75    6.8   7.5   9.9       490       486       21299    1.0

CmdStan 2.40.0, CSV:

name,Mean,MCSE,StdDev,MAD,5%,50%,95%,ESS_bulk,ESS_tail,ESS_bulk/s,R_hat
"lp__",-7.2719,0.0365168,0.768874,0.303688,-8.98426,-6.97009,-6.75007,519.29,503.309,22577.8,1.00141
"accept_stat__",0.896992,0.00460736,0.14756,0.0641484,0.569228,0.956732,1,1283.51,940.719,55804.6,0.99917
"stepsize__",1.02703,nan,2.44371e-15,0,1.02703,1.02703,1.02703,nan,nan,nan,nan
"treedepth__",1.384,0.0162641,0.486601,0,1,1,2,nan,nan,nan,nan
"n_leapfrog__",2.31,0.0340858,0.963758,0,1,3,3,790.601,786.719,34374,0.999969
"divergent__",0,nan,0,0,0,0,0,nan,nan,nan,nan
"energy__",7.78428,0.0508815,1.0314,0.745859,6.80565,7.46758,9.8864,489.874,486.438,21298.9,1.00495

Both versions already print nan in some sampler rows: stepsize__ and divergent__ in both, and treedepth__ in 2.40.0, where one split chain has a constant tree depth. With #1363, accept_stat__, stepsize__, treedepth__, n_leapfrog__ and divergent__ show nan in the MCSE, ESS and R_hat columns. I will change #1363 so that energy__ keeps these values, as lp__ does: like lp__, energy__ should have the same stationary distribution in all chains.

Would you prefer blank console cells for these rows? In the CSV file I would keep nan, because programs read that file and it already has nan in these rows.

[nhuurre]

For parameters that have a different stationary distribution for each chain, you can still compute meaningful single-chain summaries and combine them as if they were the components of a mixture distribution. From that perspective, the full ESS is the sum of single-chain ESSes; R-hat and MCSE are root-mean-squares of the single-chain R-hats and MCSEs.
Step-size is of course going to yield NaN regardless. (In my opinion, ESS for a constant chain should equal the number of draws but IIRC that was rejected.)

[avehtari]

. (In my opinion, ESS for a constant chain should equal the number of draws but IIRC that was rejected.)

It was rejected because in cases when the chain is constant due to being stuck by the failure of MCMC reporting ESS equal to the number of draws would be bad. At the point when we are computing ESS, we can't know whether the values are constant because 1) the quantity is truly structurally constant (ESS infinite make sense), 2) the chain is stuck because of failure (ESS 1 makes sense), or 3) a discretized quantity has low probability of visiting other discrete states and just hasn't yet visited them (e.g. when using indicator function x<=threshold, ESS not statistically well identified). Since we can't know at the point when we are computing which case holds, we decided to go with not available, and in R we would use NA, but as that is not available in Stan, we settled for nan.

For parameters that have a different stationary distribution for each chain, you can still compute meaningful single-chain summaries and combine them as if they were the components of a mixture distribution. From that perspective, the full ESS is the sum of single-chain ESSes; R-hat and MCSE are root-mean-squares of the single-chain R-hats and MCSEs.

I don't think that is relevant for our diagnostic quantities. For example, we would often prefer treedepth__ and n_leapfrog__ to have low variance as that would indicate low-varying curvature, and I really can't understand why would you care about ESS for these. Unless someone can demonstrate how a user could use ESS and Rhat of above mentioned sampler diagnostics for practical convergence diagnostic, I would write nan in csv and then display nan or empty (R style NA would be better than nan)

[nhuurre]

You write this PR "closes the remaining part of #903" but the reason those parts were left unclosed at the time is that Michael objected:

the evaluation of a fit they can be quite informative — large Rhats indicate
differing algorithmic behavior between the chains; low effective sample sizes
indicate large autocorrelations in algorithmic behavior that might arise,
for example, when the sampler transitions from one region of the posterior to
another.

My reading of this is that even though different chains theoretically converge to different distributions, in practice detectable differences often indicate poor sampling. Though he also wrote:

the sampler parameters (importantly the mean, variance, and quantiles) are informative of the fitting behavior [...] Rhat and ESS are incidental

so I guess he agrees that Rhat and ESS are best ignored here.

"sum of single chain ESS" is not interesting unless multiple chains are partially stuck in which case it correctly signals the problem. In contrast, NaN is the generic "calculation error" indicator but the indicated error is always just the theoretical non-stationarity.

---

when the chain is constant due to being stuck by the failure of MCMC reporting ESS equal to the number of draws would be bad.

Convergence is a global property. The important thing is that some diagnostic notices stuck chains; it doesn't matter much whether any individual ESS value looks problematic. The question should be whether a constant output is the only symptom of a stuck sampler.

When a chain is not moving it could be because:

• every transition diverges: diagnostics already warn
• very low accept_stat (but no divergence): no diagnostic but should be pretty obvious if you look
• extremely small step size: no specific diagnostic but probably hits max treedepth
• infinitesimal metric in some direction: no diagnostic, difficult to detect

Okay, that last one is quite nasty and closely corresponds to a constant parameter. I'd still prefer a more specific diagnostic.

1. the quantity is truly structurally constant (ESS infinite make sense)

Both MCMC variance and true variance are zero. Any ESS is justifiable.

2. the chain is stuck because of failure (ESS 1 makes sense)

MCSE is calculated as sqrt(sample_var / ess) so the ESS of a completely stuck chain (sample_var = 0) must be either 0 or NaN.

3. a discretized quantity has low probability of visiting other discrete states and just hasn't yet visited them

Any distribution, whether discrete or continuous, could have distant low probability modes that have not been visited yet. Such modes (if they exist) cause arbitrarily large bias in estimates of variance and ESS. That has never stopped us from estimating ESS; those estimates are just conditioned on the assumption of no unvisited modes.
Alternatively, one may employ a more handwavy definition of ESS like "inferences you make from this sample are about as good as inferences you could make with ESS-many independent draws from the true distribution." This definition only requires that the total probability mass of unvisited modes is less than 1/ESS and makes no claims about the total variance.

[avehtari]

The current output has columns ESS_bulk ESS_tail ESS_bulk/s R_hat. And I'm asking whether they should be shown as nan or empty. I'm also saying that specifically ESS_bulk ESS_tail ESS_bulk/s R_hat as they are defined are not sensible for the specific sampler diagnostic quantities. @nhuurre is arguing some other convergence diagnostics could be useful, bu then those should have their own columns. Also the earlier discussion did not show concrete proof for benefits in showing some alternative variant of ESS and Rhat for the sampler diagnostic. I'm happy to include such alternative variants for the sampler diagnostics, but I hope we are not again getting stalled because there might somehting but no-one actually demonstrates them. I suggest a new issue will be opened for discussion about diagnostics for sampler diagnostics, and the discussion here would focus on ESS_bulk ESS_tail ESS_bulk/s R_hat as they are currently defined.

[bob-carpenter]

I'm OK with blank or "n/a", though either one seems like a rendering nightmare for a table. I also don't think we should hold up the PR for final fine tuning of columns.