Python exemplarisch - Data Mining
deutsch     english

Regressionsanalyse

 

Alle Programme können von hier heruntergeladen werden.

 

Problembeschreibung
 

Bei der Regressionsanalyse studiert man, wie eine bestimmte Grösse (ein Ausgabewert) mit anderen variablen Grössen (Eingabewerten) zusammenhängt, beispielsweise der Wiederverkaufswert z eines bestimmten Autotyps vom Alter x und Kilometerstand y. Dabei versucht man, den Ausgabewert z in eine eindeutige Beziehung zu den Eingabegrössen x, y zu bringen, mathematisch ausgedrückt, eine Funktion z = f(x, y) mit zwei Variablen x, y zu finden, welche den Wiederverkaufswert "am besten" wiedergibt.


 

Eindimensionale lineare Regression
 

m einfachsten Fall hat die Funktion nur eine einzige Variablen x und die Beziehung eine lineare Funktion y = f(x) = a * x + b mit den zwei Parametern a und b. Die Funktion kann als Gerade in einem x-y-Koordinatensystem dargestellt werden.

Ausgehend von einem Trainigs-Set stellen wir uns die Aufgabe, die Parameter optimal zu bestimmen. Was versteht man darunter?

Wir betrachten ein Beispiel: Die Trainigs-Daten sind als sechs [x, y]-Wertepaare gegeben und werden in einer Liste X zusammengefasst.

X = [[1, 2], [3, 3], [4, 5], [5, 6], [6, 5], [8, 7]]

Wir können sie als Punkte in einem GPanel einzeichnen, wobei wir zur besseren Sichtbarkeit kleine Kreise zeichnen.

Programm: [►]

# Regression1.py

from gpanel import *

samples = [[1, 2], [3, 3], [4, 5], [5, 6], [6, 5], [8, 7]]

def drawData():
    for sample in samples:
        pos(sample)
        fillCircle(0.07)
 
makeGPanel(-1, 11, -1, 11, mousePressed = onMousePressed)
drawGrid(0, 10, 0, 10, "gray")
title("Data Points")
drawData()
Programmcode markieren (Ctrl+C kopieren, Ctrl+V einfügen)

regression1

Bemerkungen:
Die Punkte wurden absichtlich ziemlich "verstreut" ausgewählt.

 

Optimale Funktion -> Kleinste Fehlersumme
 

Denken wir, dass die beiden Parameter a (Steigung) und b (Offset) einen bestimmten Wert haben: a = 0.5, b = 3 und tragen wir die Gerade z = f(x) = 0.5 * x + 3 im Koordinatensystem ein. Diese Gerade ist offenbar eine schlechte Näherung für die gegebenen Datenpunkte, da sie weit daneben verläuft. Zu jedem Datenpunkt können wir einen Fehlerwert e als vertikalen Abstand des Datenpunkts von der Gerade angeben: e = abs(f(x) - y).

Programm: [►]

# Regression2.py

from gpanel import *

samples = [[1, 2], [3, 3], [4, 5], [5, 6], [6, 5], [8, 7]]
a = 0.5
b = 3

def f(x):
    return a * x + b

def drawData():
    for sample in samples:
        pos(sample)
        fillCircle(0.07)
    
def drawLine():
    line(0, f(0), 10, f(10))
    
def showErrors():
    errorSum = 0
    setColor("red")
    for sample in samples:
        x = sample[0]
        y = sample[1]
        e = abs(y - f(x))
        errorSum += e
        line(x, y, x, f(x)) 
        text(x + 0.1, y, str(e))
    return errorSum  

makeGPanel(-1, 11, -1, 11, mousePressed = onMousePressed)
drawGrid(0, 10, 0, 10, "gray")
title("Data Points")
lineWidth(2)
drawData()
drawLine()
E = showErrors()
title("Error sum: " + str(E))
Programmcode markieren (Ctrl+C kopieren, Ctrl+V einfügen)

regression2

Bemerkungen:
Durch Veränderung der Werte von a und b kann man versuchen, die Fehlersumme klein zu machen. Die optimale Gerade mit der kleinsten Fehlersumme kann man aber durch dieses Verfahren kaum finden.

 

Kleinste Summe der Quadrate
  Um die Gerade mit der kleinesten Fehlersumme mit einem Algorithmus oder. einem mathematischen Verfahren zu finden, ist es einfacher, statt der vertikalen Abstände
regression3
die Quadrate
 regression4
zu verwenden.

Die Aufgabe besteht dann darin, dann die "Kleinste Summe der Quadrate" zu finden, also im Fall einer Geraden y = a*x + b, in der Fehlersumme
                                             regression5
die Grössen a und b so zu bestimmen, dass E einen minimalen Wert hat.

