[avehtari]

Submission Checklist

• Run unit tests: ./runTests.py src/test/unit/analyze/mcmc, src/test/unit/mcmc/chainset_test.cpp, src/test/unit/mcmc/chains_test.cpp
• Run cpplint: make cpplint
• Declare copyright holder and open-source license: see below

Summary

Fixes #3453. Four commits; the first fixes the issue, the others remove the remaining differences from the R package posterior. Each commit has its own tests.

Commit 1: constant chains (#3453)

Current behavior (develop): autocorrelation() divides by the lag-0 autocovariance, which is 0 for a constant chain. The NaN reaches ess(), the guard !std::isnan(rho_hat_even + rho_hat_odd) stops the Geyer loop at once, and ess() returns exactly N/2. This happens for every matrix that has a constant column, and for a constant matrix. In split_rank_normalized_ess(), the tail indicator of a split chain that never enters a tail is constant, so the tail-ESS of chains in different modes is N/2. Separately, is_finite_and_varies() returns false if one (split) chain is constant, so R-hat and ESS are NaN when one chain does not move.

New behavior:

• autocovariance() returns 0 at every lag for a constant sequence, as posterior::autocovariance() does. A constant chain then contributes zero within-chain variance and zero autocovariance, and ess() gives a finite value.
• ess() returns NaN if all draws are equal or a draw is not finite (var_plus is not positive).
• is_finite_and_varies() returns false only if a value is not finite or every chain is constant. posterior returns NA only if all draws are equal. This commit also returns NaN when every chain is constant at different values. That case includes stepsize__, which is constant in each chain and different between chains, so stansummary output for stepsize__ does not change.

Commit 2: end of the Geyer sum

Current behavior (develop): the loop over pairs of autocorrelations stops at the first negative pair or at the lag limit N - 4. At a negative pair, ess() does not store the pair, and the result agrees with posterior to 1e-14. At the lag limit the last pair is positive and stored. rho_hat_t.head(max_t) includes its even lag, and rho_hat_t(max_t + 1) adds that lag again, so it counts three times instead of once. The loop reaches the lag limit when the autocorrelations stay positive, which happens for chains that do not mix and for short chains. The ESS then differs from posterior (it is lower in 1342 of 1847 random cases), and the relative difference is about 1/N (median over random AR(1) inputs: 0.13 at N = 10, 0.020 at N = 50, 0.0018 at N = 500).

New behavior: the truncation lag is the even lag of the last pair computed, both at a negative pair and at the lag limit, as in posterior. At the lag limit that lag counts once. The result at a negative pair does not change.

Commit 3: tail-ESS when one tail is NaN

Current behavior (develop): the tail indicators are x <= q05 and x >= q95. If more than 95% of the draws are at one value, the indicator of one tail is constant and its ESS is NaN. split_rank_normalized_ess() then reports the ESS of the other tail as the tail-ESS.

New behavior: the tail-ESS is NaN if the ESS of either tail is NaN, as posterior::ess_tail() does. The indicators do not change.

In addition, stan-dev/posterior#NNNN changes the upper-tail indicator of posterior from x <= q95 to x >= q95, as in Stan. Both then give the same tail-ESS for discrete x and -x, and a finite tail-ESS for discrete draws with more than 5% at the maximum.

Commit 4: MCSE

Current behavior (develop): mcse_mean() uses the ESS of the unsplit chains, while posterior (ess_mean()) and every other diagnostic in this namespace use split chains. mcse_sd() uses the ESS of the split chains of |x - mean|, as an old posterior version did, and computes the variance of the variance as mean(d^4) - mean(d^2)^2.

New behavior: mcse_mean() uses ess(split_chains(x)). mcse_sd() uses the ESS of the split chains of d2 = (x - mean)^2, the quantity in the variance estimate, as posterior now does. It computes the variance of the variance as mean((d2 - mean(d2))^2), which avoids a cancellation.

Intended Effect

Example of #3453, 4 chains of 1000 draws, chains 1 to 3 iid N(0, 1):

chain 4	quantity	develop	this PR	posterior
N(6, 1)	tail-ESS	2000	28.25	28.25
N(6, 1)	bulk-ESS	7.21	7.23	7.23
constant 0.3	R-hat	NaN	1.526	1.526
constant 0.3	tail-ESS	NaN	818.9	818.9
2nd half constant	R-hat	NaN	1.122	1.122
N(6, 1)	MCSE of mean	1.865	1.304	1.304

After the four commits, R-hat, bulk-ESS, tail-ESS and both MCSEs equal the values of posterior (master plus stan-dev/posterior#NNNN) to the printed digits for all test inputs, except when every chain is constant (NaN here, Inf R-hat in posterior).

How to Verify

./runTests.py src/test/unit/analyze/mcmc and ./runTests.py src/test/unit/mcmc. The reference values in the new tests are from posterior. The inputs are the bern*.csv theta draws with a change to chain 4: constant, second half constant, or shifted by 1 so that it never enters either 5% tail.

• Commit 1: new tests in autocovariance_test, check_chains_test, ess_basic_test, split_rank_normalized_ess_test and split_rank_normalized_rhat_test. With the develop headers, 9 tests (32 expectations) fail. Three ESS values at the lag limit have a tolerance of 1% in this commit.
• Commit 2: the tolerance becomes 1e-8, and a new test has autocorrelations that stay positive to the lag limit. With the headers of commit 1, 4 tests (5 expectations) fail.
• Commit 3: a new test with binary draws (46% and 96% ones, and the mirror image). With the headers of commit 2, 1 test (2 expectations) fails.
• Commit 4: mcse_test reference values from posterior with a tolerance of 1e-10 (it was 1e-4), and a constant-chain test. With the headers of commit 3, 2 tests (6 expectations) fail.

Every commit passes all of src/test/unit/analyze/mcmc, src/test/unit/mcmc/chainset_test and src/test/unit/mcmc/chains_test (13 suites, 86 tests after commit 4). Tested with develop c3db9b995 and Stan Math 579253050, g++ 14.3.0, Linux x86_64.

Side Effects

CmdStan stansummary output changes:

• R-hat, ESS and MCSE are no longer NaN for a quantity with some constant chains, for example a parameter with one stuck chain.
• ESS and MCSE change slightly for chains that do not mix (commit 2) and for all quantities (MCSE, commit 4).
• Tail-ESS is NaN for a quantity with more than 95% of the draws at one value (commit 3).
• The sampler diagnostic columns would change too. stan-dev/cmdstan#NNNN stops computing MCSE, ESS and R-hat for those columns (stansummary: separate sections for summaries of sampler parameters and model parameters cmdstan#903). Its second commit updates the stansummary_test values that this PR changes: the MCSE of lp__, theta, theta[6] and theta[7] in check_csv_output, check_csv_output_no_percentiles and check_csv_output_include_param. With CmdStan develop, these 3 tests fail with this PR; with that CmdStan branch, all 21 stansummary_test tests pass.

The deprecated compute_effective_sample_size() and compute_split_effective_sample_size() also call autocovariance(), so their results change for inputs with a constant chain; their tests pass unchanged.

Documentation

Doc comments of ess(), autocovariance(), is_finite_and_varies(), mcse_mean() and mcse_sd().

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]

@avehtari can you confirm that the failures in cmdstan are expected? The existence of stan-dev/cmdstan#1363 seems to imply so, but wanted to be sure

[avehtari]

Yes. With this PR, 3 of the 21 tests in src/test/interface/stansummary_test.cpp fail in CmdStan develop: check_csv_output, check_csv_output_no_percentiles and check_csv_output_include_param. These tests compare whole CSV lines. mcse_mean now uses split chains, which changes the MCSE values in the 3rd or 4th significant digit. The second commit of stan-dev/cmdstan#1363 updates the expected lines.

CSV output (--csv_filename) for src/test/interface/example_output/bernoulli_chain_1.csv, CmdStan develop with Stan develop:

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
"theta",0.256552,0.00610844,0.119654,0.120965,0.0802982,0.24996,0.47034,361.506,395.736,15717.6,1.00186

CmdStan develop with this PR:

name,Mean,MCSE,StdDev,MAD,5%,50%,95%,ESS_bulk,ESS_tail,ESS_bulk/s,R_hat
"lp__",-7.2719,0.0365073,0.768874,0.303688,-8.98426,-6.97009,-6.75007,519.29,503.309,22577.8,1.00141
"accept_stat__",0.896992,0.00460677,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.0162008,0.486601,0,1,1,2,902.141,902.141,39223.5,0.999067
"n_leapfrog__",2.31,0.0342761,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.0507653,1.0314,0.745859,6.80565,7.46758,9.8864,489.874,486.438,21298.9,1.00495
"theta",0.256552,0.00609368,0.119654,0.120965,0.0802982,0.24996,0.47034,361.506,395.736,15717.6,1.00186

For example, the MCSE of lp__ changes from 0.0365168 to 0.0365073 here, and from 0.97319 to 0.952577 in eight_schools_output.csv, which check_csv_output_include_param uses.

The console output shows 2 significant digits, so check_console_output passes. The only change there is the treedepth__ row. One split chain has a constant tree depth, which gave nan before this PR; now ESS and R_hat are numbers. stan-dev/cmdstan#1363 stops reporting these diagnostics for the sampler rows.

Console output, CmdStan develop with Stan develop:

                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

theta           0.26  6.1e-03     0.12   0.12  0.080  0.25  0.47       362       396       15718    1.0

CmdStan develop with this PR:

                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       902       902       39224   1.00
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

theta           0.26  6.1e-03     0.12   0.12  0.080  0.25  0.47       362       396       15718    1.0

[WardBrian]

The code here looks good but I ask that you don't merge until the cmdstan PR is also approved to prevent there being a period of time where CI is broken. Thanks!