Importance Sampling

Importance Sampling (im Deutschen manchmal auch Stichprobenentnahme nach Wichtigkeit, oder Stichprobenziehung nach Wichtigkeit<ref>International Statistical Institute: Glossary of statistical terms.</ref> genannt) ist ein Begriff aus der Statistik, der die Technik zur Erzeugung von Stichproben anhand einer Wahrscheinlichkeitsverteilung beschreibt. Importance Sampling ist eine von mehreren Möglichkeiten zur Varianzreduktion, also zur Steigerung der Effizienz von Monte-Carlo-Simulationen.

Motivation

Monte-Carlo-Simulationen werden oft benutzt, um den Erwartungswert

<math>\operatorname E_X[f(X)] = \begin{cases}\displaystyle \sum_{x \in \Omega} f(x) p(x) & \text{(diskreter Fall)}\\

                                 \displaystyle \int_\Omega f(x) p(x) \mathrm{d} x & \text{(stetiger Fall)} 
                                 \end{cases}</math>

einer reellen Zufallsvariablen <math>f(X)</math> zu berechnen, wobei die Wahrscheinlichkeitsverteilung der Zufallsvariablen <math>X</math> und die Funktion <math>f</math> bekannt sind. Die Zufallsvariable <math>X</math> nimmt Werte in der Ergebnismenge <math>\Omega</math> an. Für eine Realisierung <math>x \in \Omega</math> der Zufallsvariablen <math>X</math> ist <math>f(x) \in \R</math> eine Realisierung der Zufallsvariablen <math>f(X)</math>.

Im diskreten Fall ist <math>\Omega \subset \R^d</math> eine abzählbare Menge und die Wahrscheinlichkeitsverteilung von <math>X</math> ist durch die Wahrscheinlichkeitsfunktion <math>p</math> mit <math>\sum_{x \in \Omega}p(x)=1</math> gegeben. Im stetigen Fall ist <math>\Omega \subseteq \R^d</math> typischerweise ein <math>d</math>-dimensionales Intervall und die Wahrscheinlichkeitsverteilung von <math>X</math> ist durch eine Wahrscheinlichkeitsdichte <math>p</math> mit der Eigenschaft <math>\int_\Omega p(x) \mathrm{d}x =1</math> gegeben. Für eine Realisierung <math>x \in \Omega</math> der Zufallsvariablen <math>X</math> ist <math>f(x) \in \R</math> eine Realisierung der Zufallsvariablen <math>f(X)</math>. Da die Ergebnismenge <math>\Omega</math> im Allgemeinen hochdimensional sein kann, kann die Berechnung des Erwartungswertes sehr schwierig oder zeitaufwendig sein, so dass dann eine Approximation durch Monte-Carlo-Simulation sinnvoll ist.

Bei Anwendungen im Bereich der Physik ist z. B. die Ergebnismenge <math>\Omega</math> der Phasenraum der Teilchen im System, die Zufallsvariablen <math>X</math> und <math>f(X)</math> sind interessierende Größen und <math>p(x)</math> ist proportional zum Boltzmann-Faktor.

Statt den Erwartungswert analytisch zu berechnen oder numerisch zu approximieren, berechnet man im einfachsten Fall den Monte-Carlo-Schätzwert

<math>\operatorname\hat{E}_X[f(X)] = \frac{1}{n} \sum_{i=1}^n f(x_i), </math>

für den Erwartungswert <math>\operatorname E_X[f(X)]</math>, wobei die Zufallszahlen <math>x_1,\dots, x_n</math> Realisierungen von stochastisch unabhängigen und identisch verteilten Zufallsvariablen <math>X_1,\dots, X_n</math> mit der Verteilung von <math>X\sim p</math> sind. In statistischer Terminologie ist <math>(X_1,\dots,X_n)</math> eine einfache Zufallsstichprobe mit Stichprobenumfang <math>n</math>. Der Schätzwert <math>\operatorname\hat E_X[f(X)]</math> ist eine Realisierung des zugehörigen Schätzers

<math>\operatorname\tilde E_X[f(X)] = \frac{1}{n} \sum_{i=1}^n f(X_i)</math>,

der eine Zufallsvariable ist, die nach dem starken Gesetz der großen Zahlen mit Wahrscheinlichkeit Eins gegen den zu schätzenden Erwartungswert konvergiert. Da die Berechnung durch Monte-Carlo-Simulation eine statistische Schätzung ist, gibt es einen Schätzfehler der z. B. durch die Varianz des Schätzer beschrieben werden kann. Beim Importance Sampling wird die zuvor beschriebene einfache Monte-Carlo-Schätzmethode mit dem Ziel modifiziert, die Varianz des Schätzer zu reduzieren bzw. die Genauigkeit bei gegebenem Stichprobenumfang zu erhöhen.

