annotate corebio/moremath.py @ 8:5149eb3a89c2

Uploaded
author davidmurphy
date Fri, 20 Jan 2012 09:03:40 -0500
parents c55bdc2fb9fa
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
1 #!/usr/bin/env python
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
2
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
3 # Copyright (c) 2005 Gavin E. Crooks <gec@threeplusone.com>
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
4 #
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
5 # This software is distributed under the MIT Open Source License.
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
6 # <http://www.opensource.org/licenses/mit-license.html>
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
7 #
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
8 # Permission is hereby granted, free of charge, to any person obtaining a
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
9 # copy of this software and associated documentation files (the "Software"),
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
10 # to deal in the Software without restriction, including without limitation
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
11 # the rights to use, copy, modify, merge, publish, distribute, sublicense,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
12 # and/or sell copies of the Software, and to permit persons to whom the
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
13 # Software is furnished to do so, subject to the following conditions:
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
14 #
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
15 # The above copyright notice and this permission notice shall be included
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
16 # in all copies or substantial portions of the Software.
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
17 #
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
18 # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
19 # IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
20 # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
21 # AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
22 # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
23 # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
24 # THE SOFTWARE.
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
25 #
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
26
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
27
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
28 """ Various bits of useful math not in the standard python library.
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
29
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
30 Constants :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
31
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
32 - euler_gamma = 0.577215...
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
33 - catalan = 0.915965...
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
34 - golden_ratio = 1.618033...
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
35 - bits_per_nat = log2(e) = 1/log(2)
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
36 - sqrt_2pi = 2.50662...
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
37
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
38 Special Functions :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
39
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
40 - gamma() -- Gamma function.
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
41 - lngamma() -- Logarithm of the gamma function
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
42 - factorial() -- The factorial function.
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
43 - digamma() -- Digamma function (logarithmic derivative of gamma).
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
44 - trigamma() -- Trigamma function (derivative of digamma).
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
45 - entropy() -- The entropy of a probability vector
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
46 - incomplete_gamma() -- The 'upper' incomplete gamma function.
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
47 - normalized_incomplete_gamma() --
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
48 - lg() -- Base 2 logarithms.
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
49
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
50
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
51 Vector Operations :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
52
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
53 - rmsd() -- Root mean squared deviation of two point vectors
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
54 - minimize_rmsd() -- Find the rigid transformation that minimized the
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
55 RMSD between two vectors of points.
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
56
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
57 Minimization :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
58
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
59 - find_root() -- 1d root finding
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
60
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
61 Probability Distributions :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
62 - Gamma
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
63 - Dirichlet
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
64 - Multinomial
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
65 - Gaussian
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
66
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
67 """
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
68
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
69
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
70
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
71 __all__ = ('euler_gamma', 'catalan', 'golden_ratio', 'bits_per_nat', 'sqrt_2pi',
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
72 'gamma', 'lngamma', 'factorial', 'digamma', 'trigamma',
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
73 'entropy', 'log2',
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
74 'incomplete_gamma', 'normalized_incomplete_gamma',
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
75 # 'integrate',
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
76 # 'rmsd', 'minimize_rmsd', 'find_root',
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
77 # 'Gamma', 'Dirichlet',
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
78 # 'decompose_log_odds_array',
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
79 'argmax', 'argmin'
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
80 )
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
81
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
82 from math import *
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
83 import random
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
84 from itertools import izip, count
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
85
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
86 # Some mathematical constants
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
87 euler_gamma = 0.57721566490153286060651
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
88 catalan = 0.91596559417721901505460
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
89 golden_ratio = 1.6180339887498948482046
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
90 bits_per_nat = 1.44269504088896340735992468100 # = log_2(e) = 1/log(2)
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
91 sqrt_2pi = 2.5066282746310005024157652848110
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
92
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
93
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
94
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
95
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
96
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
97 # The Lanczos approximation for the gamma function is
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
98 #
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
99 # -(z + g + 1/2) (z + 1/2)
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
100 # Gamma(z+1) = e * (z + g + 1/2) * Sqrt(2Pi) * C
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
101 #
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
102 #
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
103 # c[1] c[2] c[3]
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
104 # C = [c[0] + ----- + ----- + ----- + ... ]
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
105 # z + 1 z + 2 z + 3
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
106 #
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
107 #
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
108 # To calculate digamma and trigamma functions we take an analytic derivative
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
109 # of the Lanczos approximation.
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
110 #
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
111 # Gamma(z) = Gamma(z+1)/z
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
112 # Digamma(z) = D ln Gamma(z)
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
113 # Trigamma(z) = D Digamma(z)
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
114
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
115 # These Lanczos constants are from
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
116 # "A note on the computation of the convergent
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
117 # Lanczos complex Gamma approximation." Paul Godfrey (2001)
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
118 # http://my.fit.edu/~gabdo/gamma.txt
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
119
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
120
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
121 __lanczos_gamma = 607./128.
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
122 __lanczos_coefficients = (
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
123 0.99999999999999709182,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
124 57.156235665862923517,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
125 -59.597960355475491248,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
126 14.136097974741747174,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
127 -0.49191381609762019978,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
128 .33994649984811888699e-4,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
129 .46523628927048575665e-4,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
130 -.98374475304879564677e-4,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
131 .15808870322491248884e-3,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
132 -.21026444172410488319e-3,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
133 .21743961811521264320e-3,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
134 -.16431810653676389022e-3,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
135 .84418223983852743293e-4,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
136 -.26190838401581408670e-4,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
137 .36899182659531622704e-5)
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
138
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
139 __factorial =(
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
140 1.,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
141 1.,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
142 2.,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
143 6.,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
144 24.,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
145 120.,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
146 720.,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
147 5040.,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
148 40320.,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
149 362880.,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
150 3628800.,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
151 39916800.,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
152 479001600.,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
153 6227020800.,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
154 87178291200.,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
155 1307674368000.,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
156 20922789888000.,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
157 355687428096000.,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
158 6402373705728000.,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
159 121645100408832000.,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
160 2432902008176640000.,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
161 51090942171709440000.,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
162 1124000727777607680000.,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
163 25852016738884976640000.,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
164 620448401733239439360000.,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
165 15511210043330985984000000.,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
166 403291461126605635584000000.,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
167 10888869450418352160768000000.,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
168 304888344611713860501504000000.,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
169 8841761993739701954543616000000.,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
170 265252859812191058636308480000000.,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
171 8222838654177922817725562880000000.,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
172 263130836933693530167218012160000000. )
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
173
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
174 def gamma(z) :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
175 """The gamma function. Returns exact results for small integers. Will
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
176 overflow for modest sized arguments. Use lngamma(z) instead.
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
177
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
178 See: Eric W. Weisstein. "Gamma Function." From MathWorld, A Wolfram Web Resource.
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
179 http://mathworld.wolfram.com/GammaFunction.html
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
180
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
181 """
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
182
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
183 n = floor(z)
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
184 if n == z :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
185 if z <= 0 :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
186 return 1.0/0.0 # Infinity
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
187 elif n <= len(__factorial) :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
188 return __factorial[int(n)-1]
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
189
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
190 zz = z
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
191 if z < 0.5 :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
192 zz = 1-z
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
193
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
194
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
195 g = __lanczos_gamma
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
196 c = __lanczos_coefficients
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
197
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
198 zz = zz - 1.
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
199 zh = zz + 0.5
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
200 zgh = zh + g
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
201 zp = zgh** (zh*0.5) # trick for avoiding FP overflow above z=141
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
202
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
203 ss = 0.0
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
204 for k in range(len(c)-1,0,-1):
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
205 ss += c[k]/(zz+k)
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
206
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
207 f = (sqrt_2pi*(c[0]+ss)) * (( zp*exp(-zgh)) *zp)
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
208
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
209 if z<0.5 :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
210 f = pi /( sin(pi*z) *f)
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
211
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
212 return f
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
213
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
214
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
215 def lngamma(z) :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
216 """The logarithm of the gamma function.
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
217 """
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
218
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
219 # common case optimization
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
220
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
221 n = floor(z)
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
222 if n == z :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
223 if z <= 0 :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
224 return 1.0/0.0 # Infinity
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
225 elif n <= len(__factorial) :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
226 return __factorial[int(n)-1]
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
227
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
228 zz = z
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
229 if z < 0.5 :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
230 zz = 1-z
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
231
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
232
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
233 g = __lanczos_gamma
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
234 c = __lanczos_coefficients
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
235
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
236 zz = zz - 1.
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
237 zh = zz + 0.5
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
238 zgh = zh + g
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
239 zp = zgh** (zh*0.5) # trick for avoiding FP overflow above z=141
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
240
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
241 ss = 0.0
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
242 for k in range(len(c)-1,0,-1):
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
243 ss += c[k]/(zz+k)
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
244
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
245 f = (sqrt_2pi*(c[0]+ss)) * (( zp*exp(-zgh)) *zp)
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
246
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
247 if z<0.5 :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
248 f = pi /( sin(pi*z) *f)
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
249
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
250 return log(f)
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
251
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
252
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
253 def factorial(z) :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
254 """ The factorial function.
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
255 factorial(z) == gamma(z+1)
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
256 """
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
257 return gamma(z+1)
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
258
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
259
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
260 def digamma(z) :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
261 """The digamma function, the logarithmic derivative of the gamma function.
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
262 digamma(z) = d/dz ln( gamma(z) )
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
263
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
264 See: Eric W. Weisstein. "Digamma Function." From MathWorld--
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
265 A Wolfram Web Resource. http://mathworld.wolfram.com/DigammaFunction.html
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
266 """
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
267
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
268 g = __lanczos_gamma
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
269 c = __lanczos_coefficients
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
270
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
271 zz = z
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
272 if z < 0.5 :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
273 zz = 1 -z
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
274
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
275 n=0.
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
276 d=0.
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
277 for k in range(len(c)-1,0,-1):
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
278 dz =1./(zz+(k+1)-2);
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
279 dd =c[k] * dz
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
280 d = d + dd
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
281 n = n - dd * dz
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
282
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
283 d = d + c[0]
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
284 gg = zz + g - 0.5
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
285 f = log(gg) + (n/d - g/gg)
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
286
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
287 if z<0.5 :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
288 f -= pi / tan( pi * z)
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
289
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
290 return f
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
291
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
292
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
293 def trigamma(z) :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
294 """The trigamma function, the derivative of the digamma function.
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
295 trigamma(z) = d/dz digamma(z) = d/dz d/dz ln( gamma(z) )
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
296
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
297 See: Eric W. Weisstein. "Digamma Function." From MathWorld--
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
298 A Wolfram Web Resource. http://mathworld.wolfram.com/TrigammaFunction.html
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
299 """
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
300
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
301 g = __lanczos_gamma
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
302 c = __lanczos_coefficients
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
303
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
304 t1=0.
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
305 t2=0.
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
306 t3=0.
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
307 for k in range(len(c)-1,0,-1):
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
308 dz =1./(z+k);
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
309 dd1 = c[k]* dz
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
310 t1 += dd1
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
311 dd2 = dd1 * dz
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
312 t2 += dd2
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
313 t3 += dd2 * dz
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
314
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
315 t1 += c[0]
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
316 c = - (t2*t2)/(t1*t1) +2*t3/t1
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
317
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
318 result = 1./(z*z)
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
319 gg = z + g + 0.5
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
320 result += - (z+0.5)/ (gg*gg)
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
321 result += 2./gg
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
322
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
323 result += c
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
324
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
325 return result
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
326
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
327 def incomplete_gamma(a,x) :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
328 """The 'upper' incomplete gamma function:
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
329
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
330 oo
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
331 -
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
332 | -t a-1
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
333 incomplete_gamma(a,x) = | e t dt.
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
334 |
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
335 -
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
336 x
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
337
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
338 In Mathematica, Gamma[a,x].
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
339
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
340 Note that, very confusingly, the phrase 'incomplete gamma fucntion'
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
341 can also refer to the same integral between 0 and x, (the 'lower'
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
342 incomplete gamma function) or to the normalized versions,
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
343 normalized_incomplete_gamma() )
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
344
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
345
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
346 See: Eric W. Weisstein. "Gamma Function." From MathWorld, A Wolfram Web Resource.
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
347 http://mathworld.wolfram.com/IncompleteGammaFunction.html
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
348
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
349 Bugs :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
350 This implentation is not very accurate for some arguments.
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
351 """
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
352 return normalized_incomplete_gamma(a,x) * gamma(a)
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
353
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
354
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
355 def normalized_incomplete_gamma(a,x) :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
356 """The upper, incomplete gamma function normalized so that the limiting
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
357 values are zero and one.
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
358
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
359 Q(a,x) = incomplete_gamma(a,x) / gamma(a)
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
360
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
361 See:
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
362 incomplete_gamma()
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
363 Bugs :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
364 This implentation is not very accurate for some arguments.
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
365 """
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
366 maxiter = 100
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
367 epsilon = 1.48e-8
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
368 small = 1e-30
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
369
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
370
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
371 if a<=0 or x<0 :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
372 raise ValueError("Invalid arguments")
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
373 if x == 0.0 : return 1.0
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
374
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
375 if x<= a+1 :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
376 # Use the series representation
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
377 term = 1./a
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
378 total = term
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
379 for n in range(1,maxiter) :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
380 term *= x/(a+n)
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
381 total += term
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
382 if abs(term/total) < epsilon :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
383 return 1. - total * exp(-x+a*log(x) - lngamma(a) )
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
384 raise RuntimeError(
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
385 "Failed to converge after %d iterations." % (maxiter) )
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
386 else :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
387 # Use the continued fraction representation
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
388 total = 1.0
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
389 b = x + 1. -a
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
390 c = 1./small
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
391 d = 1./b
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
392 h = d
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
393 for i in range(1, maxiter) :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
394 an = -i * (i-a)
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
395 b = b+2.
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
396 d = an * d + b
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
397 if abs(d) < small : d = small
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
398 c = b + an /c
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
399 if abs(c) < small : c= small
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
400 d = 1./d
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
401 term = d * c
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
402 h = h * term
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
403 if abs( term-1.) < epsilon :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
404 return h * exp(-x+a*log(x) - lngamma(a) )
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
405 raise RuntimeError(
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
406 "Failed to converge after %d iterations." % (maxiter) )
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
407
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
408
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
409
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
410 def log2( x) :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
411 """ Return the base 2 logarithm of x """
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
412 return log(x,2)
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
413
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
414
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
415 def entropy( pvec, base= exp(1) ) :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
416 """ The entropy S = -Sum_i p_i ln p_i
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
417 pvec is a frequency vector, not necessarily normalized.
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
418 """
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
419 # TODO: Optimize
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
420 if len(pvec) ==0 :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
421 raise ValueError("Zero length vector")
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
422
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
423
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
424 total = 0.0
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
425 ent = 0.0
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
426 for p in pvec:
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
427 if p>0 : # 0 log(0) =0
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
428 total += p
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
429 ent += - log(float(p)) *p
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
430 elif p<0:
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
431 raise ValueError("Negative probability")
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
432
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
433
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
434 ent = (ent/total) + log(total)
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
435 ent /= log(base)
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
436
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
437 return ent
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
438
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
439
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
440
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
441
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
442
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
443 def argmax( alist) :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
444 """Return the index of the last occurance of the maximum value in the list."""
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
445 return max(izip(alist, count() ))[1]
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
446
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
447 def argmin( alist) :
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
448 """Return the index of the first occurance of the minimum value in the list."""
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
449 return min(izip(alist, count() ))[1]
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
450
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
451
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
452
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
453
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
454
c55bdc2fb9fa Uploaded
davidmurphy
parents:
diff changeset
455