DIY anomaly detection

Een van de interessante toepassingsgebieden van machine learning is, denk ik, anomaly detection: Vaststellen of iets (sterk) afwijkt van 'het gewone'. Bijvoorbeeld als een bepaald vergelijkbaar patroon 'ineens' af gaat wijken. Vooral re-current netwerken lijken mij daarvoor geschikt omdat ze zaken 'in te tijd' kunnen bekijken. Of te wel bijvoorbeeld een opeenvolging van vaste gebeurtenissen.

Denk hierbij bijvoorbeeld aan het herkennen van dronkenmans gedrag tussen 'normale' wachtende op de trein.

Het schiet mij te binnen dat een 'normaal' neuraal netwerk ook een soort 'anomaly detection' zou kunnen doen. Het mnist cijferherkennings convolutional netwerk (wow wat een zin! ), uit de eerste oefeningen met NN, zou een betrouwbaarheids indicatie mee kunnen geven. Is dit bij een resultaat erg laag dan kan dat wellicht een 'anomalie' aangegeven.

Even in de herinnering: 60.000 handgeschreven cijfers van 0 t/m 9 kunnen worden 'aangeleerd'. De uitkomst staat in een 10 cijfers grote matrix waarvan de index correspondeert met het verwachte getal. We kregen daarbij een behoorlijk hoge betrouwbaarheid van 95,5%.  

Voorbeeld:

Deze uitkomst geeft aan dat het gelezen getal waarschijnlijk een 6 is (Op de zesde positie staat het hoogste getal : 0.18488647):

[ 0.13371149  0.06534619  0.09990948  0.0668817   0.08472356  0.13552423   0.18488647  0.04352311  0.11824396  0.06724983]

Na wat experimenteren lukt het mij om deze getallen afgedrukt te krijgen. Maar hoe bepaal je nou de betrouwbaarheid? Als ik probeer het hoogste getal te delen door de som van de getallen blijkt deze som steevast op 1 uit te komen. Ook na 100 minibatches (leerloops).  Mmm ... Dat had ik eerder nog niet echt door. Ik snap ook nog niet waar dat in het programma precies wordt gewaarborgd. 

Maar goed. Maar eens kijken of een lage hoogste waarde van deze 10 inderdaad duidelijk de anomaly cijfers laat zien. Als ik de waardes onder de 0.15 selecteer zijn de 'anomaly' cijfers wel wat bijzonder maar niet spectaculair. Eens even kijken wat bijvoorbeeld 0.12 doet. (Het theoretisch minimum zit op 0.1)

0.12 crashed. Dat wil waarschijnlijk zeggen dat er geen of iig minder dan 10 cijfers met een lagere score dan 0.12 zijn te vinden.   Bij 0.13 lijkt hij wel resultaten te geven Nog steeds zijn de gevonden cijfers niet spectaculair afwijkend:

Okeee, oeps ... dat blijkt toch een programmeerfoutje. Ook boven de .30 geeft hij dezelfde cijfers terug. Dit is wel het echte resultaat bij < 0.13: 

Lijkt inderdaad veel beter. Aardig voorstelbaar als 'anomaly'. Met stift geschreven? 

Dit zijn bijvoorbeeld cijfers waar 'hij' (begint toch een beetje levend te worden) wel behoorlijk zeker van is (>0.50) :

Tja toch weer wat verrassend. Vooral de 6 hier. Maar dat komt waarschijnlijk omdat wij ons eigen beeld hebben van wat een 'juiste schrijfwijze' is. Bijzonder hoe een klein streepje bij de nul hem heel veel meer onzeker kan maken.  

Met dit algoritme zou je iig dus op deze wijze bijvoorbeeld bij de brievenpost alleen de 'afwijkende cijfers' nog eens extra door een mens kunnen laten beoordelen. Aardig eerste 'anomaly detection' voorbeeld. Kijken wat we later nog meer kunnen verzinnen. De afwijking tussen het beste getal en het 2e is misschien een betere indicatie voor de onzekerheid. Voor een andere keer. 

Oh ja, hier weer de code. 

#export PYTHONPATH="${PYTHONPATH}:/Users/DWW"

#execfile('/Users/DWW/Documents/5_con_anomally.py')

import theano

from theano import tensor as T

from theano.sandbox.rng_mrg import MRG_RandomStreams as RandomStreams

import numpy as np

from load import mnist

from theano.tensor.nnet.conv import conv2d

from theano.tensor.signal.downsample import max_pool_2d

import matplotlib.pyplot as plt

from mpl_toolkits.axes_grid1 import ImageGrid

import time

srng = RandomStreams()

def floatX(X):

    return np.asarray(X, dtype=theano.config.floatX)

def init_weights(shape):

    return theano.shared(floatX(np.random.randn(*shape) * 0.01))

