Vorhersageintervalle für XGBoost-Prognosen: Conformal Prediction in R
- Gegeben: Ein XGBoost-Modell prognostiziert die Restlaufzeit von Anlagen. Aus der Prognose soll eine Wartungs- und Alarmentscheidung werden.
- Problem: Die Punktprognose sagt nichts über ihre eigene Verlässlichkeit, und klassische Vorhersageintervalle setzen Modellannahmen voraus, die Gradient Boosting nicht erfüllt.
- Verfahren: Split Conformal Prediction. Ein vom Training unberührtes Kalibrierset liefert das Quantil, das zur Intervallbreite wird. Verteilungsfrei und modellunabhängig.
- Ergebnis im Beispiel: ±30,2 Tage bei 91,0 Prozent Überdeckung über alle Anlagen, im kritischen Segment C aber nur 50,8 Prozent. Gruppenweise Kalibrierung bringt alle drei Segmente in den Zielbereich (0,885 / 0,915 / 0,942).
- Grenze: Die Garantie hängt an der Austauschbarkeit von Kalibrier- und Testdaten. Bei Zeitreihen und zensierten Restlaufzeiten gilt sie nicht ohne Zusatzannahmen.
Ein XGBoost-Modell sagt eine Restlaufzeit von 41 Tagen vorher. Ob 41 auch 12 heißen könnte, steht in dieser Zahl nicht, und für die Wartungsentscheidung ist genau das die relevante Information. Keine zwanzig Zeilen R liefern das fehlende Stück: Conformal Prediction macht aus der Punktprognose ein Vorhersageintervall mit nachprüfbarer Überdeckung. Der Nebeneffekt ist der eigentliche Gewinn. In dem hier gerechneten Beispiel hält das Intervall über alle Anlagen die zugesagten 90 Prozent und deckt im kritischen Segment nur rund die Hälfte der Fälle ab: 50,8 Prozent in diesem Lauf, im Mittel über 200 Wiederholungen 48,6. Sichtbar wird das erst, wenn man je Gruppe nachmisst.
Das Praxisproblem: eine Punktprognose ohne Vorhersageintervall
Prognoseprojekte enden fast immer an derselben Stelle. Das Modell läuft, die Güte ist dokumentiert, die Prognose liegt vor. Dann fragt der Fachbereich, ob man sich darauf verlassen kann. Ein Punktwert kann diese Frage nicht beantworten. Wer eine Instandhaltung terminiert, ein Ersatzteil disponiert oder eine Anlage in der Priorisierung nach oben zieht, entscheidet anhand des ungünstigen Rands und nicht anhand des erwarteten Werts.
Klassische Vorhersageintervalle kommen aus der Modellwelt der linearen Regression: normalverteilte, homoskedastische Fehler. Sobald Gradient Boosting, Random Forest oder ein gestapeltes Modell im Einsatz ist, sind diese Annahmen nicht erfüllt, und die üblichen Intervalle sind entweder nicht verfügbar oder falsch kalibriert. Die verbreitete Notlösung, Quantile aus der Streuung der Residuen im Trainingsset, ist zu optimistisch, weil sie den Fehler dort schätzt, wo das Modell angepasst wurde.
Split Conformal Prediction in einem Satz
Split-Conformal-Prediction kalibriert das Intervall auf einem separaten, vom Training unberührten Kalibrierdatensatz. Das Quantil der dort beobachteten absoluten Fehler wird zur Intervallbreite: verteilungsfrei, modell-agnostisch, mit garantierter durchschnittlicher Überdeckung. Überdeckung ist dabei die Trefferquote des Intervalls. Von 100 Prognosen mit einem 90-Prozent-Intervall sollen 90 den späteren tatsächlichen Wert einschließen, und genau das lässt sich im Nachhinein nachzählen, anders als bei einem Punktwert.
Der Mechanismus, Zeile für Zeile: Kalibrierquantil in R
Die Daten werden dreigeteilt: Training, Kalibrierung, Test. Das Modell sieht nur das Trainingsset. Auf dem Kalibrierset werden die absoluten Abweichungen berechnet, im Sprachgebrauch des Verfahrens die nonconformity scores. Das empirische (1 − α)-Quantil dieser Scores ist die halbe Intervallbreite.
library(dplyr) library(xgboost) # >= 3.0.0, neue x/y-Schnittstelle set.seed(42) idx <- sample(nrow(df)) d <- df |> slice(idx) |> # zufaellige Reihenfolge mutate(rolle = case_when( row_number() <= floor(.6 * n()) ~ "train", # nur das sieht das Modell row_number() <= floor(.8 * n()) ~ "kalib", # vom Training unberuehrt .default = "test")) train <- d |> filter(rolle == "train") kalib <- d |> filter(rolle == "kalib") test <- d |> filter(rolle == "test") fit <- xgboost(x = train |> select(all_of(x_vars)), y = train$y, nrounds = 300, max_depth = 4, learning_rate = .05, nthreads = 1, verbosity = 0, seed = 42) prognose <- function(neu) predict(fit, neu |> select(all_of(x_vars))) s <- kalib |> mutate(score = abs(y - prognose(kalib))) |> # nonconformity scores filter(!is.na(score)) |> # NA wuerde sort() still verschieben pull(score) alpha <- 0.10 k <- ceiling((length(s) + 1) * (1 - alpha)) # endliche Korrektur, nicht quantile(s, .9) q <- if (k > length(s)) Inf else sort(s)[k] # zu kleines Kalibrierset -> kein endliches Quantil out <- test |> mutate(punkt = prognose(test), lo = punkt - q, hi = punkt + q) out |> filter(!is.na(y), !is.na(punkt)) |> # NA sichtbar machen, nicht wegmitteln summarise(ueberdeckung = mean(y >= lo & y <= hi), # Ziel 0.90 n = n(), ausgeschlossen = nrow(out) - n())
Die Zeile mit k ist der Punkt, an dem selbstgebaute Varianten meist schiefgehen: Korrekt ist das ⌈(n+1)(1−α)⌉-kleinste Score und nicht das Stichprobenquantil. Diese Korrektur um den einen Rang ist der Grund, warum die Garantie auch bei kleinen Kalibriersets hält. Sie erklärt zugleich, warum k größer als die Anzahl der Scores werden kann: Bei α = 0,10 braucht man mindestens neun Kalibrierfälle, sonst existiert kein endliches Quantil und die ehrliche Antwort ist ein unendlich breites Intervall. Wer nicht selbst rechnen will, nimmt probably::int_conformal_split() aus dem tidymodels-Ökosystem; die Logik ist dieselbe, eingebettet in Workflows.
Was der Lauf zeigt: empirische Überdeckung gegen die Zusage
Das Beispiel rechnet mit 4.000 simulierten Anlagen, 2.400 für das Training und je 800 für Kalibrierung und Test, deren Restlaufzeit von Beanspruchung, Alter, Historienlänge und Temperatur abhängt. Bei 800 Kalibrierfällen und α = 0,10 ist k = 721; das zugehörige Score beträgt 30,2 Tage. Jede Prognose bekommt also ein Intervall von ± 30,2 Tagen, und auf dem Testset liegen 91,0 Prozent der tatsächlichen Restlaufzeiten darin. Garantiert ist damit die Überdeckungswahrscheinlichkeit: Solange die Scores keine Bindungen aufweisen, beträgt sie exakt k/(n+1) = 721/801 = 0,9001. Bei gerundeten Zielgrößen wie Tagen sind Bindungen realistisch, und dann gilt nur noch die untere Schranke von 0,90. Die auf einem Testset beobachtete Quote streut um diesen Wert. Die 91,0 Prozent liegen innerhalb dieser Streuung. Der konservative Überschuss von 1/(n+1) ist bei 800 Kalibrierfällen vernachlässigbar, bei 50 wäre er es nicht.
Die 91,0 Prozent sind das Ergebnis eines Splits. Über 200 zufällige Wiederholungen desselben Aufbaus liegt die empirische Überdeckung im Mittel bei 0,899 mit einer Streuung von 0,015; 90 Prozent der Läufe landen zwischen 0,872 und 0,924. Die Garantie gilt im Mittel über Splits und nicht für den einzelnen. Wer aus einem Lauf mit 0,88 schließt, das Verfahren sei kaputt, hat die Streuung des Schätzers vor sich und keinen Fehler.
Der Mittelwert verdeckt das Segment: Überdeckung je Gruppe prüfen
Interessant wird es bei den Teilgruppen. Die simulierten Anlagen gehören zu drei Segmenten mit unterschiedlich dichter Sensorhistorie; Segment C ist am schlechtesten prognostizierbar. Mit einem einzigen Quantil für alle liegt die Überdeckung in Segment A bei 0,995, in Segment B bei 0,955 und in Segment C bei 0,508. Über 200 Wiederholungen sind es in Segment C im Mittel 0,486, im ungünstigsten Lauf 0,300. Die marginalen 90 Prozent sind eingehalten; die Anlagen, bei denen die Entscheidung wehtut, sind trotzdem systematisch zu eng eingerahmt.
Fallabhängige Intervallbreite: normierte Scores und conformalized quantile regression
Das Basisverfahren liefert für jeden Fall dieselbe Intervallbreite. Das ist ehrlich, aber informationsarm: Bei einer Anlage mit dichter Sensorhistorie ist die Prognose besser als bei einer neu angeschlossenen, und das sollte das Intervall zeigen. Zwei Erweiterungen führen zu fallabhängiger Breite. Die erste normiert den Score an einer geschätzten lokalen Streuung: Ein zweites Modell prognostiziert den absoluten Fehler, und der Score wird durch diese Prognose geteilt. Das Intervall wird dort breit, wo das Fehlermodell hohe Unsicherheit erwartet. Die zweite und heute üblichere ist conformalized quantile regression: Statt eines Mittelwertmodells werden zwei Quantilmodelle geschätzt, etwa das 5- und das 95-Prozent-Quantil, und die Konformalisierung korrigiert lediglich, wie weit diese Quantile daneben liegen. Das Ergebnis erbt die Form der bedingten Verteilung und behält die Überdeckungsgarantie. In R steht sie als probably::int_conformal_quantile() bereit; wenn kein separates Kalibrierset abzugeben ist, kommt man mit Resampling statt mit einem Split aus, in R über probably::int_conformal_cv(). Diese Familie (CV+/Jackknife+) nutzt die Daten besser aus, ihre theoretische Schranke ist aber schwächer als die des Split-Verfahrens; welche Garantie im konkreten Fall gilt, gehört vor dem Einsatz in der Dokumentation nachgesehen und in jedem Fall empirisch nachgemessen.
Für Gruppenstrukturen wie Anlagentypen, Baureihen oder Standorte gibt es zusätzlich die gruppenweise Kalibrierung, in der Literatur Mondrian Conformal Prediction: je Gruppe ein eigenes Quantil. Ein einziges Quantil für alle Anlagen ist wie eine Konfektionsgröße für die gesamte Belegschaft: Im Mittel passt sie, an den Rändern passt sie niemandem. Im Beispiel liefert das 15,9 Tage für Segment A, 26,9 für B und 77,1 für C, und damit Überdeckungen von 0,885, 0,915 und 0,942 bei 0,903 insgesamt. Das kostet Kalibrierfälle pro Gruppe und hat eine harte Bedingung, die in der Praxis oft übergangen wird: Die Gruppeneinteilung muss eine Funktion der Merkmale sein, gegebenenfalls auch der Ausgabe des auf dem Trainingsset gefitteten Modells, und darf nicht von den beobachteten Zielwerten des Kalibrier- oder Testsets abhängen. Wer die Gruppen nach dem Blick auf die Überdeckung bildet, verliert die Garantie, die er gerade prüfen wollte.
Zwei Annahmen von Conformal Prediction: Austauschbarkeit und ein unberührtes Kalibrierset
Conformal Prediction ist annahmearm, aber nicht annahmefrei. Vorausgesetzt wird Austauschbarkeit von Kalibrier- und Testdaten. Das ist schwächer als „unabhängig und identisch verteilt", aber nicht beliebig schwächer: Verlangt wird, dass die gemeinsame Verteilung von Kalibrier- und Testfällen invariant gegenüber ihrer Reihenfolge ist. Die einzelnen Fälle haben damit dieselbe Randverteilung; verzichtet wird nur auf die Unabhängigkeit der Ziehung. Abhängige Fälle sind vereinbar; eine Verschiebung der Datenlage zwischen Kalibrierung und Anwendung, ein Trend oder ein Regimewechsel nicht. Bei Zeitreihen ist die Annahme systematisch verletzt. Zeitliche Splits statt zufälliger stellen die Austauschbarkeit nicht wieder her. Sie machen die Verletzung sichtbar, statt sie wegzumischen. Es braucht sie trotzdem, dazu ein rollierendes Kalibrierfenster und im Zweifel online kalibrierende Varianten, die die Intervallbreite anhand der zuletzt beobachteten Überdeckung nachziehen (conformalForecast, AdaptiveConformal).
Zweitens darf das Kalibrierset das Modell nie gesehen haben. Jede Form von Leckage macht die Intervalle zu eng: Hyperparameter-Tuning über alle Daten, Feature-Selektion vor dem Split, dieselbe Anlage in Trainings- und Kalibrierset bei genesteten Daten. Bei genesteten Daten wird deshalb auf der Ebene der Einheit gesplittet und nicht auf der Ebene der Beobachtung.
Ein dritter Punkt betrifft speziell Restlaufzeiten: Wenn Anlagen im Beobachtungsfenster noch nicht ausgefallen sind, ist die Zielgröße zensiert, und die beobachteten Fehler auf dem Kalibrierset sind systematisch zu klein. Zensierte Fälle gehören entweder aus dem Kalibrierset heraus oder in ein Verfahren, das die Zensierung modelliert. Standard-Conformal auf zensierten Daten liefert Intervalle, die man nicht verteidigen kann. Und die untere Grenze wird an der Null abgeschnitten: Eine negative Restlaufzeit ist keine Aussage. Ein Sonderfall davon entsteht durch die Alarmregel selbst: Wo sie wirkt und ein Eingriff den Ausfall verhindert, wird die Restlaufzeit zensiert. Die nächste Kalibrierung sieht also Daten, die die Regel mit erzeugt hat. Wer sie produktiv betreibt, hält die Fälle mit Eingriff gesondert fest.
So geht das in den Betrieb: einseitige Kalibrierung für die Alarmregel
In einer Instandhaltungslogik entscheidet die Intervalluntergrenze und nicht die Punktprognose. Ein Alarm wird ausgelöst, wenn die untere Grenze der prognostizierten Restlaufzeit unter den Vorlauf fällt, den Disposition und Sperrung brauchen. Für diese Regel ist das zweiseitige Intervall das falsche Werkzeug: Relevant ist allein die Unterschreitung, also wird einseitig kalibriert, und zwar mit dem vorzeichenbehafteten Score statt mit dem absoluten. Und aus demselben Grund wie oben je Segment: Ein global kalibriertes einseitiges Quantil erbt das Problem des globalen zweiseitigen Intervalls und setzt die Alarmschwelle in Segment C systematisch zu spät, also dort, wo die Entscheidung wehtut. Der Code unten zeigt die Mechanik an einer Gruppe; produktiv läuft er je Gruppe.
# Einseitige Kalibrierung fuer die Alarmregel: nur die Unterschreitung zaehlt s1 <- kalib |> mutate(score = prognose(kalib) - y) |> # positiv = Modell hat ueberschaetzt filter(!is.na(score)) |> pull(score) k1 <- ceiling((length(s1) + 1) * (1 - 0.10)) # 10 % Unterschreitung zugelassen q1 <- if (k1 > length(s1)) Inf else sort(s1)[k1] alarme <- test |> mutate(lo = pmax(prognose(test) - q1, 0), # Restlaufzeit ist nicht negativ alarm = lo < vorlauf_tage)
Der Unterschied ist beträchtlich: Das zweiseitige 90-Prozent-Intervall setzt die untere Grenze 30,2 Tage unter die Prognose. Einseitig gerechnet reichen dafür 17,3 Tage, weil das Risiko nur noch auf einer Seite liegt. Ein Teil dieser Differenz ist allerdings schlicht ein anderes Risikoniveau: Das zweiseitige Intervall lässt bei symmetrischen Fehlern nach unten nur rund 5 Prozent zu, die einseitige Regel 10. Vergleichbar wäre das zweiseitige 80-Prozent-Intervall. Wer die zweiseitige 90-Prozent-Grenze für eine einseitige Frage benutzt, plant trotzdem Vorlauf ein, den die Fragestellung nicht verlangt. In diesem Split lagen 12,9 Prozent der tatsächlichen Restlaufzeiten unter der einseitigen Grenze statt der nominellen 10, also dieselbe Split-Streuung wie oben. Genau das ist der Grund, warum die Überdeckung eine laufend gemessene Kennzahl sein muss und keine einmalige Abnahme.
Die Prognoseunsicherheit wird damit Teil der Entscheidungsregel. In einer R-Shiny-Anwendung für den Fachbereich heißt das: Balken für die Punktprognose, Band für das Intervall, Schwellenlinie darüber, und ein sichtbares Feld für die aktuelle empirische Überdeckung.
Der Kalibrierschritt gehört in die Pipeline und nicht ins Skript: Neu-Kalibrierung bei jedem Lauf, Überdeckung als überwachte Kennzahl über die Zeit, gruppenweise ausgewiesen. Driftet sie unter das Ziel, hat sich meist die Datenlage verschoben. Das ist eine Information, die man haben will, bevor der Fachbereich sie bemerkt.
Checkliste: was zu Conformal Prediction standardmäßig dazugehört
- Dreiteilung Training / Kalibrierung / Test. Das Kalibrierset bleibt vom Training und vom Tuning unberührt.
- Quantil mit der ⌈(n+1)(1−α)⌉-Korrektur, Mindest-Kalibrierumfang geprüft (bei α = 0,10 mindestens 9 Fälle), NA vorher entfernt.
- Empirische Überdeckung auf dem Testset ausgewiesen und nicht nur die nominale; der konservative Überschuss bis 1 − α + 1/(n+1) ist mitgedacht.
- Überdeckung je relevanter Teilgruppe geprüft. Gruppenweise Kalibrierung (Mondrian) nur mit vorab festgelegten Gruppen.
- Bei heteroskedastischen Daten, in denen das Modell bereichsweise deutlich unsicherer ist, fallabhängige Intervallbreite über normierte Scores oder conformalized quantile regression.
- Bei Zeitreihen zeitlicher Split, rollierende Kalibrierung, Austauschbarkeit explizit diskutiert.
- Bei genesteten Daten Split auf Einheitenebene statt auf Beobachtungsebene. Bei Restlaufzeiten Zensierung geklärt und die untere Grenze bei Null abgeschnitten.
- Neu-Kalibrierung als Pipeline-Schritt, Überdeckung als laufende Kennzahl.
Drei Anlässe für ein Gespräch über Prognoseunsicherheit
Drei Anlässe, bei denen dieser Schritt üblicherweise fällig wird: Ein Prognosemodell ist im Einsatz, und der Fachbereich fragt nach der Verlässlichkeit einzelner Werte, es gibt aber nur Gütemaße über alle Fälle. Aus einer Prognose soll eine Schwelle oder eine Alarmregel werden, und die Frage nach dem Puffer wird derzeit über Erfahrungswerte beantwortet. Oder ein bestehendes Modell liefert im Mittel gute Zahlen, in einem Segment aber sichtbar schlechte Entscheidungen, und die Überdeckung ist nie nach Gruppen geprüft worden.
Mundlak-Zerlegung in R: Within- und Between-Effekt in Paneldaten trennen. Wie sich in Paneldaten Wirkung von Selektion trennen lässt.
Ein Prognosemodell, das eine Entscheidung tragen soll: Im 15-minütigen Erstgespräch kläre ich, welche Prüfgröße dazugehört und wie sie in den Betrieb kommt.
15-MIN-TERMIN BUCHEN →Wie Prognosemodelle bei sove.it gebaut, geprüft und in Betrieb genommen werden: Predictive Analytics und produktive R-Anwendungen.
Alle Zahlen stammen aus einem tatsächlich gelaufenen Beispiel mit simulierten Daten (4.000 Anlagen, drei Segmente mit unterschiedlicher Fehlerstreuung, xgboost 3.2, set.seed(42)): keine Projektergebnisse und keine illustrativen Werte. Kalibrierquantil und Überdeckungsrechnung sind der oben gezeigte R-Code; die Aufteilung ist dieselbe Permutation wie im Lauf, in dplyr-Schreibweise. Die Gruppenauswertung, die 200 Wiederholungen des Splits sowie die Quantile je Segment stammen aus derselben Auswertung, sind hier aber nicht abgedruckt; gezeigt ist der Kern und nicht die Schleife darüber.
Quellen
- Vovk, V., Gammerman, A. & Shafer, G. (2005): Algorithmic Learning in a Random World. Springer.
- Vovk, V. (2012): Conditional Validity of Inductive Conformal Predictors. Proceedings of ACML, PMLR 25, 475–490. (Mondrian-Kalibrierung)
- Lei, J., G'Sell, M., Rinaldo, A., Tibshirani, R. J. & Wasserman, L. (2018): Distribution-Free Predictive Inference for Regression. JASA 113(523), 1094–1111.
- Romano, Y., Patterson, E. & Candès, E. J. (2019): Conformalized Quantile Regression. NeurIPS 32.
- Barber, R. F., Candès, E. J., Ramdas, A. & Tibshirani, R. J. (2021): Predictive Inference with the Jackknife+. Annals of Statistics 49(1), 486–507.
- Angelopoulos, A. N. & Bates, S. (2023): Conformal Prediction: A Gentle Introduction. Foundations and Trends in Machine Learning 16(4), 494–591.