#!/usr/bin/env python3 # Predict_then_Optimize.py """ Kapitel Prognose: Die bessere Prognose trifft die schlechtere Entscheidung. Eine Baeckerei muss jeden Abend entscheiden, wie viel sie fuer den naechsten Tag ansetzt. Zu wenig kostet die Marge des entgangenen Verkaufs, zu viel kostet den Einkaufspreis der Retoure. Das ist das Newsvendor-Problem aus dem Kapitel Unsicherheit - nur dass die Nachfrage diesmal nicht aus einer Verteilung kommt, sondern PROGNOSTIZIERT werden muss: aus Wochentag, Temperatur und Aktionstagen. Damit zerfaellt die Aufgabe in zwei Schritte, und genau an der Naht entsteht der Fehler, um den es hier geht: PREDICT ein Modell schaetzt die Nachfrage OPTIMIZE daraus wird eine Bestellmenge Der Prognostiker optimiert seinen Modellfehler, meist den MSE. Der Planer traegt die Kosten. Beide messen etwas anderes - und die beiden Masse widersprechen einander. Das Programm zeigt: 1. Vier Verfahren, verglichen nach MSE UND nach Entscheidungskosten. Das Verfahren mit dem BESTEN MSE hat die HOECHSTEN Kosten. 2. Warum ein pauschaler Sicherheitszuschlag zu kurz greift - die Streuung der Nachfrage haengt selbst von den Merkmalen ab. 3. Eine Messfalle, in die der Autor dieses Programms zuerst selbst getappt ist: Bei kurzen Testzeitraeumen ist der MSE-Vergleich nicht stabil. Benoetigt: numpy, scipy, scikit-learn """ from __future__ import annotations import numpy as np from scipy.stats import norm from sklearn.linear_model import LinearRegression, QuantileRegressor VERKAUFSPREIS = 9.0 EINKAUFSPREIS = 3.0 KOSTEN_FEHLMENGE = VERKAUFSPREIS - EINKAUFSPREIS # entgangene Marge: 6 EUR KOSTEN_UEBERHANG = EINKAUFSPREIS # Retoure: 3 EUR KRITISCHES_VERHAELTNIS = KOSTEN_FEHLMENGE / (KOSTEN_FEHLMENGE + KOSTEN_UEBERHANG) TAGE = 5000 # Simulation, siehe Hinweis unten TRAINING = 1000 SAAT = 11 def erzeuge_daten(tage: int = TAGE, saat: int = SAAT): """Taegliche Nachfrage mit Wochentag, Temperatur und Aktionstagen. Die entscheidende Eigenschaft steckt in 'streuung': An Aktionstagen ist die Nachfrage nicht nur hoeher, sondern auch viel UNSICHERER. Solche heteroskedastischen Daten sind der Normalfall - und der Grund, warum ein pauschaler Sicherheitszuschlag nicht genuegt. """ rng = np.random.default_rng(saat) wochentag = np.arange(tage) % 7 temperatur = (12 + 10 * np.sin(2 * np.pi * np.arange(tage) / 365) + rng.normal(0, 3, tage)) aktion = (rng.random(tage) < 0.15).astype(float) merkmale = np.column_stack([np.eye(7)[wochentag][:, 1:], temperatur, aktion]) erwartung = (120 + np.eye(7)[wochentag] @ np.array([0, 10, 12, 14, 18, 35, -40]) + 1.8 * temperatur + 45 * aktion) streuung = 8 + 22 * aktion nachfrage = np.maximum(0.0, erwartung + rng.normal(0, 1, tage) * streuung) return merkmale, nachfrage, aktion def tageskosten(bestellung: np.ndarray, nachfrage: np.ndarray) -> float: """Die Zahl, auf die es ankommt - und die kein Prognosemass kennt.""" fehlmenge = np.maximum(0.0, nachfrage - bestellung) ueberhang = np.maximum(0.0, bestellung - nachfrage) return float((KOSTEN_FEHLMENGE * fehlmenge + KOSTEN_UEBERHANG * ueberhang).mean()) if __name__ == "__main__": merkmale, nachfrage, aktion = erzeuge_daten() lernen = slice(0, TRAINING) pruefen = slice(TRAINING, TAGE) print("=" * 84) print(" DIE BESSERE PROGNOSE TRIFFT DIE SCHLECHTERE ENTSCHEIDUNG") print("=" * 84) print(f"Verkaufspreis {VERKAUFSPREIS:.0f} EUR, Einkauf {EINKAUFSPREIS:.0f} EUR.") print(f"Fehlmenge kostet {KOSTEN_FEHLMENGE:.0f} EUR, Ueberhang " f"{KOSTEN_UEBERHANG:.0f} EUR je Stueck.") print(f"Kritisches Verhaeltnis: {KRITISCHES_VERHAELTNIS:.3f} - der Planer sollte " f"also das") print(f"{KRITISCHES_VERHAELTNIS:.1%}-Quantil der Nachfrage bestellen, nicht ihren " f"Erwartungswert.") print(f"\nTraining: Tag 1 bis {TRAINING}. Bewertung: die restlichen " f"{TAGE - TRAINING} Tage.\n") # --- Die vier Verfahren ---------------------------------------------- kleinste_quadrate = LinearRegression().fit(merkmale[lernen], nachfrage[lernen]) punktprognose = kleinste_quadrate.predict(merkmale[pruefen]) restfehler = nachfrage[lernen] - kleinste_quadrate.predict(merkmale[lernen]) pauschalzuschlag = norm.ppf(KRITISCHES_VERHAELTNIS) * restfehler.std() # Ein Zuschlag, der nicht aus der Normalverteilung kommt, sondern direkt # auf den Trainingsdaten die Kosten minimiert. kandidaten = np.linspace(-10.0, 30.0, 401) trainingsprognose = kleinste_quadrate.predict(merkmale[lernen]) kostenzuschlag = float(kandidaten[np.argmin( [tageskosten(trainingsprognose + z, nachfrage[lernen]) for z in kandidaten])]) # Und das Verfahren, das von vornherein das richtige Quantil schaetzt. quantilmodell = QuantileRegressor(quantile=KRITISCHES_VERHAELTNIS, alpha=0.0, solver="highs") quantilmodell.fit(merkmale[lernen], nachfrage[lernen]) quantilprognose = quantilmodell.predict(merkmale[pruefen]) verfahren = [ ("bestelle die Punktprognose", punktprognose, punktprognose), ("+ Zuschlag aus der Normalverteilung", punktprognose, punktprognose + pauschalzuschlag), ("+ Zuschlag auf Kosten trainiert", punktprognose, punktprognose + kostenzuschlag), ("Quantilregression aufs kritische Quantil", quantilprognose, quantilprognose), ] print(f" {'Verfahren':<42} {'MSE':>9} {'Kosten/Tag':>12} {'gegen Zeile 1':>14}") print(" " + "-" * 80) ergebnisse = {} for name, prognose, bestellung in verfahren: mse = float(((prognose - nachfrage[pruefen]) ** 2).mean()) kosten = tageskosten(bestellung, nachfrage[pruefen]) ergebnisse[name] = (mse, kosten) basis = ergebnisse[verfahren[0][0]][1] vergleich = "" if name == verfahren[0][0] else f"{(kosten - basis) / basis:+13.1%}" print(f" {name:<42} {mse:>9.1f} {kosten:>10.2f} EUR {vergleich:>14}") bester_mse = min(ergebnisse, key=lambda k: ergebnisse[k][0]) beste_kosten = min(ergebnisse, key=lambda k: ergebnisse[k][1]) print(f"\n bester MSE: {bester_mse}") print(f" beste Kosten: {beste_kosten}") print(f"\n Das Verfahren mit dem besten MSE hat die HOECHSTEN Kosten, und das") print(f" Verfahren mit den besten Kosten hat einen um " f"{(ergebnisse[beste_kosten][0] / ergebnisse[bester_mse][0] - 1):.0%} SCHLECHTEREN MSE.") print(" Wer Prognosemodelle nach MSE auswaehlt, waehlt hier das falsche.") # --- Warum der pauschale Zuschlag zu kurz greift --------------------- print("\n" + "-" * 84) print("Warum ein pauschaler Zuschlag nicht genuegt\n") ist_aktion = aktion[pruefen] > 0.5 print(f" {'Verfahren':<42} {'normale Tage':>14} {'Aktionstage':>14}") print(" " + "-" * 74) for name, _, bestellung in verfahren: normal = tageskosten(bestellung[~ist_aktion], nachfrage[pruefen][~ist_aktion]) aktionstag = tageskosten(bestellung[ist_aktion], nachfrage[pruefen][ist_aktion]) print(f" {name:<42} {normal:>10.2f} EUR {aktionstag:>10.2f} EUR") # Nachgerechnet statt behauptet: Welcher Zuschlag waere je Tagesart richtig? aktion_training = aktion[lernen] > 0.5 z = norm.ppf(KRITISCHES_VERHAELTNIS) richtig_normal = z * restfehler[~aktion_training].std() richtig_aktion = z * restfehler[aktion_training].std() print(f"\n Der pauschale Zuschlag betraegt {pauschalzuschlag:.1f} Stueck. Aus den") print(f" Trainingsresten getrennt nach Tagesart waere richtig:") print(f" normale Tage : {richtig_normal:5.1f} Stueck") print(f" Aktionstage : {richtig_aktion:5.1f} Stueck") print(f" Ein Zuschlag fuer alle Tage kann nur einen Mittelweg treffen - hier") print(f" ist er an normalen Tagen {pauschalzuschlag / richtig_normal:.1f}-mal zu gross und an") print(f" Aktionstagen nur {pauschalzuschlag / richtig_aktion:.0%} dessen, was noetig waere.") print(f"\n Die Quantilregression schaetzt das {KRITISCHES_VERHAELTNIS:.1%}-Quantil " f"direkt aus den") print(f" Merkmalen und darf deshalb an verschiedenen Tagen verschieden weit") print(f" ueber dem Erwartungswert liegen. Genau das ist der Unterschied") print(f" zwischen 'ein Modell und danach eine Formel' und 'ein Modell, das") print(f" weiss, wofuer es gebraucht wird'.") # --- Die Messfalle --------------------------------------------------- print("\n" + "-" * 84) print("Eine Messfalle, in die der Autor zuerst selbst getappt ist\n") print(" Der erste Entwurf dieses Programms bewertete auf 230 Testtagen - ein") print(" realistischer Zeitraum. Dort hatte die Quantilregression den BESSEREN") print(" MSE, und die ganze Aussage des Kapitels stand auf dem Kopf.") print("\n Wie oft das passiert, laesst sich ausmessen:\n") rng = np.random.default_rng(0) print(f" {'Testfenster':>14} {'QR sieht MSE-besser aus':>26}") print(" " + "-" * 42) for fenster in (180, 365, 730, 2000): treffer = 0 versuche = 400 for _ in range(versuche): start = int(rng.integers(TRAINING, TAGE - fenster)) ausschnitt = slice(start, start + fenster) mse_punkt = ((kleinste_quadrate.predict(merkmale[ausschnitt]) - nachfrage[ausschnitt]) ** 2).mean() mse_quantil = ((quantilmodell.predict(merkmale[ausschnitt]) - nachfrage[ausschnitt]) ** 2).mean() treffer += mse_quantil < mse_punkt print(f" {fenster:>10} Tage {treffer / versuche:>24.1%}") print("\n Bei einem halben Jahr Testdaten sieht das schlechtere Modell in gut") print(" jedem zehnten Fall besser aus. Das ist keine grosse Zahl - aber wer") print(" EINMAL misst, hat genau eine Ziehung aus dieser Verteilung.") print("\n Die Lehre ist nicht 'nimm 4.000 Testtage' - die hat niemand. Sie") print(" lautet: Ein Kennzahlenvergleich ohne Angabe seiner Streuung ist keine") print(" Aussage. Bei kurzen Zeitraeumen gehoert eine Kreuzvalidierung dazu.") print("\n" + "=" * 84) print(" WAS MAN DARAUS MITNIMMT") print("=" * 84) print("Der Prognostiker optimiert den MSE, der Planer traegt die Kosten - und") print("die beiden Masse zeigen hier in verschiedene Richtungen. Drei Saetze:") print() print(" 1. Sagen Sie nicht den Erwartungswert vorher, sondern die Groesse, die") print(" in die Entscheidung eingeht. Beim Newsvendor ist das das kritische") print(" Quantil - und das kann man direkt schaetzen.") print(" 2. Bewerten Sie Prognosemodelle an den ENTSCHEIDUNGSKOSTEN. Die sind") print(" in Euro und damit vergleichbar; ein MSE ist es nicht.") print(" 3. Ein pauschaler Sicherheitszuschlag ist besser als nichts und") print(" schlechter als ein Modell, das die Unsicherheit selbst aus den") print(" Merkmalen liest.") print() print("Der naechste Schritt - Prognosemodelle so zu trainieren, dass sie die") print("Entscheidungskosten direkt minimieren (Smart Predict-then-Optimize,") print("differenzierbare Optimierungsschichten) - ist Forschungsstand und") print("erfordert Bibliotheken wie cvxpylayers. Die dritte Zeile der Tabelle") print("oben ist seine einfachste denkbare Form: ein einziger Parameter, auf") print("Kosten statt auf Fehler trainiert.") print("=" * 84)