-
Notifications
You must be signed in to change notification settings - Fork 1
/
pav.py
79 lines (70 loc) · 2.13 KB
/
pav.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
import numpy as np
import bisect
import math
from adversary import Adversary
# Author : Alexandre Gramfort
# license : BSD
def pav(y):
"""
PAV uses the pair adjacent violators method to produce a monotonic
smoothing of y
translated from matlab by Sean Collins (2006) as part of the EMAP toolbox
"""
y = np.asarray(y)
assert y.ndim == 1
n_samples = len(y)
v = y.copy()
lvls = np.arange(n_samples)
lvlsets = np.c_[lvls, lvls]
flag = 1
while flag:
deriv = np.diff(v)
if np.all(deriv >= 0):
break
viol = np.where(deriv < 0)[0]
start = lvlsets[viol[0], 0]
last = lvlsets[viol[0] + 1, 1]
s = 0
n = last - start + 1
for i in range(start, last + 1):
s += v[i]
val = s / n
for i in range(start, last + 1):
v[i] = val
lvlsets[i, 0] = start
lvlsets[i, 1] = last
return v
if __name__ == '__main__':
T = 200
D = 5
R = 5
adv = Adversary(D, T, R)
E = 0
for t in range(1,T):
X, Z, Y = adv.train_data[0][:t], adv.train_data[1][:t], adv.train_data[2][:t]
sZ, sY = zip(*sorted(zip(Z,Y)))
sZ, sX = zip(*sorted(zip(Z,X)))
# fit the PAVA
pavfit = pav(sY).tolist()
# get the error of t'th point
x, z, y = adv.train_data[0][t], adv.train_data[1][t], adv.train_data[2][t]
# search z in sZ and get the adjacent point
pos = bisect.bisect_left(sZ, z)
if pos == 0:
err = math.pow(pavfit[0] - y, 2)
E = E + err
print(t, err, E)
elif pos == len(sZ):
err = math.pow(sY[pavfit(sY)-1] - y, 2)
E = E + err
print(t, err, E)
else:
z1, z2, y1, y2 = sZ[pos-1], sZ[pos], pavfit[pos-1], pavfit[pos]
err = math.pow(y1 + (y2-y1)*(z-z1)/(z2-z1) - y, 2)
E = E + err
print(t, err, E)
# import pylab as pl
# pl.close('all')
# pl.plot(sZ, sY, 'rx')
# pl.plot(sZ, pavfit, 'b')
# pl.show()