Source Code for Module pylal.kdes

  1  # -*- coding: utf-8 -*- 
  2  # 
  3  #       kdes.py: KDE utilities for density estimation in unusual 
  4  #       topologies. 
  5  # 
  6  #       Copyright 2012 
  7  #       Will M. Farr <will.farr@ligo.org> 
  8  # 
  9  #       This program is free software; you can redistribute it and/or modify 
 10  #       it under the terms of the GNU General Public License as published by 
 11  #       the Free Software Foundation; either version 2 of the License, or 
 12  #       (at your option) any later version. 
 13  # 
 14  #       This program is distributed in the hope that it will be useful, 
 15  #       but WITHOUT ANY WARRANTY; without even the implied warranty of 
 16  #       MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
 17  #       GNU General Public License for more details. 
 18  # 
 19  #       You should have received a copy of the GNU General Public License 
 20  #       along with this program; if not, write to the Free Software 
 21  #       Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, 
 22  #       MA 02110-1301, USA. 
 23   
 24   
 25  import numpy as np 
 26  from scipy.stats import gaussian_kde 
 27   
 28 -class BoundedKDE(gaussian_kde): 
 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   
 76 -    def evaluate(self, pts): 
 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   
111 -    def quantile(self, pt): 
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