Grundidee des Importance-Sampling

Der Standardansatz zur Approximation des Erwartungswertes <math>\operatorname E[X]</math> einer Zufallsvariablen <math>X</math> durch Monte-Carlo-Simulation besteht darin, Zufallszahlen <math>x_1,\dots,x_n</math> als Realisierungen stochastisch unabhängiger und identisch verteilter Zufallsvariablen <math>X_1,\dots,X_n</math> mit der Wahrscheinlichkeitsverteilung von <math>X</math> zu erzeugen und dann den gesuchten Erwartungswert <math>\operatorname E[X]</math> durch den Monte-Carlo-Schätzwert

<math>\operatorname \hat E_{\mathrm {MC}}[X] = \frac{1}{n} \sum_{i=1}^n x_i</math>

zu approximieren.

Beim Importance-Sampling werden stattdessen Zufallszahlen aus einer modifizierten Wahrscheinlichkeitsverteilung mit dem Ziel verwendet, durch eine Varianzreduktion zu rechentechnisch effizienteren Berechnung zu kommen.

Für eine diskrete Zufallsvariable <math>X</math>, die Werte <math>x \in \Omega</math> mit den Wahrscheinlichkeiten <math>p(x) = P(X=x)</math> annimmt, ist die Grundlage des Importance-Sampling eine Umdeutung des Erwartungswertes

<math> \operatorname E[X]= \sum_{x \in \Omega} xp(x) = \sum_{x \in \Omega} \frac{x p(x)}{w(x)} w(x) =

\operatorname E\left[ \frac{Yp(Y)}{w(Y)} \right]\;.</math> Dabei ist <math>Y</math> eine diskrete Zufallsvariable, die Werte <math>y \in \Omega</math> mit den positiven Wahrscheinlichkeiten <math>w(y)=P(Y=y)</math> annimmt. Die Erzeugung von <math>n</math> Zufällen <math>y_1,\dots,y_n</math> aus der Verteilung von <math>Y</math> führt dann zum Monte-Carlo-Schätzwert

<math> \frac{1}{n} \sum_{i=1}^n \frac{y_i p(y_i)}{w(y_i)} \approx \operatorname E\left[\frac{Yp(Y)}{w(Y)}\right]\;,</math>

der als Importance-Sampling-Schätzwert

<math> \operatorname \hat E_{\mathrm{IS}}[X] = \frac{1}{n} \sum_{i=1}^n \frac{y_i p(y_i)}{w(y_i)}</math>

für <math>\operatorname E[X]</math> interpretiert wird. Der Umweg über die Erzeugung aus Zufallszahlen mit einer anderen Verteilung kann lohnend sein, da durch eine geeignete Wahl der Verteilung von <math>Y</math> ein Importance-Sampling-Schätzer eine erheblich kleinere Varianz als der gewöhnliche Monte-Carlo-Schätzer haben kann.

Importance-Sampling bei unbekannter Normalisierungskonstante

In der Anwendung kommt es häufig vor, dass die Wahrscheinlichkeitsfunktion <math>p</math> nur bis auf eine Konstante bekannt ist: <math>\tilde{p} / Z = p</math>. D.h., es herrscht Kenntnis von <math>\tilde{p}</math> aber nicht von der Normalisierungskonstante <math>Z</math>. Letztere lässt sich ebenfalls durch Importance-Sampling schätzen:<math display="block">Z = \sum_{y\in Y}\tilde{p}(y) = \sum_{y\in Y}\frac{\tilde{p}(y)}{w(y)}w(y) = \operatorname E\left[\frac{\tilde{p}(Y)}{w(Y)}\right] \approx \frac{1}{n}\sum_{i=1}^n \frac{\tilde{p}(y_i)}{w(y_i)}. </math>Ersetzt man <math>Z </math> durch diesen Schätzer in <math>\operatorname \hat E_{\mathrm{IS}}[X] </math>, ergibt sich eine selbstnormalisierte Variante des Importance-Sampling-Schätzers:

<math display="block">\operatorname \hat{E}_{\mathrm{SNIS}}[X] := \frac{\sum_{i=1}^n y_i \tilde{p}(y_i) / w(y_i)}{\sum_{i=1}^n \tilde{p}(y_i) /w(y_i)}. </math>

