Improved MCSE

[sethaxen]

See corresponding issues in posterior (stan-dev/posterior#232) and ArviZ (arviz-devs/arviz#2005)

[sethaxen]

More broadly, we should update our implementations of mcse in several ways:

• accept 3D arrays as inputs, where the estimator takes into account between-chain convergence. (e.g. the imse method for a single chain should yield a similar estimate to using ess_rhat for a single chain without splitting. But a better estimate with poor mixing comes from using the ESS estimate)
• accept an optional estimator argument, like ess_rhat in Improvements to ess_rhat #53. Estimators have been worked out for several functions (mean, std, quantile, see e.g. posterior: https://github.com/stan-dev/posterior/blob/879bb0f72f67dfb67637bac2d05abcc1fa173a3b/R/convergence.R)
• update batch-means approach to work for any estimator (Edit: one must work out a specific batch-means estimator for each individual function)
• add subsampling bootstrap method (SBM; see [WIP] draft subsampling bootstrap for mcse arviz-devs/arviz#1974)

I'm not 100% certain what the API should look like. My inclination is that we should ultimately move away from our existing implementations, since they ignore between-chain variance and so likely underestimate the MCSE. The ESS approach should be used by default for the estimators that use it (accepting a method keyword forwarded to ess_rhat), and this replaces imse. I'm not certain how useful ipse is by comparison. For the fallback for generic estimators, we should benchmark bootstrapping vs batch-means use SBM.