Changeset 779
- Timestamp:
- 01/12/15 16:44:03 (16 months ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
mystic/_math/stats.py
r776 r779 16 16 17 17 def 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 34 def 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 51 def 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 69 def _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 91 def _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 188 def _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 39 222 40 223
Note: See TracChangeset
for help on using the changeset viewer.