Im Gegensatz zum klassischen Importance-Sampling-Schätzer ist dieser Schätzer verzerrt, jedoch stark konsistent, d. h. er konvergiert unter geeigneten Momentbedingungen fast sicher gegen den wahren Erwartungswert.

Beispiel

Die Grundidee des Importance-Sampling sei an einem einfachen Beispiel mit <math>\Omega=\{-10,0,10\}</math> veranschaulicht. Die Zufallsvariable <math>X</math> mit der Wahrscheinlichkeitsfunktion

<math>p(x) = P(X = x) = \begin{cases}\frac{1}{10} &\text{für }x \in \{-10,10\}\\

                               \frac{8}{10}&\text{für }x  = 0
                 \end{cases} </math>

hat den Erwartungswert <math>\operatorname E[X]=0</math> und die Varianz <math>\mathrm{Var}[X]=20.</math> Die Wahrscheinlichkeitsfunktion <math>p</math> sei als bekannt vorausgesetzt. Der Erwartungswert sei als unbekannt angenommen und soll durch Simulation bestimmt werden. Der gewöhnliche Monte-Carlo-Schätzer

<math> \operatorname \tilde E_{\mathrm{MC}}[X] = \frac{1}{n} \sum_{i=1}^{n} X_i</math>

basierend auf <math>n</math> Zufallszahlen aus der Verteilung von <math>X</math> hat dann den Erwartungswert Null und die Varianz

<math>\mathrm{Var}[\tilde E_{\mathrm{MC}}[X]] = \frac{\mathrm{Var}[X]}{n} = \frac{20}{n}\;.</math>

Für das Importance-Sampling wird der Erwartungswert <math>\operatorname E[X]</math> als

<math>\operatorname E[X] = \sum_{x \in \Omega} xp(x) = \sum_{x \in \Omega}\frac{xp(x)}{w(x)}{w(x)} = \operatorname E\left[\frac{Yp(Y)}{w(Y)}\right] </math>

geschrieben und die positiven <math>w(x)</math> werden so gewählt, dass sie sich zu Eins addieren und damit als Wahrscheinlichkeiten <math> P(Y= y) = w(y)</math> für eine Zufallsvariable <math>Y</math> interpretiert werden können. Der Importance-Sampling-Schätzer

<math> \operatorname \tilde E_{\mathrm{IS}}[X] = \frac{1}{n} \sum_{i=1}^n \frac{Y_i p(Y_i)}{w(Y_i)}, </math>

der als Funktion der Zufallsvariablen <math>Y_1,\dots,Y_n</math> selbst eine Zufallsvariable ist, hat in diesem Beispiel den Erwartungswert Null und die Varianz

<math> \mathrm{Var}[\operatorname \hat E_{\mathrm{IS}}[X]] = \frac{1}{n} \left(\frac{1}{w(-10)} + \frac{1}{w(10)}\right)</math>

Die Varianz des Schätzers <math> \operatorname \tilde E_{\mathrm{IS}}[X]</math> hängt also von den gewählten Wahrscheinlichkeiten <math>w(y)</math> ab und kann bei geeigneter Wahl erheblich kleiner als die Varianz des gewöhnlichen Monte-Carlo-Schätzers sein. Z. B. ergibt sich für die Wahl <math>w(y)=1/3</math> für <math>y \in\Omega</math> die Varianz

<math>\mathrm{Var}[\operatorname \tilde E_{\mathrm{IS}}[X]] = \frac{6}{n} < \mathrm{Var}[\operatorname \tilde E_{\mathrm{MC}}[X]] = \frac{20}{n} \;. </math>

Durch den Importance-Sampling-Schätzer mit der Verteilung von <math>Y</math> lässt sich also der Erwartungswert <math>\mathbb{E}[X]</math> mit kleinerer Varianz als durch den gewöhnlichen Monte-Carlo-Schätzer bestimmen. Damit erhält man bei gleicher Anzahl <math>n</math> von erzeugten Zufallszahlen eine größere Genauigkeit.

Simple Sampling

