def estimate_bayes_factor(traces, logp, r=0.05):
"""
Esitmate Odds ratios in a random subsample of the chains in MCMC
AstroML (see Eqn 5.127, pg 237)
"""
ndim, nsteps = traces.shape # [ndim,number of steps in chain]
# compute volume of a n-dimensional (ndim) sphere of radius r
Vr = np.pi ** (0.5 * ndim) / gamma(0.5 * ndim + 1) * (r ** ndim)
# use neighbor count within r as a density estimator
bt = BallTree(traces.T)
count = bt.query_radius(traces.T, r=r, count_only=True)
# BF = N*p/rho
bf = logp + np.log(nsteps) + np.log(Vr) - np.log(count) #log10(bf)
p25, p50, p75 = np.percentile(bf, [25, 50, 75])
return p50, 0.7413 * (p75 - p25)
########################################################################
python类BallTree()的实例源码
def buildNNDataStructure(self):
"""Builds a nearest neighbor data structure. User doesn't need to
call this unless the self.problems attribute was changed manually."""
if len(self.problemFeatures)==0 or len(self.featureNames)==0:
return
try:
from sklearn.neighbors import NearestNeighbors,BallTree
from scipy.spatial import KDTree
with self.lock:
try:
farray = self.problemFeatures.array
except AttributeError:
farray = np.array(self.problemFeatures.items)
if self.metricTransform is not None:
farray = np.dot(farray,self.metricTransform)
#self.nn = NearestNeighbors(n_neighbors=1,algorithm="auto").fit(farray)
self.nn = BallTree(farray)
#self.nn = KDTree(farray)
self.nnBuildSize = len(self.problemFeatures)
except ImportError:
print "IKDatabase: Warning, scikit-learn is not installed, queries will be much slower"
with self.lock:
self.nn = None
self.nnBuildSize = 0
return
def _point_cloud_error_balltree(src_pts, tgt_tree):
"""Find the distance from each source point to its closest target point
Uses sklearn.neighbors.BallTree for greater efficiency
Parameters
----------
src_pts : array, shape = (n, 3)
Source points.
tgt_tree : sklearn.neighbors.BallTree
BallTree of the target points.
Returns
-------
dist : array, shape = (n, )
For each point in ``src_pts``, the distance to the closest point in
``tgt_pts``.
"""
dist, _ = tgt_tree.query(src_pts)
return dist.ravel()
def test_unsupervised_inputs():
# test the types of valid input into NearestNeighbors
X = rng.random_sample((10, 3))
nbrs_fid = neighbors.NearestNeighbors(n_neighbors=1)
nbrs_fid.fit(X)
dist1, ind1 = nbrs_fid.kneighbors(X)
nbrs = neighbors.NearestNeighbors(n_neighbors=1)
for input in (nbrs_fid, neighbors.BallTree(X), neighbors.KDTree(X)):
nbrs.fit(input)
dist2, ind2 = nbrs.kneighbors(X)
assert_array_almost_equal(dist1, dist2)
assert_array_almost_equal(ind1, ind2)
def build_neighbor_index(x, leaf_size):
return sk_neighbors.BallTree(x, leaf_size=leaf_size)
def indexBallTree(X,leafSize):
tree = BallTree(X, leaf_size=leafSize)
return tree
#typeDescriptors =SURF, SIFT, OEB
#k = number of images to be retrieved
def update_tree(self, time):
print 'rebuild tree'
self.tree = BallTree(self.state[:self.items, :], leaf_size=self.size)
self.last_tree_built_time = time
print 'rebuild done'
def __init__(self, points, rho, dimension):
"""Constructor
Initializes the grid and helper structures using the provided points
and rho parameter.
Args:
points: A numpy array containing the coordinates of the particles.
rho: Needed to compute the rho-boundary of the system.
dimension: The dimension of the particle system.
"""
self.points = points
self.rho = rho
self.dimension = dimension
self.cell_size = 2.0 * rho
self.aabb_min = np.amin(points, axis=0)
self.aabb_max = np.amax(points, axis=0)
self.grid_dims = (self.aabb_max - self.aabb_min) / self.cell_size
# Regarding the + 3: 1 for left side, 1 for right side, 1 for rounding
# up
self.grid_dims = np.trunc(self.grid_dims) + 3
self.grid_dims = self.grid_dims.astype(int)
self.grid_min = self.aabb_min - self.cell_size
self.grid_max = self.grid_min + self.grid_dims * self.cell_size
self.grid_count = np.zeros(self.grid_dims, dtype=int)
self.grid_elems = np.empty(self.grid_dims, dtype=object)
self.update_grid()
self.tree = NeighborsTree(
self.points, leaf_size=10, metric='euclidean')
self.neighbor_cell_list = self.compute_neighbor_cell_list()
def test_barnes_hut_angle():
# When Barnes-Hut's angle=0 this corresponds to the exact method.
angle = 0.0
perplexity = 10
n_samples = 100
for n_components in [2, 3]:
n_features = 5
degrees_of_freedom = float(n_components - 1.0)
random_state = check_random_state(0)
distances = random_state.randn(n_samples, n_features)
distances = distances.astype(np.float32)
distances = distances.dot(distances.T)
np.fill_diagonal(distances, 0.0)
params = random_state.randn(n_samples, n_components)
P = _joint_probabilities(distances, perplexity, False)
kl, gradex = _kl_divergence(params, P, degrees_of_freedom, n_samples,
n_components)
k = n_samples - 1
bt = BallTree(distances)
distances_nn, neighbors_nn = bt.query(distances, k=k + 1)
neighbors_nn = neighbors_nn[:, 1:]
Pbh = _joint_probabilities_nn(distances, neighbors_nn,
perplexity, False)
kl, gradbh = _kl_divergence_bh(params, Pbh, neighbors_nn,
degrees_of_freedom, n_samples,
n_components, angle=angle,
skip_num_points=0, verbose=False)
assert_array_almost_equal(Pbh, P, decimal=5)
assert_array_almost_equal(gradex, gradbh, decimal=5)
