Repository: arviz-devs/arviz

Issue:
Refactor Bayes Factor computation out of az.plot_bf
The function `az.plot_bf` is very useful for visualizing Bayes Factors. However, some users may want to use Bayes Factors outside the context of the plot, so it would be useful to refactor the logic to compute the BF from `arviz.plots` to `arviz.stats` and make it part of the public API. 

This seems like a very small change, and I'm happy to do it myself. But I am new to open source and have not communicated with you all before. So before I make the PR I wanted to make sure there's not a good reason to keep it the way it is right now.

Discussion:
@aloctavodia: Hi @mpmbq2, thanks for offering help. Your suggestion is a good one and you can send the PR if you want.

I just want to let you know that we are refactoring ArviZ. As part of the refactor, we are splitting functionality into 3 sub-packages and adding most of the new features in those packages (although there is no strict rule against adding new features to the main repo). 

In the context of your proposal, the computation of the BF  should live in [arviz-stats](https://github.com/arviz-devs/arviz-stats) and the plotting should be on [arvivz-plots](https://github.com/arviz-devs/arviz-plots). Making a PR into one or both of those packages will be more work than making the changes here.

Both options will be helpful for us and the community, so please select the one that you like the most.@mpmbq2: Thanks for the heads up RE the package migration. For the sake of not over-thinking it, I decided to go ahead and make a PR to this repo. Once I get it accepted, I'll have a look at submitting PRs to the stats and plots repos.@PiyushPanwarFST: Sir, I would like to work on this issue. Could you please let me know if anyone else is already working on it?@aloctavodia: Hi @PiyushPanwarFST thanks for your help. There is a PR already open addressing this issue for this repository. Maybe you work on addressing this issue on [arviz-stats](https://github.com/arviz-devs/arviz-stats) instead.@PiyushPanwarFST: Sir can you please give me the clarity where I should contribute regarding the Bayes Factor functionality. From what I understand:

The work on the main repository (arviz) is already completed.
I should focus on arviz-stats for this functionality.
I also noticed there isn’t an open issue related to Bayes Factor computation in arviz-stats. Should I create an issue for it, or can I directly raise a PR to add the functionality? Also, should I clone the arviz-stats repository and start contributing there directly?

Looking forward to your guidance!@aloctavodia: Yes, this is being addressed in the main repo, then better to focus on arviz-stats. You can open an issue if you want. Yes clone the repo and contribute there.

Solution patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index e5965d5cca..4798a43ec7 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,7 +7,8 @@
 ### Maintenance and fixes
 - Make `arviz.data.generate_dims_coords` handle `dims` and `default_dims` consistently ([2395](https://github.com/arviz-devs/arviz/pull/2395))
 - Only emit a warning for custom groups in `InferenceData` when explicitly requested ([2401](https://github.com/arviz-devs/arviz/pull/2401))
-- Update `method="sd"` of `mcse` to not use normality assumption ([2167](https://github.com/arviz-devs/arviz/pull/2167))
+- Splits Bayes Factor computation out from `az.plot_bf` into `az.bayes_factor` ([2402](https://github.com/arviz-devs/arviz/issues/2402))
+- Update `method="sd"` of `mcse` to not use normality assumption ([2167](https://github.com/arviz-devs/arviz/pull/2167)) 
 
 ### Documentation
 - Add example of ECDF comparison plot to gallery ([2178](https://github.com/arviz-devs/arviz/pull/2178))
diff --git a/arviz/plots/bfplot.py b/arviz/plots/bfplot.py
index 8a47178ba3..b95fa3f9c7 100644
--- a/arviz/plots/bfplot.py
+++ b/arviz/plots/bfplot.py
@@ -2,11 +2,9 @@
 # pylint: disable=unbalanced-tuple-unpacking
 import logging
 
-from numpy import interp
-
 from ..data.utils import extract
 from .plot_utils import get_plotting_function
-from ..stats.density_utils import _kde_linear
+from ..stats import bayes_factor
 
 _log = logging.getLogger(__name__)
 
@@ -94,32 +92,17 @@ def plot_bf(
         >>> az.plot_bf(idata, var_name="a", ref_val=0)
 
     """
-    posterior = extract(idata, var_names=var_name).values
-
-    if ref_val > posterior.max() or ref_val < posterior.min():
-        _log.warning(
-            "The reference value is outside of the posterior. "
-            "This translate into infinite support for H1, which is most likely an overstatement."
-        )
-
-    if posterior.ndim > 1:
-        _log.warning("Posterior distribution has {posterior.ndim} dimensions")
 
     if prior is None:
         prior = extract(idata, var_names=var_name, group="prior").values
 
-    if posterior.dtype.kind == "f":
-        posterior_grid, posterior_pdf = _kde_linear(posterior)
-        prior_grid, prior_pdf = _kde_linear(prior)
-        posterior_at_ref_val = interp(ref_val, posterior_grid, posterior_pdf)
-        prior_at_ref_val = interp(ref_val, prior_grid, prior_pdf)
-
-    elif posterior.dtype.kind == "i":
-        posterior_at_ref_val = (posterior == ref_val).mean()
-        prior_at_ref_val = (prior == ref_val).mean()
+    bf, p_at_ref_val = bayes_factor(
+        idata, var_name, prior=prior, ref_val=ref_val, return_ref_vals=True
+    )
+    bf_10 = bf["BF10"]
+    bf_01 = bf["BF01"]
 
-    bf_10 = prior_at_ref_val / posterior_at_ref_val
-    bf_01 = 1 / bf_10
+    posterior = extract(idata, var_names=var_name)
 
     bfplot_kwargs = dict(
         ax=ax,
@@ -128,8 +111,8 @@ def plot_bf(
         prior=prior,
         posterior=posterior,
         ref_val=ref_val,
-        prior_at_ref_val=prior_at_ref_val,
-        posterior_at_ref_val=posterior_at_ref_val,
+        prior_at_ref_val=p_at_ref_val["prior"],
+        posterior_at_ref_val=p_at_ref_val["posterior"],
         var_name=var_name,
         colors=colors,
         figsize=figsize,
diff --git a/arviz/stats/__init__.py b/arviz/stats/__init__.py
index 25bda494f9..e3fb0c14c8 100644
--- a/arviz/stats/__init__.py
+++ b/arviz/stats/__init__.py
@@ -9,6 +9,7 @@
 
 __all__ = [
     "apply_test_function",
+    "bayes_factor",
     "bfmi",
     "compare",
     "hdi",
diff --git a/arviz/stats/stats.py b/arviz/stats/stats.py
index e4c78d3bcb..4f11230a4c 100644
--- a/arviz/stats/stats.py
+++ b/arviz/stats/stats.py
@@ -27,6 +27,7 @@
 from .density_utils import get_bins as _get_bins
 from .density_utils import histogram as _histogram
 from .density_utils import kde as _kde
+from .density_utils import _kde_linear
 from .diagnostics import _mc_error, _multichain_statistics, ess
 from .stats_utils import ELPDData, _circular_standard_deviation, smooth_data
 from .stats_utils import get_log_likelihood as _get_log_likelihood
@@ -41,6 +42,7 @@
 
 __all__ = [
     "apply_test_function",
+    "bayes_factor",
     "compare",
     "hdi",
     "loo",
@@ -2337,3 +2339,87 @@ def _cjs_dist(draws, weights):
     bound = cdf_p_int + cdf_q_int
 
     return np.sqrt((cjs_pq + cjs_qp) / bound)
+
+
+def bayes_factor(idata, var_name, ref_val=0, prior=None, return_ref_vals=False):
+    r"""Approximated Bayes Factor for comparing hypothesis of two nested models.
+
+    The Bayes factor is estimated by comparing a model (H1) against a model in which the
+    parameter of interest has been restricted to be a point-null (H0). This computation
+    assumes the models are nested and thus H0 is a special case of H1.
+
+    Notes
+    -----
+    The bayes Factor is approximated as the Savage-Dickey density ratio
+    algorithm presented in [1]_.
+
+    Parameters
+    ----------
+    idata : InferenceData
+        Any object that can be converted to an :class:`arviz.InferenceData` object
+        Refer to documentation of :func:`arviz.convert_to_dataset` for details.
+    var_name : str, optional
+        Name of variable we want to test.
+    ref_val : int, default 0
+        Point-null for Bayes factor estimation.
+    prior : numpy.array, optional
+        In case we want to use different prior, for example for sensitivity analysis.
+    return_ref_vals : bool, optional
+        Whether to return the values of the prior and posterior at the reference value.
+        Used by :func:`arviz.plot_bf` to display the distribution comparison.
+
+
+    Returns
+    -------
+    dict : A dictionary with BF10 (Bayes Factor 10 (H1/H0 ratio), and BF01 (H0/H1 ratio).
+
+    References
+    ----------
+    .. [1] Heck, D., 2019. A caveat on the Savage-Dickey density ratio:
+       The case of computing Bayes factors for regression parameters.
+
+    Examples
+    --------
+    Moderate evidence indicating that the parameter "a" is different from zero.
+
+    .. ipython::
+
+        In [1]: import numpy as np
+           ...: import arviz as az
+           ...: idata = az.from_dict(posterior={"a":np.random.normal(1, 0.5, 5000)},
+           ...: prior={"a":np.random.normal(0, 1, 5000)})
+           ...: az.bayes_factor(idata, var_name="a", ref_val=0)
+
+    """
+
+    posterior = extract(idata, var_names=var_name).values
+
+    if ref_val > posterior.max() or ref_val < posterior.min():
+        _log.warning(
+            "The reference value is outside of the posterior. "
+            "This translate into infinite support for H1, which is most likely an overstatement."
+        )
+
+    if posterior.ndim > 1:
+        _log.warning("Posterior distribution has {posterior.ndim} dimensions")
+
+    if prior is None:
+        prior = extract(idata, var_names=var_name, group="prior").values
+
+    if posterior.dtype.kind == "f":
+        posterior_grid, posterior_pdf, *_ = _kde_linear(posterior)
+        prior_grid, prior_pdf, *_ = _kde_linear(prior)
+        posterior_at_ref_val = np.interp(ref_val, posterior_grid, posterior_pdf)
+        prior_at_ref_val = np.interp(ref_val, prior_grid, prior_pdf)
+
+    elif posterior.dtype.kind == "i":
+        posterior_at_ref_val = (posterior == ref_val).mean()
+        prior_at_ref_val = (prior == ref_val).mean()
+
+    bf_10 = prior_at_ref_val / posterior_at_ref_val
+    bf = {"BF10": bf_10, "BF01": 1 / bf_10}
+
+    if return_ref_vals:
+        return (bf, {"prior": prior_at_ref_val, "posterior": posterior_at_ref_val})
+    else:
+        return bf


Test patch:
diff --git a/arviz/tests/base_tests/test_plots_matplotlib.py b/arviz/tests/base_tests/test_plots_matplotlib.py
index 5afcbe1669..c2646a73f1 100644
--- a/arviz/tests/base_tests/test_plots_matplotlib.py
+++ b/arviz/tests/base_tests/test_plots_matplotlib.py
@@ -2082,7 +2082,5 @@ def test_plot_bf():
     idata = from_dict(
         posterior={"a": np.random.normal(1, 0.5, 5000)}, prior={"a": np.random.normal(0, 1, 5000)}
     )
-    bf_dict0, _ = plot_bf(idata, var_name="a", ref_val=0)
-    bf_dict1, _ = plot_bf(idata, prior=np.random.normal(0, 10, 5000), var_name="a", ref_val=0)
-    assert bf_dict0["BF10"] > bf_dict0["BF01"]
-    assert bf_dict1["BF10"] < bf_dict1["BF01"]
+    _, bf_plot = plot_bf(idata, var_name="a", ref_val=0)
+    assert bf_plot is not None
diff --git a/arviz/tests/base_tests/test_stats.py b/arviz/tests/base_tests/test_stats.py
index 57b4bf16b5..12b19c2e96 100644
--- a/arviz/tests/base_tests/test_stats.py
+++ b/arviz/tests/base_tests/test_stats.py
@@ -18,6 +18,7 @@
 from ...rcparams import rcParams
 from ...stats import (
     apply_test_function,
+    bayes_factor,
     compare,
     ess,
     hdi,
@@ -871,3 +872,13 @@ def test_priorsens_coords(psens_data):
     assert "mu" in result
     assert "theta" in result
     assert "school" in result.theta_t.dims
+
+
+def test_bayes_factor():
+    idata = from_dict(
+        posterior={"a": np.random.normal(1, 0.5, 5000)}, prior={"a": np.random.normal(0, 1, 5000)}
+    )
+    bf_dict0 = bayes_factor(idata, var_name="a", ref_val=0)
+    bf_dict1 = bayes_factor(idata, prior=np.random.normal(0, 10, 5000), var_name="a", ref_val=0)
+    assert bf_dict0["BF10"] > bf_dict0["BF01"]
+    assert bf_dict1["BF10"] < bf_dict1["BF01"]