Potentielle Verallgemeinerungen

kri-tischen Parameterregionen analytisch belegt werden. Dazu wird angestrebt, analog zur Untersuchung des Maßwechsels von Parekh und Walrand in Glasserman und Kou [GK95]

unter Verwendung der Theorie großer Abweichungen Parameterregionen zu identifizieren, in denen unsere Maßwechsel asymptotisch optimal oder zumindest effizienter als der von Parekh und Walrand sind.

5 Adaptives Importance Sampling

Die Schwierigkeiten, den optimalen IS–Sch¨atzer oder Maßwechsel, einen zumindest asym-ptotisch optimalen oder wenigstens einen solchen, der große Varianzreduktion liefert, zu finden, motivieren den Ansatz, den Maßwechsel w¨ahrend der Simulation adaptiv zu ver-bessern. Adaptiv bedeutet anpassend oder auch lernend^∗, und ein adaptives Verfahren ist ein iteratives Verfahren, das basierend auf Zwischenergebnissen angepaßte verbesserte Ergebnisse produziert. Beim Importance Sampling soll so – zumindest f¨ur eine gegebe-ne Klasse von Strategien – der optimale Maßwechsel

”gelernt“ oder angen¨ahert werden.

Adaptives Importance Sampling basiert dabei darauf, daß innerhalb der Simulation die gesuchte Gr¨oße, die Varianz und/oder weitere Gr¨oßen, die in Bezug zur G¨ute des Sch¨atzers stehen, gesch¨atzt werden k¨onnen. Der Maßwechsel wird w¨ahrend der Simulation basierend auf den bereits durchgef¨uhrten Simulationsschritten angepaßt, um so den korrespondie-renden Sch¨atzer zu verbessern.

Eine Maßwechselstrategie ist typischerweise dadurch charakterisiert, daß Parameter des zugrundeliegenden Wahrscheinlichkeitsmaßes ver¨andert werden. Dieses Wahrscheinlich-keitsmaß ist, insbesondere f¨ur mehrdimensionale Zufallsgr¨oßen oder stochastische Prozes-se, nicht nur von einem einzelnen oder einigen wenigen Parametern abh¨angig, sondern beschrieben durch eine ganze Reihe, einen Satz von Parametern, die ver¨andert werden k¨onnen. Dies sind beispielsweise bei Markovketten alle ¨Ubergangswahrscheinlichkeiten oder -raten, wobei das IS–Maß sogar durch eine noch gr¨oßere Anzahl von Parametern als das Originalmaß bestimmt sein kann, wenn man etwa eine inhomogene Markovkette zur Importance–Sampling–Simulation einer homogenen Markovkette verwendet.

Es ist klar, daß eine solche Parametrisierung prinzipiell eine Einschr¨ankung des Maß-wechsels darstellt, in dem Sinne, daß die in Frage kommenden IS–Maße oder –Dichten auf eine parametrisierte Familie eingeschr¨ankt werden, wie etwa auch bei der Skalierung, der Translation und dem exponentiellen Maßwechsel, siehe Abschnitt 3.3.5. Je nach Art der Parametrisierung kann also der optimale Maßwechsel oder sogar jeder asymptotisch optimale Maßwechsel dadurch ausgeschlossen werden, wie es bei den in 3.3.5 beschriebe-nen Strategien auch leider h¨aufig der Fall ist. Somit kann man also auch beim adaptiven Importance Sampling nicht von einem allgemeinen Verfahren sprechen, das immer zum Erfolg f¨uhrt. Die Optimierung, nach welchem Verfahren sie auch immer durchgef¨uhrt wird, erfolgt nur im eingeschr¨ankten Rahmen der gew¨ahlten Parametrisierung. Dies ist jedoch kein spezifisches Problem adaptiver Strategien, es stellt sich ebenso bei anderen Varianten, und zudem ist das Konzept der Parametrisierung sehr flexibel. Grob gesprochen, je gr¨oßer die Anzahl der Parameter, desto geringer ist die Einschr¨ankung des Maßwechsels. Man denke zum Beispiel an den Unterschied zwischen zustandsabh¨angigen und

