Create a data frame with summary statistics.

    Parameters
    ----------
    data : obj
        Any object that can be converted to an az.InferenceData object
        Refer to documentation of az.convert_to_dataset for details
    var_names : list
        Names of variables to include in summary
    include_circ : bool
        Whether to include circular statistics
    fmt : {'wide', 'long', 'xarray'}
        Return format is either pandas.DataFrame {'wide', 'long'} or xarray.Dataset {'xarray'}.
    round_to : int
        Number of decimals used to round results. Defaults to 2.
    stat_funcs : dict
        A list of functions or a dict of functions with function names as keys used to calculate
        statistics. By default, the mean, standard deviation, simulation standard error, and
        highest posterior density intervals are included.

        The functions will be given one argument, the samples for a variable as an nD array,
        The functions should be in the style of a ufunc and return a single number. For example,
        `np.mean`, or `scipy.stats.var` would both work.
    extend : boolean
        If True, use the statistics returned by `stat_funcs` in addition to, rather than in place
        of, the default statistics. This is only meaningful when `stat_funcs` is not None.
    credible_interval : float, optional
        Credible interval to plot. Defaults to 0.94. This is only meaningful when `stat_funcs` is
        None.
    order : {"C", "F"}
        If fmt is "wide", use either C or F unpacking order. Defaults to C.

    Returns
    -------
    pandas.DataFrame
        With summary statistics for each variable. Defaults statistics are: `mean`, `sd`,
        `hpd_3%`, `hpd_97%`, `mc_error`, `ess` and `r_hat`. `ess` and `r_hat` are only computed
        for traces with 2 or more chains.

    Examples
    --------
    .. code:: ipython

        >>> az.summary(trace, ['mu'])
               mean    sd  mc_error  hpd_3  hpd_97  ess  r_hat
        mu[0]  0.10  0.06      0.00  -0.02    0.23  487.0  1.00
        mu[1] -0.04  0.06      0.00  -0.17    0.08  379.0  1.00

    Other statistics can be calculated by passing a list of functions.

    .. code:: ipython

        >>> import pandas as pd
        >>> def trace_sd(x):
        ...     return pd.Series(np.std(x, 0), name='sd')
        ...
        >>> def trace_quantiles(x):
        ...     return pd.DataFrame(pd.quantiles(x, [5, 50, 95]))
        ...
        >>> az.summary(trace, ['mu'], stat_funcs=[trace_sd, trace_quantiles])
                 sd     5    50    95
        mu[0]  0.06  0.00  0.10  0.21
        mu[1]  0.07 -0.16 -0.04  0.06

