Foute bussen - deel acht - Dataschoning

Bij het schonen van de data kom ik er eerst achter dat mijn tijdberekening niet goed werkt. Uiteraard mag je tijden niet altijd (11:06 - 10:58 ?!) gewoon van elkaar aftrekken om het verschil in minuten te krijgen. Wie heeft die 60 minuten per uur verzonnen? :-) Ok, kost even tijd om dat goed te fixen. Vervolgens kijken naar de max 100 km/uur om foute input te verwijderen. Deze fouten zijn vermoedelijk ontstaan omdat ik bij het opslaan van de data niet de datum heb meegenomen. Posities uit verschillende dagen kunnen dus bij elkaar zijn komen te staan. Lastig, maar een lering. Ik kan dus in ieder geval meten of de opeenvolgende punten niet meer dan 3,3 km uit elkaar liggen ervan uitgaande dat de chauffeurs 'normaal' deze snelheid niet halen met hun bussen.

Het levert prachtige resultaten op. Enerzijds juicht mijn statistiek:

count  1087.000000

mean      0.481143

std       0.286745

min       0.017168

25%       0.275199

50%       0.448470

75%       0.623757

max       2.952549  

Minder dan 0.5 km gemiddelde voorspellingsfout! Anderzijds blijken de 'anomalieën' beperkt. (maar 4 van de 1087 > 2.0 km afwijking) en lijken dat ook echt nog overgebleven fouten. Bijvoorbeeld:

Distance =  2.95254907837

Misschien is hier de buschauffeur razendsnel omgedraaid en met zo'n 87 (30 * 2.95...) km/uur naar huis geraced. Tja, het kan natuurlijk.

Maar ik krijg ook dit soort:

Distance =  2.08255048828

Een hele rare route die vermoedelijk dus met het eerdere probleem te maken heeft.

Wel, ik ben feitelijk waar ik wezen wilde. Het lijkt erop dat ik afwijkend gedrag van de bussen nu aardig kan waarnemen. Misschien nog een keertje het neuraal netwerk rechtstreeks koppelen aan de webservice. Of nu overstappen naar iets anders. We zien wel.

De code:

#execfile('/Users/DWW/Documents/Follow1bus outlier neural.py')

import webbrowser

import time

import numpy as np

from sklearn.externals import joblib

from sklearn import svm

import pandas as pd

Data = joblib.load('/Users/DWW/Documents/Bus_outlier.pkl')

Data = np.array(sorted(Data, key = lambda x: (x[0], x[1])))

busses = sorted(list(set([Colm[0] for Colm in Data])))

def dist(est,real): # Bereken afstand tussen werkelijke en voorspelde positie. Ook procentueel.

    if real == 0:

        perc = 0

    else:

        perc = np.absolute(100* (real-est)/real)

    dis = np.absolute(real-est) * 1.852 * 60

    return est, real, perc, dis

def dista(l1, b1, l2, b2):

    return (np.absolute(float(l1)-float(l2)) + np.absolute(float(b1)-float(b2))) * 1.852 * 60

def timdif(int_n, int_0): # Maak numerieke tijd op basis van bijv. '07:34 PM' en bereken het verschil in tijd in minuten.

    if len(int_n) == 7:

        int_n = '0' + int_n

    if len(int_0) == 7:

        int_0 = '0' + int_0

    if int_n[6:8] == 'PM':     # bereken uren

        u_n = 12 + int(int_n[0:2])

    else:

        u_n = int(int_n[0:2])

    if int_0[6:8] == 'PM':

        u_0 = 12 + int(int_0[0:2])

    else:

        u_0 = int(int_0[0:2])

    if u_n == u_0:    # kijk of uren gelijk zijn of 1 verschillen

        dif = int(int_n[3:5]) - int(int_0[3:5])

    else:

        if u_n == u_0 + 1:

            dif = int(int_n[3:5]) + (60-int(int_0[3:5]))

        #print int(int_n[3:5]) , (60-int(int_0[3:5]))

        else:

            dif = (u_n-u_0) * 60 # niet gelemaal goed maar wel onderscheidend

    return dif

# Korrecties op gps posities:

def cl(lon):

    return 100 * (float(lon) - 41.)

def rl(lon):

    return float(lon)/100 + 41.

def cb(lat):

    return 100 * (float(lat) + 87.)

def rb(lat):

    return float(lat)/100 - 87.

# Opbouwen positie bestand (4 eerdere waarnemingen en de 5e werkelijke 'doelpositie'.  

X = []

Y = []

t = 0

