Add Bayes Factor plot via Savage-Dickey density ratio

arviz/plots/bf_plot.py

@@ -0,0 +1,67 @@
+# Plotting and reporting Bayes Factor given idata, var name, prior distribution and reference value
+from scipy import stats
+import matplotlib.pyplot as plt
+import numpy as np
+import logging
+# 
+_log = logging.getLogger(__name__)
+
+def plot_bf(idata, var_name, prior, family = 'normal',  ref_val=0, xlim=None, ax=None):
+    """
+    Bayes Factor approximated as the Savage-Dickey density ratio.
+    The Bayes factor is estimated by comparing a model 
+    against a model in which the parameter of interest has been restricted to a point-null.
+
+    :idata: The "trace" of model, after sampling the
+    :var_name: [str] Name of variable we want to test.
+    :prior: In case we want to use diffent prior (for sensitivity analysis of BF), we can define one and sent it to the function.
+    :family: for now, supports only the normal distribution.
+    :ref_val: reference value for BF testing
+    :xlim: limit the x axis (might be used for visualization porpuses sometimes)
+
+    # grab trace, a variable name to compute difference and prior.
+    # ref_val is the parameter we want to compare
+    # test some elemtns
+    # varName should be string
+    if not isinstance(var_name, str):
+        print('varName is not a string')
+    # BFs based on density estimation (using kernel smoothing instead of spline)
+    """
+    post = extract(idata, var_names=var_name)
+    if prior is None:
+        # grab prior from the data in case it wasn't defined by the user
+        prior = extract(idata, var_names=var_name, group="prior")
+    if post.ndim > 1:
+        print("Posterior distribution has {post.ndim} dimensions")
+    # generate vector
+    if xlim is None:
+        x = np.linspace(np.min(prior), np.max(prior),prior.shape[0])
+    else:
+        x = np.linspace(xlim[0], xlim[1],prior.shape[0])
+    #x = np.linspace(np.min(post), np.max(post),prior.shape[0])
+    my_pdf = stats.gaussian_kde(post)
+    prior_pdf = stats.gaussian_kde(prior)
+    if ax is None:
+        fig, ax = plt.subplots()
+    ax.plot(
+        x, my_pdf(x), "--", lw=2.5, alpha=0.6, label="Posterior"
+    )  # distribution function
+    ax.plot(x, prior_pdf(x), "r-", lw=2.5, alpha=0.6, label="Prior")
+    if ref_val>np.max(post) | ref_val<np.min(post):
+        print('Reference value is out of bounds of posterior')
+    else:
+        posterior = my_pdf(ref_val) # this gives the pdf at ref_val
+    prior = prior_pdf(ref_val)
+    BF10 = posterior / prior
+    BF01 = prior / posterior
+    _log.warning("the Bayes Factor 10 is %.3f" % (BF10)"
+                ("the Bayes Factor 01 is %.3f" % (BF01)
+    )
+    ax.plot(ref_val, posterior, "ko", lw=1.5, alpha=1)
+    ax.plot(ref_val, prior, "ko", lw=1.5, alpha=1)
+    ax.set_xlabel(var_name)
+    ax.set_ylabel("Density")
+    plt.legend()
+
+    return {'BF10': BF10, 'BF01':BF01}, ax
+# end