Für eine gegebene Datensammlung von Wertepaaren (x, y)] kann man E als eine Funktion der Variablen a und b auffassen und grafisch darstellen. Es handelt sich um eine Fläche im 3D-Raum. Für die obigen Werte ergibt sich folgende Fläche:

regression6

Mit der Python-Bibliothek matplotlib kann man sie mit wenigen Zeilen Code selbst erstellen:

Programm: [►]

# RegressionPlot.py

from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import matplotlib.pyplot as plt
import numpy as np

fig = plt.figure()
ax = fig.gca(projection='3d')

a = [i for i in range(-1000,1000, 200)]
b = [i for i in range(-10000, 10000, 2000)]
a, b = np.meshgrid(a, b)
E = np.array((a * 1 + b - 2)**2 + (a * 3 + b - 3)**2 + (a * 4 + b - 5)**2
           + (a * 5 + b - 6)**2 + (a * 6 + b - 5)**2 + (a * 8 + b - 7)**2)
ax.plot_surface(a, b, E, rstride=1, cstride=1, cmap=cm.autumn, linewidth=0)

ax.set_xlabel('a')
ax.set_ylabel('b')
ax.set_zlabel('E')
plt.axis((-1000, 1000, -10000, 10000))
plt.show()
Programmcode markieren (Ctrl+C kopieren, Ctrl+V einfügen)

Um das Minimum für a und b zu finden, muss man sich auf der Falllinie (descent) nach unten bewegen, also entgegen der Richtung der grössten Steigung. Aus der Mathematik weiss man, dass man die Komponenten des Steigungsvektor (den Gradienten) so erhält, dass man die Funktion nach den beiden Variablen ableitet. Die beiden Ableitungen lauten:

  regression7 regression8

Um das Minimum zu bestimmen gehen wir im Folgenden nicht algebraisch vor, sondern verwenden einen bekannten Algorithmus, den wir auch auf komplexere Funktionen anwenden können. Dazu steigen wir in kleinen Schritten iterativ die Fläche in Richtung des steilsten Abfalls (steepest descent) hinab, bis sich der Funktionswert nicht mehr wesentlich ändert.( Das Verfahren heisst "gradient descent". Die Richtung des steilsten Abstiegs entspricht dem negativen Gradienten.)

Mit jedem Mausklick wird der nächste Iterationsschritt ausgeführt. Wir zeichnen zum Vergleich auch noch die Regressionsgerade ein, die man aus einer Library-Funktion erhält.

Programm: [►]

# Regression3.py

from gpanel import *
import time

# Set of samples
samples = [[1, 2], [3, 3], [4, 5], [5, 6], [6, 5], [8, 7]]
# Iteration step in both directions
ds = 0.001

def f(x):
    return a * x + b

def showLine():
    global y0, y10
    line(0, y0, 10, y10)  # erase old line
    # New values        
    y0 = f(0)  
    y10 = f(10)    
    line(0, y0, 10, y10) # draw new line

def gradient_a():
    # Derivative with respect to a (slope)
    za = 0 
    for sample in samples:
        xi = sample[0]
        yi = sample[1]
        za += (f(xi) - yi) * xi  
    return 2 * za

def gradient_b():
    # Derivative with respect to b (offset)
    zb = 0
    for sample in samples:
        xi = sample[0]
        yi = sample[1]
        zb += (f(xi) - yi)
    return 2 * zb

def onMousePressed(x, y):
    global iterationNb
    global a, b
    for i in range(100):  # 100 iterations per click
        iterationNb += 1
        anew = a - ds * gradient_a() # descent in a direction   
        bnew = b - ds * gradient_b() # descent in b direction
        # updata a, b
        a = anew
        b = bnew
    showLine()
    showResult()

def showResult():
    title("Iteration # " + str(iterationNb) + ": Slope = "  + \
      str(round(a, 3)) + ", Offset = " + str(round(b, 3))) 
 
makeGPanel(-1, 11, -1, 11, mousePressed = onMousePressed)
drawGrid(0, 10, 0, 10, "gray")
iterationNb = 0
title("Iteration # " + str(iterationNb))
for sample in samples:
    pos(sample)
    fillCircle(0.07)

# Linear fit using library function
xdata = [samples[i][0] for i in range(len(samples))]
ydata = [samples[i][1] for i in range(len(samples))]
a, b = linfit(xdata, ydata)  # linear fit using library function
text(0, 0.7, "Slope = " + str(round(a, 3)))
text(0, 0.3, "Offset = " + str(round(b, 3)))
setColor("green")
lineWidth(3)
line(0, f(0), 10, f(10))

# Prepare iteration
a = 0.5 # first guess for slope
b = 3   # first guess for offset
lineWidth(1)
setColor("red")
setXORMode("white")
y0 = f(0)  
y10 = f(10)    
line(0, y0, 10, y10)
showResult()
Programmcode markieren (Ctrl+C kopieren, Ctrl+V einfügen)

regression9