for i in busses:

    HData = Data[Data[:,0]==i]

    if len(HData) > 5: # Alleen voor bussen waar meer dan 5 waarnemingen van zijn.

        for j in range(len(HData)-4):

            if (timdif(HData[j+4][1],HData[j][1]) < 11 and

                dista(HData[j+0][2],HData[j+0][3],HData[j+1][2],HData[j+1][3])< 3.3 and

                dista(HData[j+1][2],HData[j+1][3],HData[j+2][2],HData[j+2][3])< 3.3 and

                dista(HData[j+2][2],HData[j+2][3],HData[j+3][2],HData[j+3][3])< 3.3 and

                dista(HData[j+3][2],HData[j+3][3],HData[j+4][2],HData[j+4][3])< 3.3): # Alleen de waarnemeningen die inderdaad in tijd bij elkaar liggen. Feitelijk maximaal 8 minuten : 11 genomen voor marge. En max 100 km / uur = 3.3 per 2 min uit elkaar.

                if X == []:

                    X = [cl(HData[j][2]), cl(HData[j+1][2]), cl(HData[j+2][2]), cl(HData[j+3][2]), cb(HData[j][3]), cb(HData[j+1][3]), cb(HData[j+2][3]), cb(HData[j+3][3])]

                    Y = [cl(HData[j+4][2]),cb(HData[j+4][3])]

                else:

                    X = np.vstack((X,[cl(HData[j][2]), cl(HData[j+1][2]), cl(HData[j+2][2]), cl(HData[j+3][2]), cb(HData[j][3]), cb(HData[j+1][3]), cb(HData[j+2][3]), cb(HData[j+3][3])]))

                    Y = np.vstack((Y,[cl(HData[j+4][2]),cb(HData[j+4][3])]))

            else:

                #print t, 'overgeslagen = ', timdif(HData[j+4][1],HData[j][1]), HData[j+4][1],HData[j][1]

                #print t,'Distance overgeslagen : ',dista(HData[j+4][2],HData[j+4][3],HData[j][2],HData[j][3])

                t += 1

si = -.1 * len(X)

print 'Aantal records in dataset (X) =', len(X), 'Aantal in testset (si) =', si

from sklearn.utils import shuffle

X, Y = shuffle(X,Y)

X = X.astype(np.float32)

Y = Y.astype(np.float32)

# Voor het neuraal nw maar eens Lasagne proberen.

from lasagne import layers

from lasagne.updates import nesterov_momentum

from nolearn.lasagne import NeuralNet

net1 = NeuralNet(

    layers=[  # three layers: one hidden layer

        ('input', layers.InputLayer),

        ('hidden1', layers.DenseLayer),

        ('output', layers.DenseLayer),

        ],

    # layer parameters:

    input_shape=(None, 8),  # 8 input values per batch

    hidden1_num_units=100,  # number of units in hidden layer

    output_nonlinearity=None,  # output layer uses identity function

    output_num_units=2,  # 2 target values - lengte , breedte

                 

    # optimization method:

    update=nesterov_momentum,

    update_learning_rate=0.000001,

    update_momentum=0.9,

                 

    regression=True,  # flag to indicate we're dealing with regression problem

    max_epochs=1600,  # we want to train this many epochs

    verbose=1,

    )

net1.fit(X[1:si], Y[1:si])

def Maakkaart(Cl1, Cl2, Cl3 , Cl4, Cb1, Cb2, Cb3, Cb4, Wl, Wb, El, Eb):

    flab = '&markers=color:blue%7Clabel:A%7C'

    la1 = str(Cl1)

    lo1 = str(Cb1)

    la2 = str(Cl2)

    lo2 = str(Cb2)

    la3 = str(Cl3)

    lo3 = str(Cb3)

    la4 = str(Cl4)

    lo4 = str(Cb4)

    la5 = str(Wl)

    lo5 = str(Wb)

    la6 = str(El)

    lo6 = str(Eb)

    lab = '&markers=color:blue%7Clabel:1%7C' + la1 + ',' + lo1 + '&markers=color:blue%7Clabel:2%7C' + la2 + ',' + lo2 + '&markers=color:blue%7Clabel:3%7C' + la3 + ',' + lo3 + '&markers=color:blue%7Clabel:4%7C' + la4 + ',' + lo4 + '&markers=color:green%7Clabel:5%7C' + la5 + ',' + lo5  + '&markers=color:red%7Clabel:E%7C' + la6 + ',' + lo6

    #print lab

    import webbrowser

    webbrowser.open('https://maps.googleapis.com/maps/api/staticmap?center=&zoom=14&size=1200x600&maptype=roadmap' + lab,1,True)

# Kijk in testset (waarden > si) of er anomalieën in zitten. In dit geval afwijking > 2 km. 

Tdist = np.array([])

printfirst = True

kaarten = 0

for i in range(len(X[si:])-1):

    j = si+i

    pred = net1.predict(X[j:j+1])

    pred = pred[0].tolist()

    lpred, bpred = pred[0], pred[1]

    l_est, l_real, l_perc, l_dis = dist(rl(lpred), rl(Y[j][0]))

    b_est, b_real, b_perc, b_dis = dist(rb(bpred), rb(Y[j][1]))

    Tdist = np.append(Tdist, [l_dis + b_dis])

    if ((l_dis + b_dis > 2.0) and kaarten < 10):

        kaarten += 1

        print 'Distance = ' , l_dis + b_dis

        Maakkaart(rl(X[j,0]),rl(X[j,1]),rl(X[j,2]),rl(X[j,3]),rb(X[j,4]),rb(X[j,5]),
        rb(X[j,6]),rb(X[j,7]),rl(Y[j,0]),rb(Y[j,1]),rl(lpred),rb(bpred))

# statistische analyse met Pandas

Td = pd.DataFrame(Tdist)

print(Td.describe())