Diagnostics

Calculate the estimated Bayesian fraction of missing information (BFMI).

    BFMI quantifies how well momentum resampling matches the marginal energy distribution. For more
    information on BFMI, see https://arxiv.org/pdf/1604.00695v1.pdf. The current advice is that
    values smaller than 0.3 indicate poor sampling. However, this threshold is
    provisional and may change. See
    `pystan_workflow <http://mc-stan.org/users/documentation/case-studies/pystan_workflow.html>`_
    for more information.

    Parameters
    ----------
    data : obj
        Any object that can be converted to an :class:`arviz.InferenceData` object.
        Refer to documentation of :func:`arviz.convert_to_dataset` for details.
        If InferenceData, energy variable needs to be found.

    Returns
    -------
    z : array
        The Bayesian fraction of missing information of the model and trace. One element per
        chain in the trace.

    See Also
    --------
    plot_energy : Plot energy transition distribution and marginal energy
                  distribution in HMC algorithms.

    Examples
    --------
    Compute the BFMI of an InferenceData object

    .. ipython::

        In [1]: import arviz as az
           ...: data = az.load_arviz_data('radon')
           ...: az.bfmi(data)

    

ess(data; kwargs...) -> Dataset

Estimate the effective sample size (ESS).

The basic ESS ($N_{\mathit{eff}}$) diagnostic is computed by^{[Vehtari2019][BDA3_ESS]}:

\[\hat{N}_{\mathit{eff}} = \frac{MN}{\hat{\tau}}\\ \hat{\tau} = -1 + 2 \sum_{t'=0}^K \hat{P}_{t'},\]

where $M$ is the number of chains, $N$ the number of draws, $\hat{\rho}_t$ is the estimated autocorrelation at lag $t$, and $K$ is the largest integer for which $\hat{P}_{K} = \hat{\rho}_{2K} + \hat{\rho}_{2K+1}$ is still positive.

The current implementation is similar to Stan's, which uses Geyer's initial monotone sequence criterion^{[StanManual_ESS]}.

See also: rhat, mcse, plot_ess, summarystats

Arguments

• data::Any: Any object that can be converted to an InferenceData object. See convert_to_inference_data for details.

Keywords

• var_names: Symbol names of variables to include in the returned Dataset. Defaults to all.
• method::Symbol=:bulk: ESS method to use. Valid options are:
  • :bulk
  • :tail (specify prob)
  • :quantile (specify prob)
  • :mean: old ESS
  • :sd
  • :median
  • :mad: mean absolute deviance
  • :z_scale
  • :folded
  • :identity
  • :local (specify prob)
• prob::Union{Real,NTuple{2,Real}}: Probability for :tail, :quantile, or :local ESS method.

Examples

Calculate the ESS using the default arguments:

julia> using ArviZ

julia> data = load_example_data("non_centered_eight");