–unabh¨angi-∗Die BegriffeAdaption undLernensind vor allem in der k¨unstlichen Intelligenz und bei neuronalen Netzen ¨ublich und weit verbreitet. Auch bei der Mustererkennung sind Methoden des ¨uberwachten und des un¨uberwachten Lernens (supervised/unsupervised learning) von großer Bedeutung.

gen Maßwechseln bei Warteschlangenmodellen oder an die oben angef¨uhrte M¨oglichkeit inhomogener Markovketten als IS–Maß f¨ur homogene Markovketten. Auch im Kontext von Tail–Wahrscheinlichkeiten ist die Freiheit beim Maßwechsel groß. Betrachten wir als einfachen Fall exponentiell verteilte reelle Zufallsvariablen. Die naheliegendste Familie f¨ur die IS–Dichte ist sicherlich die Familie der exponentiell verteilten Zufallsvariablen, also schlicht die Ver¨anderung des einen Parameters, der die Exponentialverteilung beschreibt, was in diesem Fall sowohl der Translation als auch dem exponentiellen Maßwechsel ent-spricht, jedoch, wie wir in Abschnitt 3.3.3 gesehen haben, den optimalen Maßwechsel aus-schließt. Weniger eingeschr¨ankte Maßwechsel k¨onnen hier grunds¨atzlich jede reelle Dichte erlauben, insbesondere auch solche, die von mehr als einem Parameter abh¨angen. Wir sehen also, daß eine geschickte, ausreichend allgemeine, m¨oglichst viele M¨oglichkeiten umfassende Parametrisierung das Problem der Einschr¨ankung des Maßwechsels erheblich mindern kann, allerdings selbstverst¨andlich auf Kosten eines erh¨ohten Aufwands.

Die G¨ute einer Maßwechselstrategie und damit die G¨ute des korrespondierenden Sch¨atzers h¨angt dann wesentlich von den konkret gew¨ahlten Parameterwerten f¨ur das IS–Maß oder die IS–Dichte ab. Optimale Parameter f¨ur eine Strategie sind dabei solche, die die Va-rianz f¨ur diese Strategie minimieren und somit den bestm¨oglichen Sch¨atzer liefern. Die Minimierung der Varianz kann dabei sowohl direkt erfolgen als auch indirekt, indem ein der Varianz verwandtes G¨utekriterium optimiert wird. Im Bereich der Optimierung wird ublicherweise von einer zu optimierenden Kostenfunktion oder Zielfunktion gesprochen.¨ Ist diese Zielfunktion exakt verf¨ugbar, so k¨onnen Methoden der Variationsrechnung an-gewendet werden. Ebenso kann ein exponentieller Maßwechsel optimiert werden, wenn Ratenfunktionen großer Abweichungen f¨ur das zu untersuchende System bekannt sind. In der Regel sind aber weder die Varianz oder andere G¨utekriterien f¨ur IS–Sch¨atzer noch Ratenfunktionen bekannt, und selbst, falls dies der Fall ist, ist das resultierende Optimie-rungsproblem f¨ur realistische Systeme analytisch oder numerisch meist nur sehr schwer oder gar nicht l¨osbar. Ein Ausweg ist die Optimierung von Sch¨atzungen des G¨ utekriteri-ums, also etwa die Minimierung der Stichprobenvarianz des IS–Sch¨atzers. Das bedeutet, beim adaptiven Importance Sampling werden mehrere Iterationen durchgef¨uhrt, wobei jede Iteration aus einer (großen) Anzahl von Simulationsl¨aufen – unabh¨angige Replika-tionen, Batches oder regenerative Zyklen – besteht. Aus diesen Simulationsl¨aufen werden dann sowohl die zu bestimmende Gr¨oße als auch die G¨utekriterien gesch¨atzt, basierend auf den Sch¨atzungen werden die Parameter ver¨andert, so, daß sie die Sch¨atzung des G¨ utekri-teriums optimieren, und mit diesen Parametern wird die n¨achste Iteration durchgef¨uhrt.

Die Ver¨anderung der Parameter erfordert dabei nat¨urlich einen Algorithmus, der gem¨aß der Sch¨atzungen aus einer Iteration die neuen Parameter f¨ur die n¨achste Iteration be-rechnet. Ein solcher Algorithmus kann eine einfache Heuristik sein, von der man an-nimmt, daß sie den Sch¨atzer sukzessive verbessert, oder eine systematische Optimierung der Stichprobenvarianz oder des verwandten G¨utekriteriums. Dies ist typischerweise ein

mehrdimensionales nichtlineares Optimierungsproblem, wof¨ur h¨aufig stochastische Algo-rithmen eingesetzt werden. Es existieren viele verschiedene OptimierungsalgoAlgo-rithmen, die auf statistischen Sch¨atzungen operieren, so etwa die klassischen Verfahren von Robbins und Monro [RM51] sowie Kiefer und Wolfowitz [KW52] und andere stochastische Opti-mierungsverfahren wie Simulated Annealing, stochastische Gradientenabstiegs– und New-tonverfahren. Man kann diese auf verschiedene Arten klassifizieren. Eine Unterscheidung ist die nach globalen und lokalen Optimierungsverfahren, also danach, ob globale oder eventuell nur lokale Optima bestimmt werden. Eine weitere Unterscheidung ist die nach Verfahren, die Ableitungen bzw. Gradienten einbeziehen, und solchen, in denen dies nicht der Fall ist. F¨ur ausf¨uhrliche Darstellungen stochastischer Optimierungsverfahren und Optimierung mittels Simulation verweisen wir auf die zahlreiche Literatur, insbesondere auf [And98, Duf97, Fu94, Gly89, KC78, Rub86, RSU96, RS93, Rub97, RM98]. Verfah-ren zur Sch¨atzung von Ableitungen und Gradienten, insbesondere unter Verwendung von Likelihood–Quotienten, wie sie etwa auch bei der Sensitivit¨atsanalyse eingesetzt werden, sind daneben unter anderem in [Ecu90a, Gly90, Nak95] beschrieben.

Adaptives Importance Sampling bietet den offensichtlichen Vorteil, daß nicht vorab ein Maßwechsel gew¨ahlt werden muß, unter dem dann die gesamte Importance–Sampling–

Simulation durchgef¨uhrt wird. Es ist potentiell allgemeiner anwendbar und nicht etwa auf spezielle Typen stochastischer Prozesse oder auf bestimmte Dichtefamilien zugeschnit-ten, das Prinzip ist unabh¨angig von der konkreten Problemstellung, und es ist prin-zipiell kein spezielles Wissen ¨uber das konkrete zu untersuchende System erforderlich.

Die adaptive Verbesserung erleichtert insbesondere auch die Anwendung von Importance Sampling durch Nichtexperten, da der Anwender den Maßwechsel, unter dem schließlich die Sch¨atzung erfolgt, nicht selbst angeben muß. Die Vorteile des adaptiven Importan-ce Sampling gehen auf Kosten eines zum Teil erheblichen Mehraufwands f¨ur die Simu-lation, da die eigentliche Sch¨atzung der gesuchten Gr¨oße erst mit den zuvor simulativ bestimmten Maßwechselparametern erfolgt. Wie ganz allgemein bei iterativen Verfahren stellt sich auch hier die Frage nach der Konvergenz, nach der Wahl geeigneter Startwerte und Abbruchkriterien. Diese Punkte k¨onnen je nach Art des zu untersuchenden Systems, der Parametrisierung, der G¨utekriterien zur Bewertung der iterativen Sch¨atzwerte und des verwendeten Optimierungsverfahrens sehr unterschiedlich sein. Zudem sind die ge-fundenen Parameter aufgrunddessen, daß die Optimierung auf Simulationen beruht, nur ann¨ahernd optimale Sch¨atzungen, was jedoch bei hinreichend guten Sch¨atzungen in Kauf zu nehmen ist. F¨ur realistische Systeme ist meist eine analytische oder numerische Bestim-mung eines geeigneten Maßwechsels nicht m¨oglich, auch die Theorie großer Abweichungen ist dann nicht anwendbar. Vor allem in solchen F¨allen ist adaptives Importance Sampling eine gute Alternative und stellt oft die einzige praktikable M¨oglichkeit dar.

Adaptive Importance–Sampling–Strategien unter Verwendung einfacher Heuristiken wur-den Ende der 80er Jahre im Bereich der Sicherheit physikalischer Systeme zur Bestim-mung von Fehlerraten mechanischer Strukturen entwickelt und dann auch bei der

Berech-nung hochdimensionaler Integrale [OB92, OB93] sowie zur Bestimmung kleiner Bitfeh-lerwahrscheinlichkeiten in digitalen Kommunikationskan¨alen eingesetzt [SR93]. Aktuelle-re VerfahAktuelle-ren zur Bestimmung von Bitfehlerwahrscheinlichkeiten und VerfahAktuelle-ren zur Be-stimmung von ¨Uberlauf– oder Verlustwahrscheinlichkeiten in Warteschlangenmodellen, die wir in Abschnitt 5.3 vorstellen, verfolgen in der Regel den Ansatz der Minimierung der Stichprobenvarianz des IS–Sch¨atzers, wobei stochastische Suchverfahren, Simulated Annealing, stochastische Gradientenabstiegs– oder Newtonverfahren eingesetzt werden [ADT93, DT93a, DT93b, DAT93, DAFT93, FDT94, FDAT96, RSN⁺00, Sri98, Sri02].

Eine Alternative wird in [OA91] vorgeschlagen. Hier wird die IS–Dichte adaptiv basie-rend auf einem Abstandsmaß f¨ur Wahrscheinlichkeitsverteilungen, der sogenannten f–

Divergenz oder Ali–Silvey–Distanz, bestimmt, wobei leider eine ungeschickte Parametri-sierung zum Ausschluß der optimalen IS–Dichte f¨uhrt. Dieser konkrete Ansatz wurde nach unserem Wissen nicht weiter verfolgt, ist aber dennoch bemerkenswert, da die Ver-wendung eines Abstandsmaßes gerade die zentrale Idee bei der Cross–Entropy–Methode [Rub97, LR00, Kro02] ist, die in j¨ungster Zeit in diversen Anwendungsbereichen große Bedeutung erlangt hat, und die wir in Abschnitt 5.4 genauer beschreiben. In [OK00]

werden adaptive Importance–Sampling–Verfahren zur Analyse sogenannter Bayes–Netze pr¨asentiert, wobei stochastische Gradientenabstiegsverfahren sowohl zur direkten Vari-anzminimierung als auch zur Minimierung von Abstandsmaßen angewendet werden. Die-se Arbeit ist besonders bemerkenswert, da als ein solches Abstandsmaß die Kullback–

Leibler–Divergenz, also die Cross–Entropy, vorgeschlagen wird, jedoch eben gerade unter einer anderen Bezeichnung und offensichtlich ohne Kenntnis der entsprechenden Arbeiten, die die Bezeichnung Cross–Entropy verwenden^†.

5.1 Einfache heuristische Verfahren

Die ersten einfachen heuristischen adaptiven Importance–Sampling–Strategien wurden wie bereits erw¨ahnt Ende der 80er Jahre im Bereich der Sicherheit physikalischer Sy-steme zur Bestimmung von Fehlerraten mechanischer Strukturen entwickelt und einge-setzt. Stadler und Roy [SR93] beschreiben und analysieren vier solche fr¨uhen Algorith-men und ihre Anwendung auf die Sch¨atzung kleiner Bitfehlerwahrscheinlichkeiten in di-gitalen Kommunikationskan¨alen. Die allgemeine Problemstellung in diesem Kontext ist, daß empfangene Signale verf¨alscht oder verrauscht (noisy) sind. Diese Verf¨alschung ist typischerweise gem¨aß einer mehrdimensionalen reellen Dichte verteilt, h¨aufig wird eine

†Wir sehen hier einen h¨aufig zu beobachtenden Effekt, der nicht nur beim Importance Sampling son-dern viel allgemeiner in der Forschung auftritt: verschiedene Disziplinen oder Zweige mit unterschiedlichen Anwendungsgebieten verwenden ohne voneinander Notiz zu nehmen ¨ahnliche oder gar identische Tech-niken, meist aufgrund traditioneller Festlegungen oder historischer Entwicklungen unter abweichenden Bezeichnungen und Notationen. Leider sind oft die Grenzen zwischen den Anwendungsfeldern und den dort t¨atigen Forschern so starr, daß die eine Seite nicht weiß, was auf der jeweils anderen vor sich geht, und daß Entwicklungen parallel ablaufen und analoge Arbeiten doppelt gemacht werden.

Gaußsche Verf¨alschung (Gaussian noise) angenommen. Beim Importance Sampling wird dann gerade diese Dichte modifiziert, ist also die IS–Dichte.

Bei den in [SR93] beschriebenen Strategien wird nicht ein G¨utekriterium f¨ur den IS–

Sch¨atzer direkt optimiert, sondern der Ansatz verfolgt, den Maßwechsel zu optimieren, indem adaptiv die optimale IS–Dichte approximiert wird. Dazu wird die Tatsache aus-genutzt, daß die optimale IS–Dichte f¨ur die Sch¨atzung von Wahrscheinlichkeiten mittels des Erwartungswertes der Indikatorfunktion des interessierenden Ereignisses gerade die bedingte Dichte unter der Bedingung, daß das Ereignis eintritt, ist, siehe Abschnitt 3.3.3.

Neben zwei nichtparametrischen Ans¨atzen, bei denen das Ziel ist, die gesamte optimale IS–

Dichte ohne die Voraussetzung einer vorgegebenen Dichtefamilie zu sch¨atzen, beschreiben Stadler und Roy auch zwei parametrische Ans¨atze, die beide auf der Translation basie-ren. Wir wollen die parametrischen Strategien kurz skizzieren und diskutieren; f¨ur Details sei auf [SR93] und die dort zitierten Originalarbeiten verwiesen. Betrachtet werden da-bei jeweils mehrdimensionale reelle Dichten f¨ur die Verf¨alschung, wobei keine Gaußsche Verf¨alschung vorausgesetzt wird, die Originaldichte also nicht die einer Normalverteilung sein muß.

Die erste Strategie besteht darin, in der Iteration den bedingten Erwartungswert und die bedingten Kovarianzen unter der Bedingung, daß ein Fehler auftritt, zu sch¨atzen, und dann als IS–Dichte f¨ur den n¨achsten Iterationsschritt die Dichte einer Normalverteilung mit diesen Sch¨atzungen als Parameter zu w¨ahlen. Als Algorithmus ergibt sich das Schema

1) W¨ahle Anfangsvektor µinit.

2) Verschiebe die Originaldichte so, daß der Erwartungswertµinitist, und simuliere mit dieser Dichte als IS–Dichte.

3) Sch¨atze mittels der Beobachtungen, in denen Fehler aufgetreten sind, den bedingten Erwartungswert ˆµund die Matrix ˆK der bedingten Kovarianzen.

4) W¨ahle die Dichte einer Normalverteilung mit Erwartungswert ˆµund Kovarianzma-trix ˆK als IS–Dichte, simuliere und sch¨atze so die Fehlerwahrscheinlichkeit ˆp.

5) Wiederhole 3) und 4) bis ein Abbruchkriterium erf¨ullt ist.

Es f¨allt bei der Betrachtung des Algorithmus sofort auf, daß schon die Wahl eines An-fangsvektors einen kritischen Punkt darstellt. Es muß gew¨ahrleistet sein, daß die dadurch bestimmte IS–Dichte, die durch einfache Translation der Originaldichte resultiert, ausrei-chend viele Fehler hervorruft, um die Sch¨atzung des bedingten Erwartungswertes und der Matrix der bedingten Kovarianzen zu erm¨oglichen. In [SR93] bzw. der dort zitierten Ori-ginalarbeit wird hierzu eine Sensitivit¨atsanalyse vorgesehen, also bereits zur Wahl eines Startwertes ein erheblicher Aufwand in Kauf genommen. Ein weiterer kritischer Punkt ist

dann die Sch¨atzung der bedingten Gr¨oßen selbst. Stadler und Roy zeigen, daß zu unge-naue Sch¨atzungen zur Divergenz des gesamten Algorithmus oder zur Konvergenz gegen eine vom optimalen Wert deutlich abweichende L¨osung f¨uhrt. Ein wesentlicher Nachteil des Algorithmus im Hinblick auf eine Verallgemeinerung ist die Einschr¨ankung auf Nor-malverteilungen f¨ur die IS–Dichten in der Iteration. Diese Einschr¨ankung l¨aßt sich auch nicht einfach dadurch beheben, daß schlicht eine beliebige andere Verteilungsfamilie er-laubt wird, vielmehr stellt sie eine Grundlage f¨ur die Effektivit¨at des Algorithmus dar, denn Normalverteilungen sind direkt durch die Parameter Erwartungswert und Kovari-anzmatrix gegeben, und die Dichten lassen sich somit daraus leicht bestimmen. In allge-meineren F¨allen lassen sich Dichten zu vorgegebenem Erwartungswert und Kovarianzen oder etwa auch Varianzen oft nicht auf einfache Weise bestimmen. F¨ur das genannte An-wendungsgebiet zeigt sich, daß in vielen F¨allen Normalverteilungen ad¨aquat sind, zumal h¨aufig auch wie erw¨ahnt von Gaußschen Verf¨alschungen ausgegangen wird. F¨ur ande-re Anwendungsgebiete gilt dies jedoch im allgemeinen nicht, insbesondeande-re nicht f¨ur eine Verallgemeinerung auf stochastische Prozesse.

Die zweite Strategie ist eine Modifikation des Algorithmus, bei der gerade keine Kova-rianzen (und auch keine VaKova-rianzen) betrachtet werden und an Stelle des bedingten Er-wartungswertes der Modalwert tritt. Die IS–Dichte wird dann auch nicht gem¨aß einer Normalverteilung oder einer anderen vorgegebenen Verteilungsfamilie gew¨ahlt sondern lediglich im Rahmen der Familie der Originaldichte jeweils so verschoben, daß ihr Erwar-tungswert diesem zuvor gesch¨atzten Modalwert entspricht. Die Intuition dahinter ist, daß die IS–Dichte um den Punkt in der Fehlerregion herum zentriert ist, an dem ein Fehler am wahrscheinlichsten ist. Es wird vorausgesetzt, daß der Modalwert eindeutig ist, was f¨ur Normalverteilungen und auch andere Verteilungen zutrifft. Die Sch¨atzung der beding-ten Gr¨oßen wird also vermieden, stattdessen muß hier ein Optimierungsproblem gel¨ost werden. Die verwendete simple Heuristik, die IS–Dichte in Richtung der optimalen IS–

Dichte zu verschieben, entspricht der Translation und ist wesentlich von der Einfachheit der Verschiebung abh¨angig. F¨ur Normalverteilungen entspricht das auch dem exponen-tiellen Maßwechsel. Der Algorithmus realisiert damit also die adaptive Bestimmung von optimalen Translationsparametern und f¨ur Normalverteilungen auch des optimalen expo-nentiellen Maßwechsels. Eine Verallgemeinerung auf komplexere dynamische Systeme und stochastische Prozesse ist aber auch hier kaum m¨oglich.

Es sei abschließend angemerkt, daß f¨ur Normalverteilungen zwar der Erwartungswert und der Modalwert identisch sind, dies f¨ur die bedingten Gr¨oßen jedoch nur asymptotisch gilt.

Somit sind also die von beiden Algorithmen berechneten Verschiebungsparameter auch in diesem Fall nicht notwendigerweise identisch.