Mercurial > hg > medcouple

comparison pymedcouple @ 23:29a178b23219

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Whitespace style fixes
author Jordi Gutiérrez Hermoso <jordigh@octave.org>
date Sat, 17 Jan 2015 20:04:46 -0500
parents f5b4a2ab6204
children cd3094a08e03
comparison
equal deleted inserted replaced
22:f5b4a2ab6204 23:29a178b23219
2 2
3 import random 3 import random
4 import numpy as np 4 import numpy as np
5 5
6 from itertools import tee, izip 6 from itertools import tee, izip
7
7 8
8 def partsort(L, n): 9 def partsort(L, n):
9 """This is a partial sort of the weighted list L = [(a,w) for a in A, 10 """This is a partial sort of the weighted list L = [(a,w) for a in A,
10 w in W], so that the element at position n is in the correct spot. 11 w in W], so that the element at position n is in the correct spot.
11 Also known as quickselect. The elements of W are ignored for the 12 Also known as quickselect. The elements of W are ignored for the
14 15
15 beg = 0 16 beg = 0
16 end = len(L) 17 end = len(L)
17 18
18 while True: 19 while True:
19 pivot = random.randint(beg,end-1) 20 pivot = random.randint(beg, end - 1)
20 pivotval = L[pivot][0] 21 pivotval = L[pivot][0]
21 L[pivot], L[end-1] = L[end-1], L[pivot] 22 L[pivot], L[end - 1] = L[end - 1], L[pivot]
22 23
23 idx = beg 24 idx = beg
24 for i in xrange(beg, end): 25 for i in xrange(beg, end):
25 if L[i][0] < pivotval: 26 if L[i][0] < pivotval:
26 L[idx], L[i] = L[i], L[idx] 27 L[idx], L[i] = L[i], L[idx]
29 L[idx], L[end-1] = L[end-1], L[idx] 30 L[idx], L[end-1] = L[end-1], L[idx]
30 31
31 if idx == n: 32 if idx == n:
32 return L[idx] 33 return L[idx]
33 elif idx < n: 34 elif idx < n:
34 beg = idx+1 35 beg = idx + 1
35 else: 36 else:
36 end = idx 37 end = idx
38
37 39
38 def wmedian(A, W): 40 def wmedian(A, W):
39 """This computes the weighted median of array A with corresponding 41 """This computes the weighted median of array A with corresponding
40 weights W. 42 weights W.
41 """ 43 """
42 44
43 AW = zip(A,W) 45 AW = zip(A, W)
44 n = len(AW) 46 n = len(AW)
45 47
46 wtot = sum(W) 48 wtot = sum(W)
47 49
48 beg = 0 50 beg = 0
49 end = n-1 51 end = n - 1
50 52
51 while True: 53 while True:
52 mid = (beg + end)//2 54 mid = (beg + end)//2
53 55
54 partsort(AW, mid) 56 partsort(AW, mid)
55 trial = AW[mid][0] 57 trial = AW[mid][0]
56 58
57 wleft = wright = 0 59 wleft = wright = 0
58 for (a,w) in AW: 60 for (a, w) in AW:
59 if a < trial: 61 if a < trial:
60 wleft += w 62 wleft += w
61 else: 63 else:
62 # This also includes a == trial, i.e. the "middle" 64 # This also includes a == trial, i.e. the "middle"
63 # weight. 65 # weight.
68 elif 2*wright < wtot: 70 elif 2*wright < wtot:
69 beg = mid 71 beg = mid
70 else: 72 else:
71 return trial 73 return trial
72 74
75
73 def medcouple_1d(X, eps1 = 2**-52, eps2 = 2**-1022): 76 def medcouple_1d(X, eps1 = 2**-52, eps2 = 2**-1022):
74 """Calculates the medcouple robust measure of skewness. 77 """Calculates the medcouple robust measure of skewness.
75 78
76 Parameters 79 Parameters
77 ---------- 80 ----------
92 95
93 """ 96 """
94 # FIXME: Figure out what to do about NaNs. 97 # FIXME: Figure out what to do about NaNs.
95 98
96 n = len(X) 99 n = len(X)
97 n2 = (n-1)//2 100 n2 = (n - 1)//2
98 101
99 if n < 3: 102 if n < 3:
100 return 0 103 return 0
101 104
102 Z = sorted(X, reverse=True) 105 Z = sorted(X, reverse=True)
103 106
104 if n % 2 == 1: 107 if n % 2 == 1:
105 Zmed = Z[n2] 108 Zmed = Z[n2]
106 else: 109 else:
107 Zmed = (Z[n2] + Z[n2+1])/2 110 Zmed = (Z[n2] + Z[n2 + 1])/2
108 111
109 #Check if the median is at the edges up to relative epsilon 112 #Check if the median is at the edges up to relative epsilon
110 if abs(Z[0] - Zmed) < eps1*(eps1 + abs(Zmed)): 113 if abs(Z[0] - Zmed) < eps1*(eps1 + abs(Zmed)):
111 return -1.0 114 return -1.0
112 if abs(Z[-1] - Zmed) < eps1*(eps1 + abs(Zmed)): 115 if abs(Z[-1] - Zmed) < eps1*(eps1 + abs(Zmed)):
127 Zminus = [z for z in Z if Zeps >= z] 130 Zminus = [z for z in Z if Zeps >= z]
128 131
129 n_plus = len(Zplus) 132 n_plus = len(Zplus)
130 n_minus = len(Zminus) 133 n_minus = len(Zminus)
131 134
135
132 def h_kern(i, j): 136 def h_kern(i, j):
133 """Kernel function h for the medcouple, closing over the values of 137 """Kernel function h for the medcouple, closing over the values of
134 Zplus and Zminus just defined above. 138 Zplus and Zminus just defined above.
135 139
136 In case a and be are within epsilon of the median, the kernel 140 In case a and be are within epsilon of the median, the kernel
141 b = Zminus[j] 145 b = Zminus[j]
142 146
143 if abs(a - b) <= 2*eps2: 147 if abs(a - b) <= 2*eps2:
144 h = signum(n_plus - 1 - i - j) 148 h = signum(n_plus - 1 - i - j)
145 else: 149 else:
146 h = (a+b)/(a-b) 150 h = (a + b)/(a - b)
147 151
148 return h 152 return h
149 153
150 # Init left and right borders 154 # Init left and right borders
151 L = [0]*n_plus 155 L = [0]*n_plus
162 # (Be stingy, reuse same generator) 166 # (Be stingy, reuse same generator)
163 [I1, I2] = tee(i for i in xrange(0, n_plus) if L[i] <= R[i]) 167 [I1, I2] = tee(i for i in xrange(0, n_plus) if L[i] <= R[i])
164 168
165 A = [h_kern(i, (L[i] + R[i])//2) for i in I1] 169 A = [h_kern(i, (L[i] + R[i])//2) for i in I1]
166 W = [R[i] - L[i] + 1 for i in I2] 170 W = [R[i] - L[i] + 1 for i in I2]
167 Am = wmedian(A,W) 171 Am = wmedian(A, W)
168 172
169 Am_eps = eps1*(eps1 + abs(Am)) 173 Am_eps = eps1*(eps1 + abs(Am))
170 174
171 # Compute new left and right boundaries, based on the weighted 175 # Compute new left and right boundaries, based on the weighted
172 # median 176 # median
210 for j in xrange(l, r + 1): 214 for j in xrange(l, r + 1):
211 A.append(h_kern(i, j)) 215 A.append(h_kern(i, j))
212 A.sort() 216 A.sort()
213 A.reverse() 217 A.reverse()
214 218
215 Am = A[medc_idx - Ltot] 219 Am = A[medc_idx - Ltot]
216 220
217 return Am 221 return Am
218 222
219 223
220 def signum(x): 224 def signum(x):