Mercurial > hg > medcouple

comparison pymedcouple @ 13:077261db7a58

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Whitespace cleanup
author Jordi Gutiérrez Hermoso <jordigh@octave.org>
date Fri, 16 Jan 2015 09:27:44 -0500
parents c81f6d263897
children 8e02346789b2
comparison
equal deleted inserted replaced
12:c81f6d263897 13:077261db7a58
9 """This is a partial sort of the weighted list L = [(a,w) for a in A, 9 """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. 10 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 11 Also known as quickselect. The elements of W are ignored for the
12 purposes of this partial sort. 12 purposes of this partial sort.
13 """ 13 """
14 14
15 beg = 0 15 beg = 0
16 end = len(L) 16 end = len(L)
17 17
18 while True: 18 while True:
19 pivot = random.randint(beg,end-1) 19 pivot = random.randint(beg,end-1)
20 pivotval = L[pivot][0] 20 pivotval = L[pivot][0]
21 L[pivot], L[end-1] = L[end-1], L[pivot] 21 L[pivot], L[end-1] = L[end-1], L[pivot]
22 22
23 idx = beg 23 idx = beg
24 for i in xrange(beg, end): 24 for i in xrange(beg, end):
25 if L[i][0] < pivotval: 25 if L[i][0] < pivotval:
26 L[idx], L[i] = L[i], L[idx] 26 L[idx], L[i] = L[i], L[idx]
27 idx += 1 27 idx += 1
37 37
38 def wmedian(A, W): 38 def wmedian(A, W):
39 """This computes the weighted median of array A with corresponding 39 """This computes the weighted median of array A with corresponding
40 weights W. 40 weights W.
41 """ 41 """
42 42
43 AW = zip(A,W) 43 AW = zip(A,W)
44 n = len(AW) 44 n = len(AW)
45 45
46 wtot = sum(W) 46 wtot = sum(W)
47 47
48 beg = 0 48 beg = 0
49 end = n-1 49 end = n-1
50 50
51 while True: 51 while True:
52 mid = (beg + end)//2 52 mid = (beg + end)//2
53 53
54 partsort(AW, mid) 54 partsort(AW, mid)
55 trial = AW[mid][0] 55 trial = AW[mid][0]
56 56
57 wleft = wright = 0 57 wleft = wright = 0
58 for (a,w) in AW: 58 for (a,w) in AW:
79 79
80 Returns 80 Returns
81 ------- 81 -------
82 mc : float 82 mc : float
83 The medcouple statistic 83 The medcouple statistic
84 84
85 .. [1] G. Brys, M. Hubert, and A. Struyf "A Robust Measure of 85 .. [1] G. Brys, M. Hubert, and A. Struyf "A Robust Measure of
86 Skewness." Journal of Computational and Graphical Statistics, Vol. 86 Skewness." Journal of Computational and Graphical Statistics, Vol.
87 13, No. 4 (Dec., 2004), pp. 996- 1017 87 13, No. 4 (Dec., 2004), pp. 996- 1017
88 88
89 .. [2] D. B. Johnson and T. Mizoguchi "Selecting the Kth Element 89 .. [2] D. B. Johnson and T. Mizoguchi "Selecting the Kth Element
90 in $X + Y$ and $X_1 + X_2 + \cdots X_m$". SIAM Journal of 90 in $X + Y$ and $X_1 + X_2 + \cdots X_m$". SIAM Journal of
91 Computing, Vol. 7, No. 2 (May 1978), pp. 147-53. 91 Computing, Vol. 7, No. 2 (May 1978), pp. 147-53.
92 92
93 """ 93 """
94 # FIXME: Figure out what to do about NaNs. 94 # FIXME: Figure out what to do about NaNs.
95 95
96 n = len(X) 96 n = len(X)
97 n2 = (n-1)//2 97 n2 = (n-1)//2
98 98
99 if n < 3: 99 if n < 3:
100 return 0 100 return 0
103 103
104 if n % 2 == 1: 104 if n % 2 == 1:
105 Zmed = Z[n2] 105 Zmed = Z[n2]
106 else: 106 else:
107 Zmed = (Z[n2] + Z[n2+1])/2 107 Zmed = (Z[n2] + Z[n2+1])/2
108 108
109 #Check if the median is at the edges up to relative epsilon 109 #Check if the median is at the edges up to relative epsilon
110 if abs(Z[0] - Zmed) < eps1*(eps1 + abs(Zmed)): 110 if abs(Z[0] - Zmed) < eps1*(eps1 + abs(Zmed)):
111 return -1.0 111 return -1.0
112 if abs(Z[-1] - Zmed) < eps1*(eps1 + abs(Zmed)): 112 if abs(Z[-1] - Zmed) < eps1*(eps1 + abs(Zmed)):
113 return 1.0 113 return 1.0
117 117
118 # Scale inside [-0.5, 0.5], for greater numerical stability. 118 # Scale inside [-0.5, 0.5], for greater numerical stability.
119 Zden = 2*max(Z[0], -Z[-1]) 119 Zden = 2*max(Z[0], -Z[-1])
120 Z = [z/Zden for z in Z] 120 Z = [z/Zden for z in Z]
121 Zmed /= Zden 121 Zmed /= Zden
122 122
123 Zeps = eps1*(eps1 + abs(Zmed)) 123 Zeps = eps1*(eps1 + abs(Zmed))
124 124
125 # These overlap on the entries that are tied with the median 125 # These overlap on the entries that are tied with the median
126 Zplus = [z for z in Z if z >= -Zeps] 126 Zplus = [z for z in Z if z >= -Zeps]
127 Zminus = [z for z in Z if Zeps >= z] 127 Zminus = [z for z in Z if Zeps >= z]
155 Rtot = n_minus*n_plus 155 Rtot = n_minus*n_plus
156 medc_idx = (Rtot - 1)//2 156 medc_idx = (Rtot - 1)//2
157 157
158 # kth pair algorithm (Johnson & Mizoguchi) 158 # kth pair algorithm (Johnson & Mizoguchi)
159 while Rtot - Ltot > n_plus: 159 while Rtot - Ltot > n_plus:
160 160
161 # First, compute the median inside the given bounds 161 # First, compute the median inside the given bounds
162 # (Be stingy, reuse same generator) 162 # (Be stingy, reuse same generator)
163 [I1, I2] = tee(i for i in xrange(0, n_plus) if L[i] <= R[i]) 163 [I1, I2] = tee(i for i in xrange(0, n_plus) if L[i] <= R[i])
164 164
165 A = [h_kern(i, (L[i] + R[i])//2) for i in I1] 165 A = [h_kern(i, (L[i] + R[i])//2) for i in I1]
166 W = [R[i] - L[i] + 1 for i in I2] 166 W = [R[i] - L[i] + 1 for i in I2]
167 Am = wmedian(A,W) 167 Am = wmedian(A,W)
168 168
169 Am_eps = eps1*(eps1 + abs(Am)) 169 Am_eps = eps1*(eps1 + abs(Am))
170 170
171 # Compute new left and right boundaries, based on the weighted 171 # Compute new left and right boundaries, based on the weighted
172 # median 172 # median
173 P = [] 173 P = []
176 j = 0 176 j = 0
177 for i in xrange(n_plus - 1, -1, -1): 177 for i in xrange(n_plus - 1, -1, -1):
178 while j < n_minus - 1 and h_kern(i, j) - Am > Am_eps: 178 while j < n_minus - 1 and h_kern(i, j) - Am > Am_eps:
179 j += 1 179 j += 1
180 P.append(j - 1) 180 P.append(j - 1)
181 181
182 P.reverse() 182 P.reverse()
183 183
184 j = n_minus - 1 184 j = n_minus - 1
185 for i in xrange(0, n_plus): 185 for i in xrange(0, n_plus):
186 while j >= 0 and h_kern(i, j) - Am < -Am_eps: 186 while j >= 0 and h_kern(i, j) - Am < -Am_eps:
213 A.reverse() 213 A.reverse()
214 214
215 Am = A[medc_idx - Ltot] 215 Am = A[medc_idx - Ltot]
216 216
217 return Am 217 return Am
218 218
219 219
220 def signum(x): 220 def signum(x):
221 if x > 0: 221 if x > 0:
222 return 1 222 return 1
223 if x < 0: 223 if x < 0: