Christoph Buser: Differentialgleichungen mit zufälligen zeitvariierenden Parametern

Adviser: Hans R. Künsch

March 2003

Zusammenfassung
Biologische Prozesse werden mit Differentialgleichungen beschrieben. Die Annahme, dass die Parameter zeitlich invariant sind, erleichtert das Lösen und wird in der Praxis oft getroffen. Die dadurch entstehenden systematischen Fehler werden in Kauf genommen, solange sie nicht zu gross sind.
In unserem Beispiel haben wir drei Grössen: die Biomasse (Bakterien), das Substrat (Nahrung) und den Sauerstoff. Es handelt sich um Konzentrationen. Messbar ist nur die Sauerstoffkonzentration. Wir rekonstruieren die anderen Grössen aus diesen Messdaten. Dazu arbeiten wir mit einem Glätter, welcher Daten der Zukunft und der Vergangenheit berücksichtigt.

Wir geben die Konstanz der Parameter auf und modellieren diese mit zeitvariierenden stochastischen Prozessen, genauer gesagt mit dem mean-reverting Ornstein-Uhlenbeck Prozess. Das Modell wird flexibler. Der Ansatz ist bayesianisch. Wir suchen nicht die besten Parameter, sondern konstruieren die bedingte Verteilung der Parameter, gegeben die Sauerstoffmessdaten. Das ist nicht in geschlossener Form möglich.
Wir verwenden den Metropolis-Hastings Algorithmus und erzeugen eine Markovkette, welche asymptotisch die gewünschte Verteilung hat. Um zweidimensionale Vorschlagsdichten zu umgehen, arbeiten wir mit dem Gibbs-Sampler, der jeweils einen der beiden Parameter wählt, der neu vorgeschlagen wird.

In der ersten Simulation nehmen wir im Metropolis-Hastings Algorithmus bedingte Orn-stein-Uhlenbeck Prozesse als Vorschläge für die neuen Parameterwerte. Die Daten werden nicht in die Vorschlagsdichte einbezogen. Wir unterteilen das Zeitintervall $[0,T]$ in zufällige Intervalle gleicher Durchschnittslänge und ändern den Parameter nur auf einem solchen Intervall. Das ist notwendig, um vernünftige Akzeptierungswahrscheinlichkeiten zu erhalten.
In der zweiten Simulation benutzen wir die quadratischen Abweichungen der Sauerstoffdaten, um in einem Intervall einen Vorschlag zu konstruieren. Die zusätzliche Information reduziert die Varianz der Vorschlagsdichte. Der Rechenaufwand vergrössert sich.

Während des Verfahrens sind wir mit einem Problem konfrontiert. Solange Substrat vorhanden ist, dominiert der Wachstumsparameter den Sterbeparameter. Dieser Maskierungseffekt erhöht die Unsicherheit bei der Bestimmung des Sterbeparameters im ersten Zeitabschnitt. Die Unsicherheit überträgt sich auf die Hauptprozesse. In beiden Simulationen gelingt es meist gut bis sehr gut, die Verteilungen aller Prozesse zu bestimmen. Probleme des Filters, der nur Messwerte der Vergangenheit verwendet, werden durch den Glätter behoben. Der Glätter bringt mehr Daten in das Verfahren und ist dem Filter vorzuziehen.
Der Algorithmus ist rechenintensiv. Einerseits ist zum Erreichen der stationären Verteilung eine lange Einschwingphase erforderlich. Andererseits verringern wir die Abhängigkeiten in der Markovkette, indem wir nicht jedes Element verwenden. Daher ist eine grosse Anzahl Schritte im Algorithmus notwendig.

Es gibt Varianten der Vorschlagsdichte. Wir verzichten auf den Gibbs-Sampler und arbeiten zweidimensional. Möglicherweise wird so das Zusammenspiel der beiden Parameter besser wiedergegeben und der Maskierungseffekt kompensiert.
Ein anderer Algorithmus versucht, mehr Information aus den Sauerstoffabweichungen zu gewinnen, indem deren Vorzeichen berücksichtigt wird.

Download: Compressed Postscript (671 Kb) / PDF (1769 Kb).