Über die AGU-Konferenz in San Francisco habe ich ja nun schon ein bisschen berichtet, aber hab ich denn da auch was geschafft? Klar doch, deswegen gehts heute um mein Poster. Das Poster umfasste den Inhalt eines Paper, das gerade im Review steckt, und ich hoffe, dass die Reviewer endlich mal ihre Hüfli bewegen und es da weitergeht.
Das Paper ist ein Teil der Dinge, die ich diesen Sommer während meiner drei Monate an der UC Irvine bei Los Angeles mit Jasper Vrugt zusammengestöpselt habe. Die ganze Kiste ist etwas kompliziert und es kommen u.a. vier verschiedene Optimierungsalgorithmen zum Einsatz; aber ich versuche die Kernpunkte zu erläutern. Im Endeffekt ist die Grundidee einfach, und die Studie schön rund und hat eine klare Aussage.
Bayesian Model Averaging
Um erst einmal die Schussrichtung ein wenig einzugrenzen, wir befinden uns hier auf dem Gebiet der Ensemble-Modellierung und befassen uns mit einer Methode, sich sich Bayesian Model Averaging nennt.
Ich habe in einem Artikel in Nature einen schönen Satz von Philip England zu Modellen gefunden:
The danger in creating fully detailed models of complex systems is ending up with two things you don’t understand — the system you started with, and your model of it.
Wenn man ein System modelliert (und ich bin jetzt extra noch eine Zeit lang vage, welches System das bei uns ist, weil diese Aussagen und Methode sehr generell verwendbar sind), dann kann man versuchen die Physik oder sämtliche bestimmenden Vorgänge ganz genau und klein einzubauen. Dabei läuft man aber Gefahr, dass man nicht mehr nur sein System, sondern auch das Modell selbst nicht mehr ganz im Griff hat.
Das heißt nicht, dass komplexe Modelle nicht ihren Platz haben (und meistens arbeite ich ja auch damit), aber für viele Fragen können einfache, konzeptionelle Modelle genau so gute Antworten geben. Diese versuchen, wenige Zustandsvariablen durch möglichst einfache Gleichungen mit wenigen Parametern zu verbinden (und ich hoffe ihr bekommt keine wenn ich Zustandsvariablen ab sofort als ‘Zustände’ abkürze). Aus dem Vergleich mit Messdaten kann man das Modell (bzw. die Parameter) kalibrieren, und mit den Randbedingungen für die Zukunft laufen lassen.
Aber Modellierer heißen auch die Unsicherheiten willkommen, und versuchen sie nicht nur in Unsicherheiten der Parameter oder Zustände der Modelle zu würdigen, sondern auch in einer Nachbearbeitung der Ergebnisse. Und auch die Tatsache, dass ein einzelnes Modell niemals in jedem Bereich gültig sein wird, wird erkannt. Stattdessen versucht man oft, ein ganzes Rudel Modelle (Ensemble) auf die Wirklichkeit loszulassen. Die Modelle sind dabei alle leicht unterschiedlich, z.B. mit unterschiedlicher Struktur (z.B. lineare gg. nichtlineare Gleichungen), mit unterschiedlicher Parametrisierung oder variierten Start- oder Randbedingungen. Die Idee der Modellmittelung ist dann, dass man eine bessere Aussage gewinnen kann, wenn man mehrere Modelle in irgendeiner Form mittelt. Denn selbst wenn ein Modell insgesamt schlechtere Ergebnisse liefert, mag es z.B. in einem Zeitraum dem besten Modell Informationen hinzufügen, die die unweigerlich vorhandenen Grenzen des besten Modells ergänzen.
Eine spezielle Variante der Mittelung ist das Bayesian Model Averaging (BMA). Es folgt einem Bayes-Ansatz, der Wahrscheinlichkeiten verrechnet. Man hat bei Bayes immer eine erste Annahme einer Wahrscheinlichkeit (a-priori-Verteilung). Wollen wir z.B. einen Parameter optimieren, kann das z.B. ein – vermutlich recht breiter – geschätzter Rahmen sein, in den der Parameter fallen sollte. Der a-priori-Teig wird dann durch ein Förmchen (“likelihood”) gequetscht, das angibt, wie wahrscheinlich es ist, dass unser Modell gerade die Messungen beschreibt, genauer welcher Teil unserer a-priori-Verteilung näher an der Wahrheit ist. Hinten kommt dann eine (hoffentlich passender geformte) neue Wahrscheinlichkeitsverteilung raus, die a-posteriori-Verteilung.
Im Falle des BMA liefert jedes Modell in unserem Ensemble eine Vorhersage über die Entwicklung des Systems. Für die Mittelung machen wir daraus eine Wahrscheinlichkeitsverteilung (wie das geschieht, bzw. welche Annahme man dort macht und warum diese nicht stimmt, ist der Inhalt dieser Studie), die unsere beste Näherung an die Wirklichkeit repräsentieren würde, wäre denn gerade dieses Modell tatsächlich das beste im Ensemble. Diese Verteilung nennt man die konditionelle pdf gk des Modells k – ich verwende jetzt die englischsprachige Kurzform pdf für die Wahrscheinlichkeitsdichtefunktion. Was wir dann bestimmen müssen, sind Gewichte wk. Diese stellen die (a-posteriori-)Wahrscheinlichkeit dar, dass k auch wirklich das beste Modell ist.
Im Kurzdurchlauf haben wir viele Vorhersagen vieler Modelle die wir zu einer gemeinsamen Vorhersage kombinieren wollen. Wir geben dazu jedem Modell ein Gewicht, das sagt wie sehr dieses Modell zur gemeinsamen Vorhersage beitragen darf. Zu jedem Zeitpunkt können wir dann ein gewichtetes Mittel der einzelnen Vorhersagen bilden.
Und auch noch mal als Bild (das ist aus Vrugt&Robinson, 2007, doi:10.1029/2005WR004838 geklaut):
In gestrichelt sieht man die konditionelle pdf für die drei einzelnen Modelle. Bitte merkt euch schonmal, dass diese Gauß-Verteilungen um den vorhersagten Wert bilden. Aus dem Vergleich mit der echten Messung (über die ganze Zeitserie) kann man die Gewichte für die Modelle abschätzen, und die durchgezogene Linie bildet dann die BMA-Vorhersage aus den Einzelmodellen.
Wasserabfluss
Jetzt lasst uns mal die Modelle etwas genauer definieren. Bislang war das nicht so nötig – denn BMA kann man überall antreffen. Das gibt es als (besser?) Alternative zu p-Werten in Auswertungen medizinischer Studien, in der Finanzwelt, ich habe es auch schon gesehen als Methode, um verschiedene Studien zu kombinieren, die die abschreckende Wirkung der Todesstrafe untersucht haben. In der Geowissenschaft trifft man aber auch oft Ensemble-Modelle. Wenn der Wetterbericht z.B. von einer Regenwahrscheinlichkeit von 40 % spricht, dann heißt das nichts anderes, als dass ein Vorhersage-Ensemble läuft (in diesem Fall gleiches Modell, verschiedene Startbedingungen), und dass 40 % der Modellläufe Regen vorhersagen.
Adrian Raftery hat das BMA in der Form, mit der wir uns befassen, in die Meteorologie eingebracht, und jetzt wird es auch fleißig in der Hydrologie getestet. Ich habe hier ein wenig meinen vertrauten Bereich der Geophysik verlassen und mich in den Bereich der Abflussmodellierung von Einzugsgebieten begeben. Macht aber nichts, wir haben mit einem vertrauten historischen Datensatz vom Leaf River in Mississippi zu tun. Die Fragestellung solcher Modelle ist einfach: Wir haben das Einzugsgebiet eines Flusses, wissen wie viel es regnet und wie viel verdunstet -> wie viel Wasser fließt unten heraus.
Das könnte man jetzt komplex modellieren: Als physikalisches Modell, das den Boden modelliert, die Geländeneigung und die Oberfläche zur Atmosphäre. Diese Fragestellung lässt sich aber auch einfacher mit konzeptionellen Modellen bearbeiten. Diese haben nur ganz wenige Zustände, die quasi in einem Eimer alles sammeln, was an Wasser im Einzugsgebiet vorhanden ist. Die Gleichungen geben dann an, wie schnell das unten wieder heraus fließt (oder auch nicht); und man kann noch verschiedene Eimer für schnellen und langsamen Abfluss haben. Wenige Zustände (4 bis 8), eine oder zwei Hände voll Parameter, die man an Messungen anpasst, um dann eine Maschine zu haben die aus Niederschlagsprognosen Abflüsse vorhersagt.
Man kann also nun ein Ensemble an Modellen basteln, das ein wenig andere Maschinenteile verwendet (z.B. in dem es unterschiedliche Mengen Eimer hinstellt) oder einige Parameter festhält (und so beispielsweise nichtlineare Gleichungen ausschaltet). Gerrit Schoups hat ein flexibles Modell gebastelt, dass er hModel nennt, mit dem man so einfach ein Ensemble an Modellen erstellen und auf den gleichen Datensatz loslassen kann. Letztendlich haben wir eine leicht andere Methode verwendet für die Studie, aber ich habe beides getestet, die Aussage ist prinzipiell die gleiche. Es hat jetzt ‘historische’ Gründe, wie wir es durchgeführt haben. Und damit sind wir endlich am Beginn der Studie…
Konditionell und glockenförmig – oder?
Erinnert ihr euch noch an das Bild oben? Die konditionelle pdf war die Wahrscheinlichkeitsverteilung der vorhergesagten Abflussmenge unter der Bedingung, dass dieses (einzelne) Modell auch tatsächlich das beste wäre. Da aus so einem einfachen Modell zunächst nur ein einzelner Wert herausfällt und keine pdf (deterministische oder Punkt-Vorhersage), können wir in dem BMA-Prozess, der die Gewichte herausfindet, auch gleich noch Varianzen mitschätzen lassen, also die pdf als Gauß-Verteilung annehmen.
Diese traditionelle Annahme haben wir überprüft – und daher eine Methode eingesetzt, die mit einem einzelnen Modell nicht nur Punktwerte, sondern beliebig geformte Wahrscheinlichkeitsverteilungen ermitteln kann. Außerdem ist die Annahme auch, dass die Breite dieser Gauß-Verteilung sich nicht mit der Zeit ändert – auch dies wollen und können wir überprüfen. Wir verarbeiten dazu die Modelle in Partikelfiltern.
Aber gehen wir der Reihe nach vor – hier ist der Fahrplan der Studie:
Der erste Schritt ist die Kalibrierung. Wir wollen uns im folgenden Prozess nicht auch noch mit der Optimierung der Parameter herumschlagen müssen, daher fixieren wir diese zuerst. Und als Bonus können wir das direkt noch als alternative Methode einsetzen, ein Ensemble zu bilden.
Zur Kalibrierung setzen wir den bewährten SCE-UA-Algorithmus (Shuffled Complex Evolution – University of Arizona, Duan et al. 1992) ein, der uns ein Punkt-Set an Parametern gibt. Wir verwenden nur eine Modellvariante, setzen aber folgenden Trick ein, um an ein Ensemble zu gelangen: Wir führen die Kalibration nicht an allen Messdaten (1096 Tageswerte des Abfluss) durch, sondern nur an einem Teil. So haben wir dann vier Modelle, die nur an (i) niedrigen Abflussmengen, (ii) mittleren Mengen, (iii) hohen Mengen oder (iv) doch allen Daten kalibriert sind. Modell (i) z.B. sollte dann z.B. besonders gut passen wenn es wenig Niederschlag gibt. Das Ergebnis sehen wir hier:
Oben in der Abbildung steckt die antreibende Kraft, der Niederschlag. Die Linien darunter zeigen exemplarisch für 100 Tage in 1953 die vorhergesagten Abflussmengen für die verschiedenen Modelle, und als Punkte die tatsächlich gemessenen Werte. Man sieht, dass die meisten Modelle ziemlich bescheidene Ergebnisse liefern, aber z.B. passt das “low”-Modell gut wenn es eben nicht regnet, und könnte wertvolle ergänzende Information liefern. Es war das Ziel, ein möglichst extremes Ensemble zu testen. Genauso kann man mit verschiedenen Modellstrukturen ein “braveres” Ensemble erzeugen, aber für die Aussage der Studie ist das nicht relevant.
Der entscheidende Schritt folgt als nächstes, wenn wir mit den so kalibrierten Modellen von vorne anfangen, aber nicht einzelne Modellläufe ansehen, sondern für jedes Modell einen Partikelfilter laufen lassen.
Partikelfilter
Wie genau Partikelfilter funktionieren und wie man die Methode zur Parameterschätzung einsetzen kann, habe ich hier bereits ausführlicher beschrieben. Hier will ich nur die wichtigsten Ideen zusammenfassen. Da wir die Parameter ja schon kalibriert haben, müssen wir diese natürlich nicht mehr mitschätzen, und können uns auf die ursprüngliche Aufgabe von Filtern beschränken: Das Anpassen von Zuständen in Konfrontation mit Messdaten.
Hier ist nochmal der prinzipielle Ablauf eines Partikelfilter-Schritt. Ein Partikel kann man wie ein Rennauto ansehen, das einen Satz möglicher Werte für die Zustände unseres Systems transportiert, in unserem Fall ist ein Zustand das Maß, wie sehr ein Eimer gefüllt ist. Aus den Zuständen kann man dank des Modells eine Vorhersage ableiten, wie groß der Wasserabfluss sein wird. Im Vergleich mit den Messungen stellt sich dann heraus, welches Rennauto von besseren Zustandsschätzungen gesteuert wird: Bei besserer Übereinstimmung wird das Auto schneller die Runde beenden. Wir würdigen das, indem wir dem Partikel ein großes Gewicht geben.
Nachdem alle Partikel gewichtet sind, verteilen wir sie neu, und zwar so dass dort, wo vorher ein schweres Partikel saß, jetzt die Wahrscheinlichkeit höher ist mehrere Partikel vorzufinden, die aber alle wieder gleich schwer sind. Beide Darstellungen, die unterschiedlich schweren oder die unterschiedlich dicht hockenden Partikel, sind eine diskrete Darstellung der Wahrscheinlichkeitsverteilung für einen Zustand!
Als nächstes muss dann die nächste Runde absolviert werden, das heißt das Modell wird auf die Zustände jedes Partikels losgelassen, so lange bis am Rundenende die nächste Messung bereit steht.
Der Vorteil des Partikelfilters ist, dass wir beliebig geformte pdfs erhalten können, die sich auch zeitlich ändern dürfen. Im Vergleich zu den deterministischen Modellen, denen wir erst eine Gauß-Verteilung zuschustern müssen, liefert der Filter also direkt die gesamte konditionelle pdf. Und so schaut das dann aus:
Für jedes der Modelle gibt es jetzt eine Vorhersage mit Unsicherheitsbereich. Man sieht vor allem, dass auch die schlechten Punkt-Vorhersagen bereits stark verbessert wurden. Das ist die Stärke des Partikelfilters. Aus dem Bewegungsspielraum, den man durch die definierten Fehlermodelle gibt, schöpft er die Fähigkeit, die Zustände zu verbessern und die Vorhersagen näher zu den Messungen zu bringen. Und das ohne die Parameter zu verändern, es ist also auch eine Art der Nachbearbeitung.
Der Filter liefert folglich sofort eine konditionelle pdf, die das kann, was wir überprüfen möchten: zeitlich veränderliche, nicht-Gauß-Verteilungen. Falls die Annahme des BMA nicht stimmt, sollten wir also eine Abweichung der Filter-pdf von einer Gauß-Verteilung sehen und/oder zeitliche Veränderungen. Und genau das tun wir!
Histogramme und Gauß-Mischungen
Wir können uns einfach mal für ein Modell zu verschiedenen Zeitpunkten die Vorhersagen der verschiedenen Partikel herausgreifen und als Histogramm (blaue Balken) auftragen:
Wir sehen an der Form der Wahrscheinlichkeitsverteilungen sofort: Diese sind eben NICHT immer Gauß-verteilt und schon gar nicht unveränderlich in der Zeit. Diese klassische Annahme des BMA ist also falsch. Wie groß aber konkret der Einfluss dieser ungültigen Annahme ist, das sehen wir uns im folgenden an. Dazu müssen wir die Vorhersagen der Partikelfilter aber als konditionelle pdf in den BMA-Prozess hineinbekommen. Dazu brauchen wir eine geschlossene Form der pdf, also eine Funktion, wir haben aber eine offene Form, denn die Partikel geben uns ja nur die Werte an einigen Stellen an.
Die Lösung? Wir passen ein Mischung von Gauß-Verteilungen an das Histogramm an (mit dem dritten Optimierungsalgorithmus, einer angepassten “Cross-Entropy-Methode”, beschrieben und implementiert 2004 durch Botev und Kroese). Dabei addiert man mehrere Gauß-Verteilungen, die an unterschiedlichen Psoitionen sitzen, unterschiedlich breit sein und unterschiedlich stark gewichtet werden können. Dadurch kann man letztendlich jegliche Form einer Verteilung gut annähern, wie man an den roten Linien in obigem Bild sieht.
Dreimal Mitteln, bitte
Um jetzt letztlich den Einfluss der Methoden auseinanderbasteln zu können, haben wir drei verschiedene BMA-Methoden eingesetzt, die im obigen Fahrplan aufgeschlüsselt sind: Als “calBMA” bezeichnen wir es, nur die Punktwerte aus der Kalibration zu nehmen, mit einer Gauß-Verteilung zu versehen und so in ein BMA einfließen zu lassen. Das ist der klassische Weg. Am Ende unseres Prozesses erhalten wir direkt aus den Partikelfiltern konditionelle, zeitlich veränderliche pdfs, die wir mit Gauß-Mischungen in eine geschlossene Form für das BMA bringen, das ist als “mixBMA” bezeichnet. Um den Einfluss des Partikelfilters von der Annahme der Gauß-verteilten pdf zu trennen, gibt es auch noch den Mittelweg “medBMA”, bei dem man die mittlere Vorhersage des Filters als Punkt-Vorhersage nimmt und mit einer Gauß-Verteilung versieht.
Die Ergebnisse sind in der folgenden Tabelle zusammengefasst:
In der linken Hälfte stehen die Gewichte der einzelnen Modelle für jeden der BMA-Prozesse. Letztlich haben wir hier das (vierte) Optimierungsverfahren DREAM eingesetzt, um die Gewichte zu optimieren. Man sieht, dass außer beim medBMA nur ein Modell das meiste Gewicht vereint, wenig überraschend ist es auch das, welches auf alle Messungen kalibriert war.
Aber interessanter ist die rechte Hälfte: Hier ist der mittlere Fehler der Vorhersagen aufgeführt. Zunächst für jede Einzelvorhersage (also für die verschiedenen Arten, zur konditionellen pdf zu kommen), und darunter für die gemittelte Vorhersage aus BMA. Zunächst fällt auf, dass das BMA-Modell gar nicht unbedingt den niedrigsten Fehler aufweist. In diesem Fall liegt das aber daran, dass wir in DREAM nicht diesen Fehler optimiert haben, sondern die “log-likelihood”. Diese ist dann tatsächlich für das BMA-Modell jeweils am besten, aber der mittlere Fehler lässt sich einfach besser darstellen und mit früheren Studien vergleichen (wohl wissend dass es sich eigentlich um ein Maß mit Gaußschen Annahmen handelt…).
Uns interessieren jetzt nur die Abfolge der Fehler für “BMA” über die drei Ansätze. Die Verwendung des Partikelfilters, der ja die Zustände verbessern kann, führt zu einer deutlichen Verbesserung von 21,3 % auf 17,9 % Fehler bei medBMA. Der weitere Schritt ist dann der Einfluss durch die fehlerhafte Annahme einer Gauß-Verteilung der konditionellen pdf bei medBMA. Man sieht, dass die Annahme den Fehler von 15,2 % (mixBMA) auf die 17,9 % bei medBMA steigert, also durchaus nochmal ein relevanter Anteil.
Wir haben also gesehen, das Partikelfilter einsetzbar sind, um für Bayesian Model Averaging zeitlich variable, beliebig geformte Wahrscheinlichkeitsverteilungen zu erhalten. Diese erlauben, die ungültige Annahme von zeitlich konstanten Gauß-Verteilungen der Vorhersagen zu ersetzen und führen zu besser sitzenden BMA-Modellen.
Kommentare (2)