Changeset 779 for mystic


Ignore:
Timestamp:
01/12/15 16:44:03 (16 months ago)
Author:
mmckerns
Message:

add gamma and lgamma, update erf

File:
1 edited

Legend:

Unmodified
Added
Removed
  • mystic/_math/stats.py

    r776 r779  
    1616 
    1717def erf(x): 
    18   """Error function approximation.  
    19 Source: http://www.johndcook.com/python_erf.html 
    20 """ 
    21   # constants 
    22   a1 =  0.254829592 
    23   a2 = -0.284496736 
    24   a3 =  1.421413741 
    25   a4 = -1.453152027 
    26   a5 =  1.061405429 
    27   p  =  0.3275911 
    28  
    29   # Save the sign of x 
    30   sign = 1 
    31   if x < 0: 
    32     sign = -1 
    33   x = abs(x) 
    34  
    35   # A&S formula 7.1.26 
    36   t = 1.0/(1.0 + p*x) 
    37   y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*math.exp(-x*x) 
    38   return sign*y 
     18    "evaluate the error function at x"  
     19    try:                                                
     20        from math import erf 
     21    except ImportError: 
     22        try: 
     23            from ctypes import util, CDLL, c_double 
     24            erf = CDLL(util.find_library('m')).erf 
     25            erf.argtypes = [c_double] 
     26            erf.restype = c_double 
     27        except ImportError: 
     28            erf = _erf 
     29    if hasattr(x, '__len__'): 
     30        from numpy import vectorize 
     31        erf = vectorize(erf) 
     32    return erf(x) 
     33 
     34def gamma(x): 
     35    "evaluate the gamma function at x" 
     36    try: 
     37        from math import gamma 
     38    except ImportError: 
     39        try: 
     40            from ctypes import util, CDLL, c_double 
     41            gamma = CDLL(util.find_library('m')).gamma 
     42            gamma.argtypes = [c_double] 
     43            gamma.restype = c_double 
     44        except ImportError: 
     45            gamma = _gamma 
     46    if hasattr(x, '__len__'): 
     47        from numpy import vectorize 
     48        gamma = vectorize(gamma) 
     49    return gamma(x)  
     50 
     51def lgamma(x): 
     52    "evaluate the natual log of the abs value of the gamma function at x" 
     53    try: 
     54        from math import lgamma 
     55    except ImportError: 
     56        try: 
     57            from ctypes import util, CDLL, c_double 
     58            lgamma = CDLL(util.find_library('m')).lgamma 
     59            lgamma.argtypes = [c_double] 
     60            lgamma.restype = c_double 
     61        except ImportError: 
     62            lgamma = _lgamma 
     63    if hasattr(x, '__len__'): 
     64        from numpy import vectorize 
     65        lgamma = vectorize(lgamma) 
     66    return lgamma(x)  
     67 
     68 
     69def _erf(x): 
     70    "approximate the error function at x" 
     71    # Source: http://picomath.org/python/erf.py.html 
     72    a1 =  0.254829592 
     73    a2 = -0.284496736 
     74    a3 =  1.421413741 
     75    a4 = -1.453152027 
     76    a5 =  1.061405429 
     77    p  =  0.3275911 
     78 
     79    # Save the sign of x 
     80    sign = 1 
     81    if x < 0: 
     82        sign = -1 
     83    x = abs(x) 
     84 
     85    # A&S formula 7.1.26 
     86    t = 1.0/(1.0 + p*x) 
     87    y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*math.exp(-x*x) 
     88    return sign*y 
     89 
     90 
     91def _gamma(x): 
     92    "approximate the gamma function at x" 
     93    # Source: http://picomath.org/python/gamma.py.html 
     94    if x <= 0: 
     95        raise ValueError("Invalid input") 
     96 
     97    # Split the function domain into three intervals: 
     98    # (0, 0.001), [0.001, 12), and (12, infinity) 
     99 
     100    ########################################################################### 
     101    # First interval: (0, 0.001) 
     102    # 
     103    # For small x, 1/Gamma(x) has power series x + gamma x^2  - ... 
     104    # So in this range, 1/Gamma(x) = x + gamma x^2 with error order of x^3. 
     105    # The relative error over this interval is less than 6e-7. 
     106 
     107    gamma = 0.577215664901532860606512090 # Euler's gamma constant 
     108 
     109    if x < 0.001: 
     110        return 1.0/(x*(1.0 + gamma*x)) 
     111 
     112    ########################################################################### 
     113    # Second interval: [0.001, 12) 
     114 
     115    if x < 12.0: 
     116        # The algorithm directly approximates gamma over (1,2) and uses 
     117        # reduction identities to reduce other arguments to this interval. 
     118         
     119        y = x 
     120        n = 0 
     121        arg_was_less_than_one = (y < 1.0) 
     122 
     123        # Add or subtract integers as necessary to bring y into (1,2) 
     124        # Will correct for this below 
     125        if arg_was_less_than_one: 
     126            y += 1.0 
     127        else: 
     128            n = int(math.floor(y)) - 1  # will use n later 
     129            y -= n 
     130 
     131        # numerator coefficients for approximation over the interval (1,2) 
     132        p = [ 
     133            -1.71618513886549492533811E+0, 
     134             2.47656508055759199108314E+1, 
     135            -3.79804256470945635097577E+2, 
     136             6.29331155312818442661052E+2, 
     137             8.66966202790413211295064E+2, 
     138            -3.14512729688483675254357E+4, 
     139            -3.61444134186911729807069E+4, 
     140             6.64561438202405440627855E+4 
     141        ] 
     142 
     143        # denominator coefficients for approximation over the interval (1,2) 
     144        q = [ 
     145            -3.08402300119738975254353E+1, 
     146             3.15350626979604161529144E+2, 
     147            -1.01515636749021914166146E+3, 
     148            -3.10777167157231109440444E+3, 
     149             2.25381184209801510330112E+4, 
     150             4.75584627752788110767815E+3, 
     151            -1.34659959864969306392456E+5, 
     152            -1.15132259675553483497211E+5 
     153        ] 
     154 
     155        num = 0.0 
     156        den = 1.0 
     157 
     158        z = y - 1 
     159        for i in range(8): 
     160            num = (num + p[i])*z 
     161            den = den*z + q[i] 
     162        result = num/den + 1.0 
     163 
     164        # Apply correction if argument was not initially in (1,2) 
     165        if arg_was_less_than_one: 
     166            # Use identity gamma(z) = gamma(z+1)/z 
     167            # The variable "result" now holds gamma of the original y + 1 
     168            # Thus we use y-1 to get back the orginal y. 
     169            result /= (y-1.0) 
     170        else: 
     171            # Use the identity gamma(z+n) = z*(z+1)* ... *(z+n-1)*gamma(z) 
     172            for _ in range(n): 
     173                result *= y 
     174                y += 1 
     175 
     176        return result 
     177 
     178    ########################################################################### 
     179    # Third interval: [12, infinity) 
     180 
     181    if x > 171.624: 
     182        # Correct answer too large to display.  
     183        return 1.0/0 # float infinity 
     184 
     185    return math.exp(log_gamma(x)) 
     186 
     187 
     188def _lgamma(x): 
     189    "approximate the natual log of the abs value of the gamma function at x" 
     190    # Source: http://picomath.org/python/gamma.py.html 
     191    if x <= 0: 
     192        raise ValueError("Invalid input") 
     193 
     194    if x < 12.0: 
     195        return math.log(abs(gamma(x))) 
     196 
     197    # Abramowitz and Stegun 6.1.41 
     198    # Asymptotic series should be good to at least 11 or 12 figures 
     199    # For error analysis, see Whittiker and Watson 
     200    # A Course in Modern Analysis (1927), page 252 
     201 
     202    c = [ 
     203         1.0/12.0, 
     204        -1.0/360.0, 
     205         1.0/1260.0, 
     206        -1.0/1680.0, 
     207         1.0/1188.0, 
     208        -691.0/360360.0, 
     209         1.0/156.0, 
     210        -3617.0/122400.0 
     211    ] 
     212    z = 1.0/(x*x) 
     213    sum = c[7] 
     214    for i in range(6, -1, -1): 
     215        sum *= z 
     216        sum += c[i] 
     217    series = sum/x 
     218 
     219    halfLogTwoPi = 0.91893853320467274178032973640562 
     220    logGamma = (x - 0.5)*math.log(x) - x + halfLogTwoPi + series 
     221    return logGamma 
    39222 
    40223 
Note: See TracChangeset for help on using the changeset viewer.