Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
estimate Bayes factor
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment