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
'''
Least absolute deviations (linear case)

Sylvain, 2010.

Reference: http://en.wikipedia.org/wiki/Least_absolute_deviations
'''
#-------------------------------------------------------------------------------
def lad(X, Y):
    '''
    Trivial (brute force) LAD implementation for the linear case. 

    Given a data set consisting of the points (x_i, y_i) for i in [1, n], We
    want to find a function f(x) = a * x + b minimizing S, where: 

    S = sum(abs(y_i - f(x_i))), for i = 1, 2, ..., n.

    It returns a 5-tuple: (l, m, a, b, S) where (x_l, y_l) and (x_m, y_m) are
    the two perfect fits of the data set.

    WARNING: This is O(n^2).
    '''
    def comb_n_2(n):
        for i in xrange(n):
            for j in xrange(i+1, n):
                yield i, j
    def fit_all(X, Y):
        for i, j in comb_n_2(len(X)):
            a = (Y[j] - Y[i]) / (X[j] - X[i])
            b = Y[j] - a * X[j]
            yield i, j, a, b, sum(abs(Y - (a * X + b)))
    return min(fit_all(X, Y), key = lambda i:i[4])

#-------------------------------------------------------------------------------
# Exemplification
if __name__ == '__main__':
    # Interpolation
    from numpy import array
    X_737 = array([ -4.30, -15.90, -10.70,  2.10])
    Y_737 = array([  1.92,   3.47,   3.47,  3.05])
    X_739 = array([-15.30, -26.90, -21.70, -8.90])
    Y_739 = array([  2.81,   3.85,   3.40,  1.99])

    a1, b1   = lad(X_737, Y_737)[2:4]
    YFIT_737 = a1 * X_737 + b1

    a2, b2   = lad(X_739, Y_739)[2:4]
    YFIT_739 = a2 * X_739 + b2

    # Calcul de la consommation
    from numpy import frompyfunc
    positivate = frompyfunc(lambda x: max(x, 0), 1, 1)
    c = sum(positivate(a1 * X_737)) * 1000

    # Production du graphe
    from matplotlib import pyplot
    f = pyplot.figure()
    p = f.add_subplot(111)

    p.set_title (u'Consommation \xe9lectrique, hiver 2009-2010')
    p.set_ylabel(u'Consommation (MWh)')
    p.set_xlabel(u'Gradient de temp\xe9rature moyen t (\xb0C)')
    p.plot(X_737, Y_737   , 'ro')
    p.plot(X_737, YFIT_737, 'r', label = '737 Avenue Outremont')
    p.plot(X_739, Y_739   , 'bo')
    p.plot(X_739, YFIT_739, 'b', label = '739 Avenue Outremont') 

    p.text( -8, 3.3, u'%.3ft + %.3f' % (a1, b1))
    p.text(-27,   3, u'%.3ft + %.3f' % (a2, b2))
    p.text(-14,  .1, u'Chauffage estim\xe9 du garage: %d kWh' % c)
    p.set_ylim(0, 5)
    p.legend()
    pyplot.show()

#-------------------------------------------------------------------------------
# That's all, folks!