estimate Bayes factor

bayes_factor.py

import numpy as np
from astroML import stats

def estimate_bayes_factor(traces, logp, r=0.05, return_list=False,
                          old_version=False, normalize_space=True):
    """Estimate the bayes factor using the local density of points"""
    D, N = traces.shape
    
    if normalize_space:
        traces = traces.copy()
        for i in range(traces.shape[0]):
            #traces[i] /= traces[i].std()
            traces[i] /= sigmaG(traces[i])

    if old_version:
        # use neighbor count within r as a density estimator
        bt = BallTree(traces.T)
        count = bt.query_radius(traces.T, r=r, count_only=True)

        # compute volume of a D-dimensional sphere of radius r
        Vr = np.pi ** (0.5 * D) / gamma(0.5 * D + 1) * (r ** D)

        BF = logp + np.log(N) + np.log(Vr) - np.log(count)
    else:
        bt = BallTree(traces.T)
        log_density = bt.kernel_density(traces.T, r, return_log=True)

        BF = logp + np.log(N) - log_density

    if return_list:
        return BF
    else:
        p25, p50, p75 = np.percentile(BF, [25, 50, 75])
        return p50, 0.7413 * (p75 - p25)