Im einfachsten Fall (einfache Stichprobenentnahme, {{#invoke:Vorlage:lang|full|CODE=en|SCRIPTING=Latn|SERVICE=englisch}}) werden aus dem Ergebnisraum gleichverteilt zufällig Zustände <math>Y_i \sim \text{uniform}(\Omega)</math> für die Stichprobe <math>S</math> ausgewählt. Dann ergibt sich für den geschätzten Mittelwert der selbstnormalisierte Importance-Sampling Schätzer:

<math>\operatorname \hat{E}_X[f(X)] = \frac{\sum_{y \in S} p(y) \, f(y)}{\sum_{y \in S} p(y)},</math>

wobei die Summation über die zufälligen Realisierungen <math>y_i \sim Y_i</math>in der Stichprobe läuft. <math>p(y)</math> ist die ursprüngliche Wahrscheinlichkeit(sdichte) für die – durch Simple Sampling erzeugte – Realisierung <math>y</math>. <math>f(y)</math> ist eine Funktionsauswertung

Definition Importance Sampling

Die Methode des Simple Sampling ist meistens nicht sehr effizient, da oft nur wenige relevante Zustände in die Mittelwertbildung eingehen. Um dieses Problem zu umgehen und so die Standardabweichung des gemessenen Mittelwertes bei gleichem Stichprobenumfang zu reduzieren, versucht man Zustände mit einem größeren Gewicht häufiger in die Mittelwertbildung eingehen zu lassen als Zustände mit einem geringeren Gewicht: Der obigen Schätzer des Simple Sampling kann durch Erweitern mit <math>1=W(y)/W(y)</math> auch wie folgt ausgedrückt werden:

<math>\operatorname E_X[f(X)] = \frac{\sum_{y \in S} p(y)/W(y) \, f(y) W(y)}{\sum_{y \in S} p(y)/W(y) W(y)}.</math>

Werden Realisierungen <math>y</math> mit der Wahrscheinlichkeit <math>W(y)</math> erzeugt (Stichprobenentnahme nach Wichtigkeit {{#invoke:Vorlage:lang|full|CODE=en|SCRIPTING=Latn|SERVICE=englisch}}, <math>y \sim W</math>), wird also eine andere Stichprobe S' erzeugt, so berechnet sich der geschätzte Mittelwert in der Folge einfach mithilfe von

<math>\operatorname E_X[f(X)] = \frac{\sum_{y \in S'} p(y)/W(y) f(y)}{\sum_{y \in S'} p(y)/W(y)},</math> für <math>y \sim W</math>

Die Wahrscheinlichkeitsdichte <math>W</math> wird auch als „biased distribution“, „proposal distribution“ oder „sample distribution“ bezeichnet. Da die Stichprobe aus der Verteilung <math>W</math> gezogen wurde, müssen die Erwartungswerte durch entsprechendes Reweighting mit <math>W</math> errechnet werden (siehe Formel) und nicht einfach als arithmetischen Mittel. Dieses Reweighting wird beispielsweise beim „Temperature Reweighting“ von Monte-Carlo Simulationen molekularer Systeme genutzt, bei denen Aussagen über angrenzende Temperaturen gemacht werden sollen.<ref>Bachmann, M. (2014). Thermodynamics and Statistical Mechanics of Macromolecular Systems. Vereinigtes Königreich: Cambridge University Press. Seiten 104, 105 Google Books</ref>

Um eine Stichprobenentnahme nach Wichtigkeit in der Praxis zu erreichen, geht man von einer Startkonfiguration aus und erzeugt mithilfe des Metropolisalgorithmus eine Markow-Kette aus Systemzuständen.

Schlussfolgerungen

Ziehen von Stichproben, ohne dass die Normierungskonstante einer Verteilung berechnet werden kann

Werden die Realisierungen <math>y</math> mit einer Wahrscheinlichkeit <math>W(y)</math> proportional zu <math>p(y)</math> vorgeschlagen (das ist gerade die Metropoliswahl), so ergibt sich (wie beim Gesetz der großen Zahlen)

<math>\operatorname \hat{E}_X[f(X)] = \frac{1}{N}\sum_{y \in S'} f(y).</math>

Gerade, dass hier nur die Proportionalität <math>W(y) \propto p(y)</math> erforderlich ist, ist ein Vorteil der Methode, da die Zustandssumme nicht ausgewertet werden muss.

Simple Sampling

Simple sampling erhält man für den Fall, dass die Vorschlagsdichte konstant ist: <math>W(y)=\text{const}</math>.

Wang-Landau Sampling

Das multikanonische Ensemble kann mit dem Wang-Landau-Algorithmus simuliert werden, indem die Wahl <math>W(x)=1/D(\operatorname E(x))</math> getroffen wird (wobei <math>D(\operatorname E(x))</math> die Zustandsdichte der Energie ist, welche dem Zustand <math>x</math> zugeordnet ist).

Literatur

{{#invoke:Vorlage:Literatur|f}}
{{#invoke:Vorlage:Literatur|f}}
{{#invoke:Vorlage:Literatur|f}}
{{#invoke:Vorlage:Literatur|f}}
{{#invoke:Vorlage:Literatur|f}}
{{#invoke:Vorlage:Literatur|f}}

Einzelnachweise