mcse

cultivars.bayes.chains.mcse(*chains)

Monte Carlo standard error of the posterior mean.

\[\widehat{\operatorname{se}}(\bar x) = \frac{\hat\sigma}{\sqrt{\widehat{S}_{\text{eff}}}},\]

with \(\hat\sigma\) the standard deviation of the pooled draws and \(\widehat{S}_{\text{eff}}\) the Geyer effective sample size of the raw split chains rather than the rank-normalized ones, since it is the untransformed draws whose average is being reported. The number says how far the reported mean is likely to sit from the mean of the distribution the sampler targets; it says nothing about how far that distribution sits from the truth, which is the posterior standard deviation’s job.

Parameters:

*chains (ArrayLike) – One or more (S,) or (S, d) chains of the same shape, with S at least eight and finite throughout.

Return type:

float | NDArray[double]

Returns:

A float for (S,) chains, else (d,), in the units of the quantity. nan for a quantity that never moved.

Raises:

• SpecificationError – If no chain is given, or a chain is shorter than eight draws.

• DimensionError – If a chain is not one- or two-dimensional, or the chains differ in shape.

• NumericalError – If any draw is not finite.

Note

A standard error a tenth of the posterior standard deviation is the usual working target: beyond that, more draws sharpen the reported mean to digits the posterior does not support.

See also

ess_bulk(): The rank-normalized size, which is the one to

read for a convergence verdict.

References

Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., & Bürkner,

P.-C. (2021). Rank-normalization, folding, and localization: An improved \(\widehat{R}\) for assessing convergence of MCMC. Bayesian Analysis, 16(2), 667-718.

Example

Ten thousand independent standard normal draws pin the mean to about \(1 / \sqrt{10000} = 0.01\):

>>> import numpy as np
>>> rng = np.random.default_rng(4)
>>> round(float(mcse(rng.standard_normal(10_000))), 1)
0.0