Regression

Abschnitt 3

Regression mit Python

0 von 7 Übungen erledigt

Sobald es um grössere Datenmengen oder um schöne Grafiken geht, sind Taschenrechner und Tabellenkalkulation schnell am Ende. In der Praxis werden Regressionen deshalb meist programmiert. Wir benutzen dazu Python mit zwei Bibliotheken:

Die wichtigsten Befehle

np.array([2, 3, 5])
erzeugt ein Array, also eine Liste von Zahlen, mit der man rechnen kann. Es gilt zum Beispiel 2 * np.array([2, 3, 5]) = np.array([4, 6, 10]).
np.polyfit(x, y, 1)
berechnet die Regressionsgerade zu den Datenpunkten und liefert die Liste [m, b] mit der Steigung \(m\) und dem y-Achsenabschnitt \(b\) zurück. Die 1 steht für den Grad des Polynoms: Mit einer 2 erhalten wir eine quadratische Regression (also die Liste [a, b, c] zu \(y = a \cdot x^2 + b \cdot x + c\)).
np.corrcoef(x, y)[0, 1]
liefert den empirischen Korrelationskoeffizienten \(r\) (siehe Abschnitt 4).
np.linspace(0, 5, 100)
liefert 100 gleichmässig verteilte Werte zwischen 0 und 5. Damit zeichnen wir Kurven.
plt.scatter(x, y)
zeichnet ein Streudiagramm, plt.plot(x, y) verbindet die Punkte zu einer Linie.

Lineare Regression

Wir bestimmen die Regressionsgerade zu den Daten aus Abschnitt 1 (Körpergrösse gegen Körpermasse) noch einmal, diesmal mit Python:

Beispiel. Lineare Regression mit numpy

Wir erhalten also genau dieselbe Regressionsgerade \(y \approx 0.9095 \cdot x - 89.5941\) und denselben empirischen Korrelationskoeffizienten \(r \approx 0.9035\) wie zuvor mit dem Taschenrechner. Das ist wenig überraschend: Der Rechner und Python lösen dasselbe Minimierungsproblem, nur eben mit unterschiedlichen Knöpfen.

Nicht-lineare Regression durch Linearisierung

Für eine exponentielle Regression suchen wir eine Funktion der Form \(h(t) = a \cdot e^{k \cdot t}\). Logarithmieren wir diese Gleichung, so erhalten wir

\[\ln(h(t)) = k \cdot t + \ln(a),\]

also einen linearen Zusammenhang zwischen t und ln(h). Damit lässt sich das Problem auf eine lineare Regression zurückführen: Wir führen die lineare Regression nicht mit den Daten \((t_i, h_i)\), sondern mit den Daten \((t_i, \ln(h_i))\) durch. Genauso rechnet auch der Taschenrechner.

Beispiel. Exponentielle Regression durch Linearisierung — die Schaumkronen-Daten aus Abschnitt 1

Derselbe Trick für eine Potenzregression

Für \(y = a \cdot x^n\) (Potenzregression) gilt \(\ln(y) = n \cdot \ln(x) + \ln(a)\). Wir führen die lineare Regression also mit den Daten \((\ln(x_i), \ln(y_i))\) durch, in Python mit np.polyfit(np.log(x), np.log(y), 1).

Beliebige Modelle

Nicht jedes Modell lässt sich linearisieren (zum Beispiel eine sinusförmige oder eine logistische Funktion). Für solche Fälle gibt es in der Bibliothek scipy die Funktion curve_fit: Wir geben die Gestalt der gesuchten Funktion als Python-Funktion vor, und curve_fit bestimmt diejenigen Parameter, für welche die Summe der Residuenquadrate minimal wird.

Beispiel. Logistische Regression mit curve_fit

Die Grösse \(G \approx 10\,880.01\) ist hier die sogenannte Kapazitätsgrenze: Mehr als etwa 10 880 Kranke sind in diesem Modell nie möglich — vergleichen Sie mit Übung 1.07 weiter unten.

Ohne Startwerte geht es meist nicht

Anders als np.polyfit sucht curve_fit die Parameter iterativ, ausgehend von den Startwerten p0. Ohne sinnvolle Startwerte findet curve_fit die Parameter oft gar nicht (oder ein falsches Minimum). Eine grobe Schätzung «von Auge» reicht meistens.

Übungen

Für die Wetterstation Basel/Binningen sind die folgenden monatlichen Temperaturdurchschnitte bekannt:

