Vorhersageintervalle für XGBoost-Prognosen: Conformal Prediction in R

Ein 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. Achtzehn Zeilen R liefern sie: ein Intervall 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 49 Prozent der Fälle ab. Sichtbar wird das erst, wenn man je Gruppe nachmisst.

Das Praxisproblem

Prognoseprojekte enden fast immer an derselben Stelle. Das Modell läuft, die Güte ist dokumentiert, die Prognose liegt vor — und 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, 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.

Das Verfahren 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.

Der Mechanismus, Zeile für Zeile

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(xgboost)   # >= 3.0.0, neue x/y-Schnittstelle

set.seed(42)
n   <- nrow(df); idx <- sample(n)
tr  <- idx[seq_len(floor(.6 * n))]                # Training
cal <- idx[(floor(.6 * n) + 1):floor(.8 * n)]     # Kalibrierung — vom Training unberuehrt
te  <- idx[(floor(.8 * n) + 1):n]                 # Test

fit <- xgboost(x = df[tr, x_vars, drop = FALSE], y = df$y[tr],
               nrounds = 300, max_depth = 4, learning_rate = .05,
               nthreads = 1, verbosity = 0, seed = 42)

pred <- function(rows) predict(fit, df[rows, x_vars, drop = FALSE])

s <- abs(df$y[cal] - pred(cal))                   # nonconformity scores
s <- s[!is.na(s)]                                 # NA wuerde sort() still verschieben
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

fc  <- pred(te)
out <- data.frame(y = df$y[te], fit = fc, lo = fc - q, hi = fc + q)

mean(out$y >= out$lo & out$y <= out$hi)           # empirische Ueberdeckung, Ziel 0.90

Die Zeile mit k ist der Punkt, an dem selbstgebaute Varianten meist schiefgehen: Korrekt ist das ⌈(n+1)(1−α)⌉-kleinste Score, nicht das Stichprobenquantil. Diese Korrektur um den einen Rang ist der Grund, warum die Garantie auch bei kleinen Kalibriersets hält — und 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

Das Beispiel unten rechnet mit 4.000 simulierten Anlagen — 2.400 für das Training, 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. Die Garantie sagt für diese Größenordnung eine Überdeckung zwischen 0,9000 und 0,9012 voraus — 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, 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, nicht einen Fehler.

Der Mittelwert verdeckt das Segment, auf das es ankommt

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.

Konstante Breite ist selten die Wahrheit

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, leistet probably::int_conformal_cv() dasselbe über Resampling statt über einen Split.

Für Gruppenstrukturen — Anlagentypen, Baureihen, Standorte — gibt es zusätzlich die gruppenweise Kalibrierung, in der Literatur Mondrian Conformal Prediction: je Gruppe ein eigenes Quantil. 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 vorab feststehen und darf nicht aus den Daten oder aus der Modellausgabe abgeleitet werden. Wer die Gruppen nach dem Blick auf die Überdeckung bildet, verliert die Garantie, die er gerade prüfen wollte.

Zwei Annahmen, die man kennen muss

Conformal Prediction ist annahmearm, aber nicht annahmefrei. Vorausgesetzt wird Austauschbarkeit von Kalibrier- und Testdaten. Das ist schwächer als „identisch verteilt und unabhängig“ und heißt nicht, dass beide aus derselben Verteilung stammen müssen: Verlangt wird, dass die gemeinsame Verteilung der Fälle invariant gegenüber ihrer Reihenfolge ist — jede Permutation ist gleich wahrscheinlich. Abhängige Fälle sind damit vereinbar, ein Trend oder ein Regimewechsel nicht. Bei Zeitreihen ist die Annahme systematisch verletzt; hier braucht es zeitliche Splits statt zufälliger, 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 — Hyperparameter-Tuning über alle Daten, Feature-Selektion vor dem Split, dieselbe Anlage in Trainings- und Kalibrierset bei genesteten Daten — macht die Intervalle zu eng. Bei genesteten Daten wird deshalb auf der Ebene der Einheit gesplittet, 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.

So nutzen wir das produktiv

In einer Instandhaltungslogik entscheidet die Intervalluntergrenze, 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 — mit dem vorzeichenbehafteten Score, nicht mit dem absoluten.

# Einseitige Kalibrierung fuer die Alarmregel: nur die Unterschreitung zaehlt
s1 <- pred(cal) - df$y[cal]                       # positiv = Modell hat ueberschaetzt
s1 <- s1[!is.na(s1)]
k1 <- ceiling((length(s1) + 1) * (1 - 0.10))      # 10 %% Unterschreitung zugelassen
q1 <- if (k1 > length(s1)) Inf else sort(s1)[k1]

lo <- pmax(pred(te) - 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. Wer die zweiseitige Grenze für eine einseitige Frage benutzt, plant also gut zwei Wochen Vorlauf ein, die die Fragestellung gar nicht verlangt. In diesem Split lagen 12,9 Prozent der tatsächlichen Restlaufzeiten unter der einseitigen Grenze statt der nominellen 10 — dieselbe Split-Streuung wie oben, und der Grund, warum die Überdeckung eine laufend gemessene Kennzahl sein muss und keine einmalige Abnahme.

Die Prognoseunsicherheit wird damit Teil der Entscheidungsregel, statt in einer Fußnote zu stehen. In einer Shiny-Anwendung 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.

Was 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 prüfen (bei α = 0,10 mindestens 9 Fälle), NA vorher entfernen.
• Empirische Überdeckung auf dem Testset ausweisen, nicht nur die nominale — und den konservativen Überschuss bis 1 − α + 1/(n+1) mitdenken.
• Überdeckung je relevanter Teilgruppe prüfen; gruppenweise (Mondrian) nur mit vorab festgelegten Gruppen kalibrieren.
• Bei heteroskedastischen Daten fallabhängige Breite (normierte Scores oder conformalized quantile regression).
• Bei Zeitreihen zeitlich splitten, rollierend kalibrieren, Austauschbarkeit explizit diskutieren.
• Bei genesteten Daten auf Einheitenebene splitten; bei Restlaufzeiten Zensierung klären und die untere Grenze bei Null abschneiden.
• Neu-Kalibrierung als Pipeline-Schritt, Überdeckung als laufende Kennzahl.

Wann das ein Thema für ein Gespräch ist

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.

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, Überdeckungsrechnung und Gruppenauswertung sind der oben gezeigte R-Code.