Use posterior package for convergence diagnostics

[avehtari]

Currently, convergence diagnostic functions are defined in monitor.R. These are mostly the same as what is in the posterior package, but there are some updates in posterior package that have not been added to monitor.R. To reduce amount of duplicate code and to unify the diagnostics between rstan and cmdstanr, it would be good to switch to use the functions directly from posterior package

[VisruthSK]

Should the posterior functions be wrapped? Or is it ok to completely remove some of the currently exported implementations and just use the posterior functions in monitor.R.

e.g. should

#' @export
ess_bulk <- function(sims) {
  ess_rfun(z_scale(split_chains(sims)))
}

go to

#' @export
ess_bulk <- function(sims) {
  posterior::ess_bulk(sims)
}

or can it be deleted completely? I'm replacing every usage internally with posterior::ess_bulk() anyway.

[jgabry]

I think if it's currently exported then we should still export it and not delete it (just replace internals). For things that aren't exported we can probably get rid of them if they're not needed anymore.

[VisruthSK]

I've been running into issues running tests, could you see if any fail? I somehow got the test suite to run once but it skipped a bunch of test.

Error:
! ! System command 'Rcmd.exe' failed

[jgabry]

I didn't get any test failures, but I think for some reason monitor is not tested, I'm not sure why. It may be tested elsewhere. Anyway, when I try running monitor on your branch, I get an error:

Error: 'quantile' is not an exported object from 'namespace:posterior'

I think it's coming from here:

quan2 <- posterior::quantile(sims_i, probs = c(0.05, 0.5, 0.95))

because posterior exports quantile2 but not quantile.

I guess for testing, I would just run some example models and call monitor on them and make sure there's no error. It's probably enough to just use the two models in the Examples section of help("stan", package = "rstan"), that's what I tried and got the error.

[VisruthSK]

Oop I totally missed that, thanks. Cool, I'll run those examples and fix anything.

[jgabry]

Actually, just use this example model, it should be enough.

Here's an example model with monitor results using the CRAN version of rstan. You can compare these to the results of using monitor on the same model but using your branch to make sure the results aren't meaningfully different (they may be very slightly different due to Monte Carlo error because even setting the same seed doesn't produce 100% identical results on with hardware/compiler/etc.):

excode <- '
  transformed data {
    array[20] real y;
    y[1] = 0.5796;  y[2] = 0.2276;   y[3]  = -0.2959;
    y[4] = -0.3742; y[5] = 0.3885;   y[6]  = -2.1585;
    y[7] = 0.7111;  y[8] = 1.4424;   y[9]  = 2.5430;
    y[10] = 0.3746; y[11] = 0.4773;  y[12] = 0.1803;
    y[13] = 0.5215; y[14] = -1.6044; y[15] = -0.6703;
    y[16] = 0.9459; y[17] = -0.382;  y[18] = 0.7619;
    y[19] = 0.1006; y[20] = -1.7461;
  }
  parameters {
    real mu;
    real<lower=0, upper=10> sigma;
    matrix[3,2] z;
    real<lower=0> alpha;
  }
  model {
    y ~ normal(mu, sigma);
    for (i in 1:3)
      z[i] ~ normal(0, 1);
    alpha ~ exponential(2);
  }
'

# you can run for whatever number of iterations you want, 
# just make sure it's the same between the CRAN version and your version. 
# My results below are with iter = 500
exfit <- stan(model_code = excode, iter = 500)
monitor(exfit)

results using monitor() from CRAN

          Q5   Q50   Q95  Mean  SD  Rhat Bulk_ESS Tail_ESS
mu      -0.4   0.1   0.5   0.1 0.3  1.01     1277      689
sigma    0.9   1.1   1.6   1.2 0.2  1.00     1551      765
z[1,1]  -1.7   0.0   1.6   0.0 1.0  1.01     1626      855
z[2,1]  -1.7   0.1   1.8   0.1 1.0  1.00     1397      643
z[3,1]  -1.7   0.0   1.7   0.0 1.0  1.01     1737      756
z[1,2]  -1.7   0.0   1.7   0.0 1.0  1.00     1864      689
z[2,2]  -1.8   0.0   1.7   0.0 1.0  1.01     1874      768
z[3,2]  -1.6   0.0   1.6   0.0 1.0  1.00     1390      739
alpha    0.0   0.4   1.4   0.5 0.5  1.00     1120      668
lp__   -22.1 -17.4 -14.8 -17.8 2.3  1.00      504      438