Changeset 779

Timestamp:

01/12/15 16:44:03 (16 months ago)

Author:

mmckerns

Message:

add gamma and lgamma, update erf

File:

• mystic/_math/stats.py (modified) (1 diff)

mystic/_math/stats.py

@@ -16 +16 @@
 
 def erf(x):
-  """Error function approximation. 
-Source: http://www.johndcook.com/python_erf.html
-"""
-  # constants
-  a1 =  0.254829592
-  a2 = -0.284496736
-  a3 =  1.421413741
-  a4 = -1.453152027
-  a5 =  1.061405429
-  p  =  0.3275911
-
-  # Save the sign of x
-  sign = 1
-  if x < 0:
-    sign = -1
-  x = abs(x)
-
-  # A&S formula 7.1.26
-  t = 1.0/(1.0 + p*x)
-  y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*math.exp(-x*x)
-  return sign*y
+    "evaluate the error function at x" 
+    try:                                               
+        from math import erf
+    except ImportError:
+        try:
+            from ctypes import util, CDLL, c_double
+            erf = CDLL(util.find_library('m')).erf
+            erf.argtypes = [c_double]
+            erf.restype = c_double
+        except ImportError:
+            erf = _erf
+    if hasattr(x, '__len__'):
+        from numpy import vectorize
+        erf = vectorize(erf)
+    return erf(x)
+
+def gamma(x):
+    "evaluate the gamma function at x"
+    try:
+        from math import gamma
+    except ImportError:
+        try:
+            from ctypes import util, CDLL, c_double
+            gamma = CDLL(util.find_library('m')).gamma
+            gamma.argtypes = [c_double]
+            gamma.restype = c_double
+        except ImportError:
+            gamma = _gamma
+    if hasattr(x, '__len__'):
+        from numpy import vectorize
+        gamma = vectorize(gamma)
+    return gamma(x) 
+
+def lgamma(x):
+    "evaluate the natual log of the abs value of the gamma function at x"
+    try:
+        from math import lgamma
+    except ImportError:
+        try:
+            from ctypes import util, CDLL, c_double
+            lgamma = CDLL(util.find_library('m')).lgamma
+            lgamma.argtypes = [c_double]
+            lgamma.restype = c_double
+        except ImportError:
+            lgamma = _lgamma
+    if hasattr(x, '__len__'):
+        from numpy import vectorize
+        lgamma = vectorize(lgamma)
+    return lgamma(x) 
+
+
+def _erf(x):
+    "approximate the error function at x"
+    # Source: http://picomath.org/python/erf.py.html
+    a1 =  0.254829592
+    a2 = -0.284496736
+    a3 =  1.421413741
+    a4 = -1.453152027
+    a5 =  1.061405429
+    p  =  0.3275911
+
+    # Save the sign of x
+    sign = 1
+    if x < 0:
+        sign = -1
+    x = abs(x)
+
+    # A&S formula 7.1.26
+    t = 1.0/(1.0 + p*x)
+    y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*math.exp(-x*x)
+    return sign*y
+
+
+def _gamma(x):
+    "approximate the gamma function at x"
+    # Source: http://picomath.org/python/gamma.py.html
+    if x <= 0:
+        raise ValueError("Invalid input")
+
+    # Split the function domain into three intervals:
+    # (0, 0.001), [0.001, 12), and (12, infinity)
+
+    ###########################################################################
+    # First interval: (0, 0.001)
+    #
+    # For small x, 1/Gamma(x) has power series x + gamma x^2  - ...
+    # So in this range, 1/Gamma(x) = x + gamma x^2 with error order of x^3.
+    # The relative error over this interval is less than 6e-7.
+
+    gamma = 0.577215664901532860606512090 # Euler's gamma constant
+
+    if x < 0.001:
+        return 1.0/(x*(1.0 + gamma*x))
+
+    ###########################################################################
+    # Second interval: [0.001, 12)
+
+    if x < 12.0:
+        # The algorithm directly approximates gamma over (1,2) and uses
+        # reduction identities to reduce other arguments to this interval.
+        
+        y = x
+        n = 0
+        arg_was_less_than_one = (y < 1.0)
+
+        # Add or subtract integers as necessary to bring y into (1,2)
+        # Will correct for this below
+        if arg_was_less_than_one:
+            y += 1.0
+        else:
+            n = int(math.floor(y)) - 1  # will use n later
+            y -= n
+
+        # numerator coefficients for approximation over the interval (1,2)
+        p = [
+            -1.71618513886549492533811E+0,
+             2.47656508055759199108314E+1,
+            -3.79804256470945635097577E+2,
+             6.29331155312818442661052E+2,
+             8.66966202790413211295064E+2,
+            -3.14512729688483675254357E+4,
+            -3.61444134186911729807069E+4,
+             6.64561438202405440627855E+4
+        ]
+
+        # denominator coefficients for approximation over the interval (1,2)
+        q = [
+            -3.08402300119738975254353E+1,
+             3.15350626979604161529144E+2,
+            -1.01515636749021914166146E+3,
+            -3.10777167157231109440444E+3,
+             2.25381184209801510330112E+4,
+             4.75584627752788110767815E+3,
+            -1.34659959864969306392456E+5,
+            -1.15132259675553483497211E+5
+        ]
+
+        num = 0.0
+        den = 1.0
+
+        z = y - 1
+        for i in range(8):
+            num = (num + p[i])*z
+            den = den*z + q[i]
+        result = num/den + 1.0
+
+        # Apply correction if argument was not initially in (1,2)
+        if arg_was_less_than_one:
+            # Use identity gamma(z) = gamma(z+1)/z
+            # The variable "result" now holds gamma of the original y + 1
+            # Thus we use y-1 to get back the orginal y.
+            result /= (y-1.0)
+        else:
+            # Use the identity gamma(z+n) = z*(z+1)* ... *(z+n-1)*gamma(z)
+            for _ in range(n):
+                result *= y
+                y += 1
+
+        return result
+
+    ###########################################################################
+    # Third interval: [12, infinity)
+
+    if x > 171.624:
+        # Correct answer too large to display. 
+        return 1.0/0 # float infinity
+
+    return math.exp(log_gamma(x))
+
+
+def _lgamma(x):
+    "approximate the natual log of the abs value of the gamma function at x"
+    # Source: http://picomath.org/python/gamma.py.html
+    if x <= 0:
+        raise ValueError("Invalid input")
+
+    if x < 12.0:
+        return math.log(abs(gamma(x)))
+
+    # Abramowitz and Stegun 6.1.41
+    # Asymptotic series should be good to at least 11 or 12 figures
+    # For error analysis, see Whittiker and Watson
+    # A Course in Modern Analysis (1927), page 252
+
+    c = [
+         1.0/12.0,
+        -1.0/360.0,
+         1.0/1260.0,
+        -1.0/1680.0,
+         1.0/1188.0,
+        -691.0/360360.0,
+         1.0/156.0,
+        -3617.0/122400.0
+    ]
+    z = 1.0/(x*x)
+    sum = c[7]
+    for i in range(6, -1, -1):
+        sum *= z
+        sum += c[i]
+    series = sum/x
+
+    halfLogTwoPi = 0.91893853320467274178032973640562
+    logGamma = (x - 0.5)*math.log(x) - x + halfLogTwoPi + series
+    return logGamma