def summary(
    data,
    var_names=None,
    fmt="wide",
    round_to=2,
    include_circ=None,
    stat_funcs=None,
    extend=True,
    credible_interval=0.94,
    order="C",
):
    """Create a data frame with summary statistics.

    Parameters
    ----------
    data : obj
        Any object that can be converted to an az.InferenceData object
        Refer to documentation of az.convert_to_dataset for details
    var_names : list
        Names of variables to include in summary
    include_circ : bool
        Whether to include circular statistics
    fmt : {'wide', 'long', 'xarray'}
        Return format is either pandas.DataFrame {'wide', 'long'} or xarray.Dataset {'xarray'}.
    round_to : int
        Number of decimals used to round results. Defaults to 2.
    stat_funcs : dict
        A list of functions or a dict of functions with function names as keys used to calculate
        statistics. By default, the mean, standard deviation, simulation standard error, and
        highest posterior density intervals are included.

        The functions will be given one argument, the samples for a variable as an nD array,
        The functions should be in the style of a ufunc and return a single number. For example,
        `np.mean`, or `scipy.stats.var` would both work.
    extend : boolean
        If True, use the statistics returned by `stat_funcs` in addition to, rather than in place
        of, the default statistics. This is only meaningful when `stat_funcs` is not None.
    credible_interval : float, optional
        Credible interval to plot. Defaults to 0.94. This is only meaningful when `stat_funcs` is
        None.
    order : {"C", "F"}
        If fmt is "wide", use either C or F unpacking order. Defaults to C.

    Returns
    -------
    pandas.DataFrame
        With summary statistics for each variable. Defaults statistics are: `mean`, `sd`,
        `hpd_3%`, `hpd_97%`, `mc_error`, `ess` and `r_hat`. `ess` and `r_hat` are only computed
        for traces with 2 or more chains.

    Examples
    --------
    .. code:: ipython

        >>> az.summary(trace, ['mu'])
               mean    sd  mc_error  hpd_3  hpd_97  ess  r_hat
        mu[0]  0.10  0.06      0.00  -0.02    0.23  487.0  1.00
        mu[1] -0.04  0.06      0.00  -0.17    0.08  379.0  1.00

    Other statistics can be calculated by passing a list of functions.

    .. code:: ipython

        >>> import pandas as pd
        >>> def trace_sd(x):
        ...     return pd.Series(np.std(x, 0), name='sd')
        ...
        >>> def trace_quantiles(x):
        ...     return pd.DataFrame(pd.quantiles(x, [5, 50, 95]))
        ...
        >>> az.summary(trace, ['mu'], stat_funcs=[trace_sd, trace_quantiles])
                 sd     5    50    95
        mu[0]  0.06  0.00  0.10  0.21
        mu[1]  0.07 -0.16 -0.04  0.06
    """
    posterior = convert_to_dataset(data, group="posterior")
    var_names = _var_names(var_names, posterior)
    posterior = posterior if var_names is None else posterior[var_names]

    fmt_group = ("wide", "long", "xarray")
    if not isinstance(fmt, str) or (fmt.lower() not in fmt_group):
        raise TypeError("Invalid format: '{}'! Formatting options are: {}".format(fmt, fmt_group))

    unpack_order_group = ("C", "F")
    if not isinstance(order, str) or (order.upper() not in unpack_order_group):
        raise TypeError(
            "Invalid order: '{}'! Unpacking options are: {}".format(order, unpack_order_group)
        )

    alpha = 1 - credible_interval

    metrics = []
    metric_names = []

    if stat_funcs is not None:
        if isinstance(stat_funcs, dict):
            for stat_func_name, stat_func in stat_funcs.items():
                metrics.append(
                    xr.apply_ufunc(
                        _make_ufunc(stat_func), posterior, input_core_dims=(("chain", "draw"),)
                    )
                )
                metric_names.append(stat_func_name)
        else:
            for stat_func in stat_funcs:
                metrics.append(
                    xr.apply_ufunc(
                        _make_ufunc(stat_func), posterior, input_core_dims=(("chain", "draw"),)
                    )
                )
                metric_names.append(stat_func.__name__)

    if extend:
        metrics.append(posterior.mean(dim=("chain", "draw")))
        metric_names.append("mean")

        metrics.append(posterior.std(dim=("chain", "draw")))
        metric_names.append("sd")

        metrics.append(
            xr.apply_ufunc(_make_ufunc(_mc_error), posterior, input_core_dims=(("chain", "draw"),))
        )
        metric_names.append("mc error")

        metrics.append(
            xr.apply_ufunc(
                _make_ufunc(hpd, index=0, credible_interval=credible_interval),
                posterior,
                input_core_dims=(("chain", "draw"),),
            )
        )
        metric_names.append("hpd {:g}%".format(100 * alpha / 2))

        metrics.append(
            xr.apply_ufunc(
                _make_ufunc(hpd, index=1, credible_interval=credible_interval),
                posterior,
                input_core_dims=(("chain", "draw"),),
            )
        )
        metric_names.append("hpd {:g}%".format(100 * (1 - alpha / 2)))

    if include_circ:
        metrics.append(
            xr.apply_ufunc(
                _make_ufunc(st.circmean, high=np.pi, low=-np.pi),
                posterior,
                input_core_dims=(("chain", "draw"),),
            )
        )
        metric_names.append("circular mean")

        metrics.append(
            xr.apply_ufunc(
                _make_ufunc(st.circstd, high=np.pi, low=-np.pi),
                posterior,
                input_core_dims=(("chain", "draw"),),
            )
        )
        metric_names.append("circular standard deviation")

        metrics.append(
            xr.apply_ufunc(
                _make_ufunc(_mc_error, circular=True),
                posterior,
                input_core_dims=(("chain", "draw"),),
            )
        )
        metric_names.append("circular mc error")

        metrics.append(
            xr.apply_ufunc(
                _make_ufunc(hpd, index=0, credible_interval=credible_interval, circular=True),
                posterior,
                input_core_dims=(("chain", "draw"),),
            )
        )
        metric_names.append("circular hpd {:.2%}".format(alpha / 2))

        metrics.append(
            xr.apply_ufunc(
                _make_ufunc(hpd, index=1, credible_interval=credible_interval, circular=True),
                posterior,
                input_core_dims=(("chain", "draw"),),
            )
        )
        metric_names.append("circular hpd {:.2%}".format(1 - alpha / 2))

    if len(posterior.chain) > 1:
        metrics.append(effective_sample_size(posterior, var_names=var_names))
        metric_names.append("ess")

        metrics.append(rhat(posterior, var_names=var_names))
        metric_names.append("r_hat")

    joined = xr.concat(metrics, dim="metric").assign_coords(metric=metric_names)

    if fmt.lower() == "wide":
        dfs = []
        for var_name, values in joined.data_vars.items():
            if len(values.shape[1:]):
                metric = list(values.metric.values)
                data_dict = OrderedDict()
                for idx in np.ndindex(values.shape[1:] if order == "C" else values.shape[1:][::-1]):
                    if order == "F":
                        idx = tuple(idx[::-1])
                    ser = pd.Series(values[(Ellipsis, *idx)].values, index=metric)
                    key = "{}[{}]".format(var_name, ",".join(map(str, idx)))
                    data_dict[key] = ser
                df = pd.DataFrame.from_dict(data_dict, orient="index")
                df = df.loc[list(data_dict.keys())]
            else:
                df = values.to_dataframe()
                df.index = list(df.index)
                df = df.T
            dfs.append(df)
        summary_df = pd.concat(dfs, sort=False)
    elif fmt.lower() == "long":
        df = joined.to_dataframe().reset_index().set_index("metric")
        df.index = list(df.index)
        summary_df = df
    else:
        summary_df = joined
    return summary_df.round(round_to)