Bemerkungen:
Die Anfangswerte a, b dürfen nicht zu weit vom Minimum gewählt werden. Wenn die Funktion mehrere "Täler" hat, ist man nicht sicher, das global tiefte Tal zu finden.

 

Multidimensionale lineare Regression
 

Wir packen nun das eingangs erwähnte Problem an, eine Voraussage über den Wiederverkaufswert von Autos zu machen. Davon gehen wir von folgendem Trainings-Set aus.

Alter
(in Jahren)		 KM-Stand
(in 1000 km)		 Verkaufspreis
(in 1000 SFr)
5 57 8.5
4 40 10.3
6 77 7
5 60 8,2
5 49 8.9
5 47 9.8
6 58 6.6
6 39 9.5
2 8 16.9
7 69 7
4 89 4.8

Wir machen einen linearen Modellansatz: z = f(x, y) = a0 + a1* x + a2 * y, wo x das Alter, y der KM-Stand und z der Verkaufspreis sind. Die Parameter a0, a1, a2 müssen mit dem Regressionsverfahren optimal angepasst werden (a0 entspricht dem Neuwert, a1 und a2 sind sicher negativ).

Die Fehlersumme und die Ableitungen nach den drei Parametern lauten:

  regression10
  regression11 regression12 regression13

Wir gehen wie oben schrittweise vor und schreiben immer nach 10'000 Iterationen die Parameter und die Fehlersumme aus. Nach 50'000 Iteration ist die Konvergenz gut genug und wir können den voraussichtlichen Occasionspreis eines siebenjährigen Autos mit einem Kilometerstand von 60'000 km bestimmen.

Programm: [►]

# Regression4.py

import time

# Set of samples [x1, x2, y]
# Used car price [years, km /1000, price SFr/1000]
samples = [[5, 57, 8.5], [4, 40, 10.3], [6, 77, 7], [5, 60, 8.2], 
           [5, 49, 8.9], [5, 47, 9.8], [6, 58, 6.6], [6, 39, 9.5], 
           [2, 8, 16.9], [7, 69, 7], [7, 89, 4.8]]
# Iteration step in both directions
ds = 0.00001

def f(x1, x2):
    return a1*x1 + a2*x2 + a0

def error():
    e = 0
    for sample in samples:
        de = (f(sample[0], sample[1]) - sample[2])**2 
        e += de
    return e
    
def gradient_a1():
    # Derivative with respect to a1
    s = 0 
    for sample in samples:
        x1i = sample[0]
        x2i = sample[1]
        yi = sample[2]
        s += (f(x1i, x2i) - yi) * x1i  
    return 2 * s

def gradient_a2():
    # Derivative with respect to a2
    s = 0 
    for sample in samples:
        x1i = sample[0]
        x2i = sample[1]
        yi = sample[2]
        s += (f(x1i, x2i) - yi) * x2i  
    return 2 * s

def gradient_a0():
    # Derivative with respect to a0
    s = 0
    for sample in samples:
        x1i = sample[0]
        x2i = sample[1]
        yi = sample[2]
        s += (f(x1i, x2i) - yi)  
    return 2 * s

# Start values
a1 = 0
a2 = 0
a0 = 0
iterationNb = 0
print "Learning now. Please wait..."
for i in range(500000): 
    iterationNb += 1
    a1new = a1 - ds * gradient_a1() # descent in a1 direction   
    a2new = a2 - ds * gradient_a2() # descent in a2 direction
    a0new = a0 - ds * gradient_a0() # descent in a0 direction   
    # updata a1, a2, a0
    a1 = a1new
    a2 = a2new
    a0 = a0new
    if iterationNb % 100000 == 0:
        print "%3d iters-> a1 a2 a0: %4.2f %4.2f %4.2f -- Error: %4.2f " \
               %(iterationNb, a1, a2, a0, error())
print "Learning process terminated"
print "Prediction for a 7-years old car with 60'000 km:", \
       int(1000 * f(7, 60)), "SFr"
Programmcode markieren (Ctrl+C kopieren, Ctrl+V einfügen)

Resultat:
Learning now. Please wait...
100000 iters-> a1 a2 a0: 0.57 -0.13 12.69 -- Error: 25.22
200000 iters-> a1 a2 a0: -0.46 -0.10 16.48 -- Error: 8.20
300000 iters-> a1 a2 a0: -0.79 -0.09 17.71 -- Error: 6.41
400000 iters-> a1 a2 a0: -0.90 -0.08 18.11 -- Error: 6.22
500000 iters-> a1 a2 a0: -0.93 -0.08 18.24 -- Error: 6.20
Learning process terminated
Prediction for a 7-years old car with 60'000 km: 6744 SFr

Bemerkungen:
Das Verfahren lässt sich sogar auf nicht-lineare Regressionen mit beliebig vielen Parametern anwenden.