julia> ess(data)
Dataset with dimensions:
  Dim{:school} Categorical{String} String[Choate, Deerfield, …, St. Paul's, Mt. Hermon] Unordered
and 4 layers:
  :mu      Float64 dims:
  :theta_t Float64 dims: Dim{:school} (8)
  :theta   Float64 dims: Dim{:school} (8)
  :tau     Float64 dims:

Calculate the ESS using the :tail method, leaving the prob argument at its default value:

julia> ess(data; method=:tail)
Dataset with dimensions:
  Dim{:school} Categorical{String} String[Choate, Deerfield, …, St. Paul's, Mt. Hermon] Unordered
and 4 layers:
  :mu      Float64 dims:
  :theta_t Float64 dims: Dim{:school} (8)
  :theta   Float64 dims: Dim{:school} (8)
  :tau     Float64 dims:

mcse(data; kwargs...) -> Dataset

Calculate the Markov Chain Standard Error statistic.

See also: ess, plot_mcse, summarystats

Arguments

• data::Any: Any object that can be converted to an InferenceData object. See convert_to_inference_data for details.

Keywords

• var_names: Symbol names of variables to include in the returned Dataset. Defaults to all.

• method::Symbol=:bulk: ESS method to use. Valid options are:

  • :mean:
  • :sd
  • :median
  • :quantile (specify prob)

• prob::Real: Probability for :quantile ESS method.

Examples

Calculate the MCSE using the default arguments:

julia> using ArviZ

julia> data = load_example_data("non_centered_eight");

julia> mcse(data)
Dataset with dimensions:
  Dim{:school} Categorical{String} String[Choate, Deerfield, …, St. Paul's, Mt. Hermon] Unordered
and 4 layers:
  :mu      Float64 dims:
  :theta_t Float64 dims: Dim{:school} (8)
  :theta   Float64 dims: Dim{:school} (8)
  :tau     Float64 dims:

Calculate the MCSE using the :quantile method:

julia> mcse(data; method=:quantile, prob=0.7)
Dataset with dimensions:
  Dim{:school} Categorical{String} String[Choate, Deerfield, …, St. Paul's, Mt. Hermon] Unordered
and 4 layers:
  :mu      Float64 dims:
  :theta_t Float64 dims: Dim{:school} (8)
  :theta   Float64 dims: Dim{:school} (8)
  :tau     Float64 dims:

Compute estimate of rank normalized splitR-hat for a set of traces.

    The rank normalized R-hat diagnostic tests for lack of convergence by comparing the variance
    between multiple chains to the variance within each chain. If convergence has been achieved,
    the between-chain and within-chain variances should be identical. To be most effective in
    detecting evidence for nonconvergence, each chain should have been initialized to starting
    values that are dispersed relative to the target distribution.

    Parameters
    ----------
    data : obj
        Any object that can be converted to an :class:`arviz.InferenceData` object.
        Refer to documentation of :func:`arviz.convert_to_dataset` for details.
        At least 2 posterior chains are needed to compute this diagnostic of one or more
        stochastic parameters.
        For ndarray: shape = (chain, draw).
        For n-dimensional ndarray transform first to dataset with ``az.convert_to_dataset``.
    var_names : list
        Names of variables to include in the rhat report
    method : str
        Select R-hat method. Valid methods are:
        - "rank"        # recommended by Vehtari et al. (2019)
        - "split"
        - "folded"
        - "z_scale"
        - "identity"
    dask_kwargs : dict, optional
        Dask related kwargs passed to :func:`~arviz.wrap_xarray_ufunc`.

    Returns
    -------
    xarray.Dataset
      Returns dataset of the potential scale reduction factors, :math:`\hat{R}`

    See Also
    --------
    ess : Calculate estimate of the effective sample size (ess).
    mcse : Calculate Markov Chain Standard Error statistic.
    plot_forest : Forest plot to compare HDI intervals from a number of distributions.

    Notes
    -----
    The diagnostic is computed by:

      .. math:: \hat{R} = \frac{\hat{V}}{W}

    where :math:`W` is the within-chain variance and :math:`\hat{V}` is the posterior variance
    estimate for the pooled rank-traces. This is the potential scale reduction factor, which
    converges to unity when each of the traces is a sample from the target posterior. Values
    greater than one indicate that one or more chains have not yet converged.

    Rank values are calculated over all the chains with ``scipy.stats.rankdata``.
    Each chain is split in two and normalized with the z-transform following Vehtari et al. (2019).

    References
    ----------
    * Vehtari et al. (2019) see https://arxiv.org/abs/1903.08008
    * Gelman et al. BDA (2014)
    * Brooks and Gelman (1998)
    * Gelman and Rubin (1992)

    Examples
    --------
    Calculate the R-hat using the default arguments:

    .. ipython::

        In [1]: import arviz as az
           ...: data = az.load_arviz_data("non_centered_eight")
           ...: az.rhat(data)

    Calculate the R-hat of some variables using the folded method:

    .. ipython::

        In [1]: az.rhat(data, var_names=["mu", "theta_t"], method="folded")