Merge pull request #38 from aimalz/issue/35/quantilespacing

Issue/35/quantilespacing

docs/notebooks/demo.ipynb

@@ -90,7 +90,8 @@
    },
    "outputs": [],
    "source": [
-    "quantiles = P.quantize(percent=10.)\n",
+    "quants = np.array([0.01,0.02,0.03,0.04,0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,0.95,0.96,0.97,0.98,0.99])\n",
+    "quantiles = P.quantize(quants=quants)#(percent=10.)\n",
     "Q = qp.PDF(quantiles=quantiles)"
    ]
   },

qp/pdf.py

@@ -4,7 +4,6 @@
 
 import qp
 
-<<<<<<< HEAD
 class PDF(object):
 
     def __init__(self, truth=None, quantiles=None, histogram=None,
@@ -37,8 +36,7 @@ def __init__(self, truth=None, quantiles=None, histogram=None,
             self.initialized = 'histogram'
         self.last = self.initialized
 
-        if vb and self.truth is None and self.quantiles is None
-              and self.histogram is None:
+        if vb and self.truth is None and self.quantiles is None and self.histogram is None:
             print 'Warning: initializing a PDF object without inputs'
         self.difs = None
         self.mids = None
@@ -96,12 +94,14 @@ def integrate(self, limits):
         """
         return None
 
-    def quantize(self, percent=1., number=None, vb=True):
+    def quantize(self, quants=None, percent=1., number=None, vb=True):
         """
         Computes an array of evenly-spaced quantiles from the truth.
 
         Parameters
         ----------
+        quants: ndarray, float
+            array of quantile locations as decimals
         percent: float
             the separation of the requested quantiles, in percent
         number: int
@@ -123,25 +123,34 @@ def quantize(self, percent=1., number=None, vb=True):
         object stored in `self.truth`. This calculates the inverse CDF.
         See `the Scipy docs <https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.rv_continuous.ppf.html#scipy.stats.rv_continuous.ppf>`_ for details.
         """
-        if number is not None:
-            # Compute the spacing of the quantiles:
-            quantum = 1.0 / float(number+1)
+        if quants is not None:
+            full_quants = np.append(np.array([0.]),quants)
+            full_quants = np.append(full_quants,np.array([1.]))
+            quantdifs = full_quants[1:]-full_quants[:-1]
+            assert np.sum(quantdifs)==1.
+            self.quantpoints = quants
         else:
-            quantum = percent/100.0
-            # Over-write the number of quantiles:
-            number = np.ceil(100.0 / percent) - 1
-            assert number > 0
-
-        points = np.linspace(0.0+quantum, 1.0-quantum, number)
-        if vb: print("Calculating quantiles: ", points)
+            if number is not None:
+                # Compute the spacing of the quantiles:
+                quantum = 1.0 / float(number+1)
+            else:
+                quantum = percent/100.0
+                # Over-write the number of quantiles:
+                number = np.ceil(100.0 / percent) - 1
+                assert number > 0
+
+            self.quantpoints = np.linspace(0.0+quantum, 1.0-quantum, number)
+        if vb: print("Calculating quantiles: ", self.quantpoints)
         if self.truth is not None:
-            self.quantiles = self.truth.ppf(points)
+            self.quantiles = self.truth.ppf(self.quantpoints)
         else:
             print('New quantiles can only be computed from a truth distribution in this version.')
             return
 
         if vb: print("Result: ", self.quantiles)
+        self.quantvals = self.evaluate(self.quantiles)
         self.last = 'quantiles'
+        self.quantiles = (self.quantiles, self.quantvals)
         return self.quantiles
 
     def histogramize(self, binends=None, nbins=10, binrange=[0., 1.], vb=True):
@@ -222,9 +231,10 @@ def interpolate(self, using=None, vb=True):
             if self.quantiles is None:
                 self.quantiles = self.quantize()
 
-            self.difs = self.quantiles[1:]-self.quantiles[:-1]
-            self.mids = (self.quantiles[1:]+self.quantiles[:-1])/2.
-            self.vals = (1.0/(len(self.quantiles)+1))/self.difs
+            #self.difs = self.quantiles[0][1:]-self.quantiles[0][:-1]
+            #self.mids = (self.quantiles[0][1:]+self.quantiles[0][:-1])/2.
+            self.mids = self.quantiles[0]
+            self.vals = self.quantiles[1]#(1.0/(len(self.quantiles)+1))/self.difs
 
         if using == 'histogram':
             # First find the histogram if none exists:
@@ -300,7 +310,7 @@ def plot(self, limits, points=None):
             plt.plot(x, self.truth.pdf(x), color='k', linestyle='-', lw=1.0, alpha=1.0, label='True PDF')
 
         if self.quantiles is not None:
-            plt.vlines(self.quantiles, 0., 1., color='b', linestyle='--', lw=1.0, alpha=1., label='Quantiles')
+            plt.vlines(self.quantiles[0], np.zeros(len(self.quantiles[1])), self.quantiles[1], color='b', linestyle='--', lw=1.0, alpha=1., label='Quantiles')
             if points is not None:
                 (grid, qinterpolated) = self.approximate(points, using='quantiles')
                 plt.plot(grid, qinterpolated, color='b', linestyle=':', lw=2.0, alpha=1.0, label='Quantile Interpolated PDF')