Monat 1 2 3 4 5 6 7 8 9 10 11 12
Temperatur (°C) 1.6 2.7 6.6 10.0 14.2 17.4 19.7 19.1 15.1 10.9 5.5 2.6
  1. Stellen Sie die Daten in einem Streudiagramm dar und passen Sie mit curve_fit eine Funktion der Form \(f(x) = a \cdot \sin(b \cdot x + c) + d\) an die Daten an. Zeichnen Sie die gefundene Kurve zusammen mit den Daten. (Hinweis: Wählen Sie als Startwerte p0=(9, 0.5, -2, 10). Ohne sinnvolle Startwerte findet curve_fit die Parameter nicht.)
  2. Unter der minimalen Wachstumszeit einer Kulturpflanze versteht man die kleinste Anzahl Tage im Jahr, an denen die Tagesmitteltemperatur über 10 °C liegen muss, damit die Pflanze erfolgreich angebaut werden kann. Ist die Region Basel für den Anbau einer Nutzpflanze mit einer minimalen Wachstumszeit von 190 Tagen geeignet? Werten Sie dazu Ihr Modell an vielen Tagen im Jahr aus und zählen Sie, an wie vielen Tagen die Temperatur über 10 °C liegt.
Lösungsvorschlag

(a) \(f(x) \approx 8.92 \cdot \sin(0.53 \cdot x - 2.16) + 10.54\)

(b) Die Daten lassen auf etwa 188 Tage im Jahr mit Tagesmitteltemperatur über 10 °C schliessen. Die Region ist für eine Pflanze mit einer minimalen Wachstumszeit von 190 Tagen also knapp nicht geeignet.

Lösungsvorschlag (wetter-sin.py)

Für die Ausbreitung einer Viruserkrankung wurde die Anzahl Kranker am Ende jeder Woche erfasst:

Woche 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
Kranke 2 4 7 17 32 60 129 254 502 979 1862 3365 5614 8076 9645 9986
  1. Wie lautet die Gleichung der entsprechenden Regressionskurve, wenn wir annehmen, dass die Ausbreitung logistisch erfolgt (siehe Beispiel oben)? Runden Sie die Parameter auf zwei Nachkommastellen.
  2. Gibt es nach diesem Modell jemals mehr als 11'000 Kranke?
Lösungsvorschlag

(a) \(y \approx 10\,880.01 / (1 + 70\,563.74 \cdot e^{-0.87 \cdot x})\)

(b) Nein: Die Anzahl Kranker kann in diesem Modell niemals \(G \approx 10\,880.01\) übersteigen, bleibt also immer unter 11'000.

Lösungsvorschlag (virus.py)

Ein Autohändler notiert für sechs gleichaltrige Occasionen des gleichen Modells den Kilometerstand x (in 1000 km) und den Verkaufspreis y (in 1000 CHF):

x (1000 km) 15 38 52 71 96 120
y (1000 CHF) 24.5 20.8 18.9 15.4 12.1 9.0
  1. Führen Sie eine lineare Regression durch und geben Sie die Gleichung der Regressionsgeraden sowie den empirischen Korrelationskoeffizienten \(r\) an.
  2. Was bedeuten die Steigung und der y-Achsenabschnitt in dieser Situation konkret?
  3. Welchen Preis erwarten Sie für ein Auto mit einem Kilometerstand von 60'000 km, welchen für ein Auto mit 150'000 km?
  4. Bei welchem Kilometerstand wäre das Auto nach diesem Modell gratis? Was sagt das über die Aussagekraft des Modells?
Lösungsvorschlag
  1. \(y \approx -0.1487 \cdot x + 26.4994\) mit \(r \approx -0.9985\)
  2. Der Preis sinkt um etwa 149 CHF pro 1000 km. Der y-Achsenabschnitt ist der Preis, den ein solches Auto mit Kilometerstand 0 hätte, also etwa 26'499 CHF.
  3. Etwa 17'576 CHF beziehungsweise etwa 4'192 CHF.
  4. Bei etwa 178'189 km. Das Modell darf also nicht beliebig weit extrapoliert werden: Ein reales Auto wird nicht gratis, und der Preis fällt erfahrungsgemäss nicht linear.

Lösungsvorschlag (auto.py)

Für ein Automodell wird der Bremsweg s (in m) bei verschiedenen Geschwindigkeiten v (in km/h) gemessen:

v (km/h) 30 50 70 90 110 130
s (m) 6 15 29 48 72 100
  1. Führen Sie eine lineare Regression durch. Wie gross ist \(r\)? Zeichnen Sie zusätzlich ein Streudiagramm. Was fällt auf?
  2. Führen Sie nun eine quadratische Regression durch. Geben Sie die Parameter auf fünf Nachkommastellen genau an.
  3. Welchen Bremsweg sagen die beiden Modelle für \(v = 160\) km/h voraus?
  4. Aus der Physik wissen wir, dass der Bremsweg proportional zum Quadrat der Geschwindigkeit ist. Welches der beiden Modelle ist also vorzuziehen? Vergleichen Sie dazu auch die Summe der Residuenquadrate beider Modelle.
