1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25 import numpy as np
26 from scipy.stats import gaussian_kde
27
29 """Density estimation using a KDE on bounded domains.
30
31 Bounds can be any combination of low or high (if no bound, set to
32 ``float('inf')`` or ``float('-inf')``), and can be periodic or
33 non-periodic. Cannot handle topologies that have
34 multi-dimensional periodicities; will only handle topologies that
35 are direct products of (arbitrary numbers of) R, [0,1], and S1.
36
37 :param pts:
38 ``(Ndim, Npts)`` shaped array of points (as in :class:`gaussian_kde`).
39
40 :param low:
41 Lower bounds; if ``None``, assume no lower bounds.
42
43 :param high:
44 Upper bounds; if ``None``, assume no upper bounds.
45
46 :param periodic:
47 Boolean array giving periodicity in each dimension; if
48 ``None`` assume no dimension is periodic.
49
50 :param bw_method: (optional)
51 Bandwidth estimation method (see :class:`gaussian_kde`)."""
52 - def __init__(self, pts, low=None, high=None, periodic=None, bw_method=None):
53
54 if bw_method is None:
55 bw_method='scott'
56
57 if low is None:
58 self._low = np.zeros(pts.shape[0])
59 self._low[:] = float('-inf')
60 else:
61 self._low = low
62
63 if high is None:
64 self._high = np.zeros(pts.shape[0])
65 self._high[:] = float('inf')
66 else:
67 self._high = high
68
69 if periodic is None:
70 self._periodic = np.zeros(pts.shape[0], dtype=np.bool)
71 else:
72 self._periodic = periodic
73
74 super(BoundedKDE, self).__init__(pts, bw_method=bw_method)
75
77 """Evaluate the KDE at the given points."""
78
79 pts_orig=pts
80 pts=np.copy(pts_orig)
81
82 den=super(BoundedKDE, self).evaluate(pts)
83
84 for i, (low,high,period) in enumerate(zip(self._low, self._high, self._periodic)):
85 if period:
86 P = high-low
87
88 pts[i, :] += P
89 den += super(BoundedKDE, self).evaluate(pts)
90
91 pts[i,:] -= 2.0*P
92 den += super(BoundedKDE, self).evaluate(pts)
93
94 pts[i,:]=pts_orig[i,:]
95
96 else:
97 if not (low == float('-inf')):
98 pts[i,:] = 2.0*low - pts[i,:]
99 den += super(BoundedKDE, self).evaluate(pts)
100 pts[i,:] = pts_orig[i,:]
101
102 if not (high == float('inf')):
103 pts[i,:] = 2.0*high - pts[i,:]
104 den += super(BoundedKDE, self).evaluate(pts)
105 pts[i,:] = pts_orig[i,:]
106
107 return den
108
109 __call__ = evaluate
110
112 """Quantile of ``pt``, evaluated by a greedy algorithm.
113
114 :param pt:
115 The point at which the quantile value is to be computed.
116
117 The quantile of ``pt`` is the fraction of points used to
118 construct the KDE that have a lower KDE density than ``pt``."""
119
120 return float(np.count_nonzero(self.evaluate(self.dataset) < self.evaluate(pt)))/float(self.n)
121