Fix problems with CDF method (#15)

CDF evaluation was previously problematic in a number of ways. This should both correct it and make things cleaner.

---------

Co-authored-by: Caden Gobat <[email protected]>

Author: mikemoss3

asymmetric_uncertainty/core.py

@@ -14,6 +14,9 @@
 import numpy as np
 import matplotlib.pyplot as plt
 import warnings
+from math import erf
+
+np_erf = np.vectorize(erf, doc="Vectorized error function at x.")  # prevents need for SciPy dependency
 
 class a_u:
     """
@@ -74,29 +77,33 @@ def _repr_latex_(self):
         else:
             return "$%f_{-%f}^{+%f}$" %(self.value,self.minus,self.plus)
 
-    def pdf(self,x):
+    def pdf(self, x):
         """
         Computes and returns the values of the probability distribution function for the specified input.
         """
-        return np.piecewise(x, [x<self.value, x>=self.value],
-                            [lambda x : np.sqrt(2)/np.sqrt(np.pi)/(self.plus+self.minus) * np.exp(-1*(x-self.value)**2 / (2*self.minus**2)),
-                             lambda x : np.sqrt(2)/np.sqrt(np.pi)/(self.plus+self.minus) * np.exp(-1*(x-self.value)**2 / (2*self.plus**2))])
+        x_arr = np.asanyarray(x).astype(float)
+        pdf_arr =  np.piecewise(x_arr, [x_arr<self.value, x_arr>=self.value],
+                                [lambda x: np.sqrt(2)/np.sqrt(np.pi)/(self.plus+self.minus) * np.exp(-1*(x-self.value)**2 / (2*self.minus**2)),
+                                 lambda x: np.sqrt(2)/np.sqrt(np.pi)/(self.plus+self.minus) * np.exp(-1*(x-self.value)**2 / (2*self.plus**2))])
+        return pdf_arr.item() if np.isscalar(x) else pdf_arr
 
-    def cdf(self,x):
+    def cdf(self, x):
         """
         Computes and returns the values of the cumulative distribution function for the specified input.
         """
-        return np.cumsum(self.pdf(x))/np.sum(self.pdf(x))
-        
+        x_arr = np.asanyarray(x).astype(float)
+        cdf_arr = np.piecewise(x_arr, [x_arr<self.value, x_arr>=self.value],
+                               [lambda x: 0.5*(1 + np_erf((x-self.value)/(self.minus*np.sqrt(2)))),
+                                lambda x: 0.5*(1 + np_erf((x-self.value)/(self.plus*np.sqrt(2))))])
+        return cdf_arr.item() if np.isscalar(x) else cdf_arr
+    
     def pdfplot(self,num_sigma=5,discretization=100,**kwargs):
         """
         Plots the associated PDF over the specified number of sigma, using 2*`discretization` points.
         `**kwargs` are passed on to `matplotlib` for configuration of the resulting plot.
         """
-        neg_x = np.linspace(self.value-(num_sigma*self.minus),self.value,discretization)
-        pos_x = np.linspace(self.value,self.value+(num_sigma*self.minus),discretization)
-        p_neg = np.sqrt(2)/np.sqrt(np.pi)/(self.plus+self.minus) * np.exp(-1*(neg_x-self.value)**2 / (2*self.minus**2))
-        p_pos = np.sqrt(2)/np.sqrt(np.pi)/(self.plus+self.minus) * np.exp(-1*(pos_x-self.value)**2 / (2*self.plus**2))
+        neg_x = np.linspace(self.value-(num_sigma*self.minus), self.value, discretization)
+        pos_x = np.linspace(self.value, self.value+(num_sigma*self.plus), discretization)
         x = np.array(list(neg_x)+list(pos_x))
         pdf = self.pdf(x)
         plt.plot(x,pdf,**kwargs)
@@ -107,13 +114,10 @@ def cdfplot(self,num_sigma=5,discretization=100,**kwargs):
         Plots the associated CDF over the specified number of sigma, using 2*`discretization` points.
         `**kwargs` are passed on to `matplotlib` for configuration of the resulting plot.
         """
-        neg_x = np.linspace(self.value-(num_sigma*self.minus),self.value,discretization)
-        pos_x = np.linspace(self.value,self.value+(num_sigma*self.minus),discretization)
-        p_neg = np.sqrt(2)/np.sqrt(np.pi)/(self.plus+self.minus) * np.exp(-1*(neg_x-self.value)**2 / (2*self.minus**2))
-        p_pos = np.sqrt(2)/np.sqrt(np.pi)/(self.plus+self.minus) * np.exp(-1*(pos_x-self.value)**2 / (2*self.plus**2))
+        neg_x = np.linspace(self.value-(num_sigma*self.minus), self.value, discretization)
+        pos_x = np.linspace(self.value, self.value+(num_sigma*self.plus), discretization)
         x = np.array(list(neg_x)+list(pos_x))
-        pdf = self.pdf(x)
-        cdf = np.cumsum(pdf)/np.sum(pdf)
+        cdf = self.cdf(x)
         plt.plot(x,cdf,**kwargs)
         plt.show()
 
@@ -426,4 +430,4 @@ def neg_errors(array):
     Stand-alone function to return an array of the negative errors of an array of `a_u` objects.
     Functional equivalent to `UncertaintyArray(array).minus`.
     """
-    return [v.minus for v in array]
+    return [v.minus for v in array]

setup.py

@@ -5,7 +5,7 @@
 
 if __name__ == '__main__':
     setup(name="asymmetric_uncertainty",
-          version="0.2.2",
+          version="0.2.3",
           author="Caden Gobat",
           author_email="[email protected]",
           packages=find_packages(),