Lösungsvorschlag
  1. \(s \approx 0.94286 \cdot v - 30.42857\) mit \(r \approx 0.98300\). Obwohl \(r\) nahe bei 1 liegt, sieht man im Streudiagramm deutlich, dass die Punkte auf einer nach oben gekrümmten Kurve und nicht auf einer Geraden liegen.
  2. \(s \approx 0.00603 \cdot v^2 - 0.02143 \cdot v + 1.11161\)
  3. Das lineare Modell liefert etwa 120.43 m, das quadratische etwa 151.97 m.
  4. Das quadratische Modell. Es passt nicht nur besser zu den Daten (Summe der Residuenquadrate \(\approx 0.1786\) statt \(\approx 217.1429\)), sondern auch zur Physik. Das lineare Modell unterschätzt den Bremsweg bei hohen Geschwindigkeiten deutlich.

Lösungsvorschlag (bremsquad.py)

Nach dem Trinken eines Espressos wird die Koffeinmenge c (in mg) im Blut stündlich gemessen:

t (h) 0 1 2 3 4 5 6
c (mg) 120 104 89 78 67 58 50
  1. Führen Sie eine exponentielle Regression durch und geben Sie die Funktionsgleichung an.
  2. Um wie viel Prozent nimmt die Koffeinmenge pro Stunde ab?
  3. Nach welcher Zeit ist nur noch die Hälfte des Koffeins im Blut (Halbwertszeit)?
  4. Wann sind weniger als 20 mg Koffein im Blut?
Lösungsvorschlag
  1. \(c(t) \approx 120.01 \cdot 0.86^t\)
  2. Um etwa 13.6 % pro Stunde.
  3. Nach etwa 4.76 h.
  4. Nach etwa 12.30 h, also nach gut einem halben Tag.

Lösungsvorschlag (koffein.py)

Für die sechs klassischen Planeten sind der mittlere Sonnenabstand a (in astronomischen Einheiten AE, wobei 1 AE der Abstand Erde–Sonne ist) und die Umlaufzeit T (in Jahren) gegeben:

Planet Merkur Venus Erde Mars Jupiter Saturn
a (AE) 0.387 0.723 1.000 1.524 5.203 9.537
T (Jahre) 0.241 0.615 1.000 1.881 11.862 29.457
  1. Führen Sie eine Potenzregression durch, also eine Regression mit einer Funktion der Form \(T = c \cdot a^n\) (siehe Hinweis zur Linearisierung weiter oben). Wie gross ist \(r\)?
  2. Der gefundene Exponent liegt sehr nahe bei einem Bruch mit kleinem Nenner. Bei welchem? Formulieren Sie damit einen Zusammenhang zwischen T und a ohne Kommazahlen. (Bemerkung: Sie haben damit das dritte Keplersche Gesetz aus Daten «wiederentdeckt».)
  3. Neptun ist 30.07 AE von der Sonne entfernt. Wie lange dauert dort ein Jahr?
Lösungsvorschlag
  1. \(T \approx 1.00 \cdot a^{1.4997}\) mit \(r \approx 1\) (genauer: \(r \approx 0.99999998\)).
  2. Der Exponent ist praktisch genau 3/2. Es gilt also \(T \approx a^{3/2}\) und damit \(T^2 \approx a^3\).
  3. Etwa 164.78 Jahre (der tatsächliche Wert ist 164.8 Jahre).

Lösungsvorschlag (kepler.py)

Eine Schülerin lernt jeden Tag Vokabeln und notiert, wie viele Wörter einer Liste sie nach x Tagen auswendig kennt:

x (Tage) 1 2 3 4 5 6 8 10
y (Wörter) 20 42 52 62 68 73 80 88
  1. Probieren Sie die vier Regressionstypen linear, logarithmisch (\(y = a + b \cdot \ln(x)\)), exponentiell und Potenz aus und notieren Sie jeweils \(r\) (nutzen Sie dazu np.corrcoef auf den passend linearisierten Daten). Welches Modell passt am besten zu den Daten?
  2. Geben Sie die Funktionsgleichung des besten Modells an.
  3. Wie viele Wörter kennt die Schülerin nach diesem Modell nach 20 Tagen? Nach wie vielen Tagen kennt sie 100 Wörter?
  4. Die Liste enthält 100 Wörter. Was sagt das über die Grenzen des Modells?
Lösungsvorschlag
  1. Man erhält \(r \approx 0.9395\) (linear), \(r \approx 0.9993\) (logarithmisch), \(r \approx 0.8466\) (exponentiell) und \(r \approx 0.9708\) (Potenz). Am besten passt die logarithmische Regression.
  2. \(y \approx 20.76 + 29.10 \cdot \ln(x)\)
  3. Etwa 108 Wörter nach 20 Tagen und 100 Wörter nach etwa 15.23 Tagen.
  4. Das Modell wächst unbeschränkt und sagt nach 20 Tagen mehr Wörter voraus, als die Liste überhaupt enthält. Ein Modell mit einer Sättigung (zum Beispiel ein logistisches) wäre hier realistischer. Ausserhalb des erhobenen Bereichs darf man einer Regression also nicht blind vertrauen.

Lösungsvorschlag (modellwahl.py)