def rectify(X):

    return T.maximum(X, 0.)

def softmax(X):

    e_x = T.exp(X - X.max(axis=1).dimshuffle(0, 'x'))

    return e_x / e_x.sum(axis=1).dimshuffle(0, 'x')

def dropout(X, p=0.):

    if p > 0:

        retain_prob = 1 - p

        X *= srng.binomial(X.shape, p=retain_prob, dtype=theano.config.floatX)

        X /= retain_prob

    return X

def RMSprop(cost, params, lr=0.001, rho=0.9, epsilon=1e-6):

    grads = T.grad(cost=cost, wrt=params)

    updates = []

    for p, g in zip(params, grads):

        acc = theano.shared(p.get_value() * 0.)

        acc_new = rho * acc + (1 - rho) * g ** 2

        gradient_scaling = T.sqrt(acc_new + epsilon)

        g = g / gradient_scaling

        updates.append((acc, acc_new))

        updates.append((p, p - lr * g))

    return updates

def model(X, w, w2, w3, w4, p_drop_conv, p_drop_hidden, p_drop):

    X = dropout(X, p_drop)

    l1a = rectify(conv2d(X, w, border_mode='full'))

    l1 = max_pool_2d(l1a, (2, 2))

    l1 = dropout(l1, p_drop_conv)

    l2a = rectify(conv2d(l1, w2))

    l2 = max_pool_2d(l2a, (2, 2))

    l2 = dropout(l2, p_drop_conv)

    l3a = rectify(conv2d(l2, w3))

    l3b = max_pool_2d(l3a, (2, 2))

    l3 = T.flatten(l3b, outdim=2)

    l3 = dropout(l3, p_drop_conv)

    l4 = rectify(T.dot(l3, w4))

    l4 = dropout(l4, p_drop_hidden)

    pyx = softmax(T.dot(l4, w_o))

    return l1, l2, l3, l4, pyx

trX, teX, trY, teY = mnist(onehot=True)

trX = trX.reshape(-1, 1, 28, 28)

teX = teX.reshape(-1, 1, 28, 28)

X = T.ftensor4()

Y = T.fmatrix()

w = init_weights((32, 1, 3, 3))

w2 = init_weights((64, 32, 3, 3))

w3 = init_weights((128, 64, 3, 3))

w4 = init_weights((128 * 3 * 3, 625))

w_o = init_weights((625, 10))

noise_l1, noise_l2, noise_l3, noise_l4, noise_py_x = model(X, w, w2, w3, w4, 0.2, 0.5, 0.8)

l1, l2, l3, l4, py_x = model(X, w, w2, w3, w4, 0., 0., 0.)

y_x = T.argmax(py_x, axis=1)

#ch = py_x #[T.argmax(py_x, axis=1)]

ch = T.max(py_x, axis=1)

ser = py_x

cost = T.mean(T.nnet.categorical_crossentropy(noise_py_x, Y))

params = [w, w2, w3, w4, w_o]

updates = RMSprop(cost, params, lr=0.001)

train = theano.function(inputs=[X, Y], outputs=cost, updates=updates, allow_input_downcast=True)

predict = theano.function(inputs=[X], outputs=y_x, allow_input_downcast=True)

# Hier de uitbreinding van het model om de betrouwbaarheid bij de uitkomst te krijgen.

change = theano.function(inputs=[X], outputs=ch, allow_input_downcast=True) # Hoogste waarde uit output matrix

serie = theano.function(inputs=[X], outputs=ser, allow_input_downcast=True) # Alle waardes uit output matrix

t0 = time.time()

for i in range(20):

    for start, end in zip(range(0, len(trX), 128), range(128, len(trX), 128)):

        cost = train(trX[start:end], trY[start:end])

    print i, time.time() - t0, "seconds", np.mean(np.argmax(teY, axis=1) == predict(teX))

    t0 = time.time()

#Print enkele voorbeelden

#for i in range(100):

#    print '>>',np.argmax(teY[i]), predict(teX[i:i+1]), change(teX[i:i+1]), serie(teX[i:i+1])

def fillim(c):

    im = teX[c]

    im.shape = 28,28

    return im

fig = plt.figure(1, (10., 10.))

grid = ImageGrid(fig, 111, # similar to subplot(111)

                 nrows_ncols = (2, 5), # creates 2x5 grid of axes

                 axes_pad=0.1, # pad between axes in inch.

                 )

c = 0

i = 0

while c < 10:

    i += 1

    x =change(teX[i:i+1])

    if x[0] < 0.13: # Dit bepaald selectie op de betrouwbaarheid van de betreffende voorspelling

        grid[c].imshow(fillim(i),cmap=plt.cm.gray)

        c += 1

plt.show()