MCMC-Verfahren

Markow-Chain-Monte-Carlo-Verfahren (kurz MCMC-Verfahren; seltener a​uch Markow-Ketten-Monte-Carlo-Verfahren) s​ind eine Klasse v​on Algorithmen, d​ie zufällige Stichproben a​us Wahrscheinlichkeitsverteilungen (Monte-Carlo-Verfahren) ziehen. Dies geschieht a​uf der Basis d​er Konstruktion e​iner Markow-Kette, welche d​ie erwünschte Verteilung a​ls ihre stationäre Verteilung aufweist. Der Zustand d​er Kette n​ach einer großen Zahl v​on Schritten w​ird dann a​ls Stichprobe d​er erwünschten Verteilung benutzt. Die Qualität d​er Stichprobe steigt m​it zunehmender Zahl d​er Schritte.

Konvergenz des Metropolis-Hastings-Algorithmus; die blaue Verteilung soll durch die orange approximiert werden.

Üblicherweise gelingt es leicht, eine Markow-Kette mit den erwünschten Eigenschaften zu konstruieren. Das Schwierigere ist, zu ermitteln, wie viele Schritte nötig sind, um Konvergenz zur stationären Verteilung mit akzeptablem Fehler zu erreichen. Also den Algorithmus so zu gestalten, dass möglichst effektiv unabhängige Systemzustände generiert werden. Eine gute Kette mit einer gut gewählten Anfangsverteilung wird schnell konvergieren, d. h. die stationäre Verteilung wird schnell erreicht. Bei typischer Anwendung von MCMC-Verfahren kann die Zielverteilung nur näherungsweise erreicht werden, da es immer einen gewissen Resteffekt der Anfangsverteilung gibt. Stichproben, welche mit MCMC-Algorithmen generiert werden, weisen typischerweise hohe Autokorrelationen auf. Daher wird der Fehler des Mittelwerteschätzers bei Verwendung des Standardfehlers unterschätzt.

Anwendungsgebiete

Häufige Anwendungen dieser Algorithmen finden s​ich bei d​er numerischen Berechnung mehrdimensionaler Integrale. Diese finden s​ich oft i​m Rahmen d​er bayesschen Statistik[1] s​owie in Computersimulationen i​n der statistischen Physik o​der zur Auswertung v​on Pfadintegralen i​n Quantenfeldtheorien[2]

MCMC-Verfahren können z. B. z​ur Simulation v​on Systemen i​m kanonischen Ensemble herangezogen werden. Simulationen n​ahe einem Phasenübergang konvergieren i. d. R. langsamer u​nd erfordern e​ine längere Equilibrierung. Eine hinreichende, a​ber nicht notwendige, Bedingung, d​ass ein MCMC-Verfahren d​en kanonischen Zustand a​ls stationäre Verteilung aufweist, i​st die Detailed-Balance-Eigenschaft.

Beispiel

Angenommen, wir können eine diskrete Gleichverteilung auf ${\displaystyle \{0,1\}}$ simulieren. Dies entspricht dem Wurf einer fairen Münze mit ${\displaystyle P(\{1\})=P(\{0\})=0{,}5}$, wobei "1=Kopf" und "0=Zahl" sein soll. Nun wollen wir aber eine Verteilung auf ${\displaystyle S=\{1,2,3\}}$ simulieren mit ${\displaystyle P_{Z}(\{1\})=0{,}2}$ und ${\displaystyle P_{Z}(\{2\})=P_{Z}(\{3\})=0{,}4}$. Dazu definieren wir eine Markow-Kette, bei der alle Übergangswahrscheinlichkeiten simuliert werden können (sprich 0,5, 0 oder 1 sind) und die die Verteilung ${\displaystyle P_{Z}}$ als stationäre Verteilung hat. Folgendes Vorgehen bildet eine Markow-Kette, welche die Voraussetzungen erfüllt: Befindet man sich im Zustand 1, gehe nach 2. Befindet man sich im Zustand 2, wirf eine Münze. Zeigt sie Kopf, gehe zu 3, ansonsten verharre in der 2. Befindet man sich im Zustand 3, wirf eine Münze. Zeigt sie Kopf, gehe zu 1, ansonsten verharre in der 3. Diese Markow-Kette wird durch die folgende Übergangsmatrix beschrieben:

${\displaystyle M={\begin{bmatrix}0&1&0\\0&0{,}5&0{,}5\\0{,}5&0&0{,}5\end{bmatrix}}=((P(i|k))_{k,i})}$.

Da alle Einträge der Matrix ${\displaystyle M^{t}>0}$ positiv sind ab ${\displaystyle t\geq 3}$, ist die Markow-Kette irreduzibel und aperiodisch. Daher konvergiert die Markow-Kette für jede beliebige Startverteilung ${\displaystyle v_{0}}$ gegen die stationäre Verteilung ${\displaystyle \pi =P_{Z}}$

Für die Simulation starten wir nun im Zustand 1. In der Tabelle findet sich jeweils der Simulationsschritt in der ersten Spalte, dann die exakte Aufenthaltswahrscheinlichkeit für den Zustand 1 in der 2. Spalte: ${\displaystyle P(X_{k}=1)=\sum _{i=1}^{3}P(1|i)P(X_{k-1}=i)}$

Diese konvergiert nach Voraussetzung gegen 0,2. In der letzten Spalte steht die relative Häufigkeit des Zustandes 1, gemittelt über eintausend Simulationen. Diese nähert sich nach mehreren Simulationsschritten der stationären Verteilung an. Lässt man also die Simulation lange laufen, so ist die so erhaltene Stichprobe ${\displaystyle P_{Z}}$ verteilt.

Schritte ${\displaystyle k}$	${\displaystyle P(X_{k}=1)}$	Relative Häufigkeit der 1 in Tausend Simulationen
0	1	1
1	0	0
2	0	0
3	0,25	0,228
4	0,25	0,271
5	0,1875	0,173
6	0,1875	0,176
7	0,2031	0,236
8	0,2031	0,180
9	0,1992	0,202
10	0,1992	0,171
11	0,2002	0,205
12	0,2002	0,198
13	0,2000	0,190
14	0,2000	0,206

Normiert man die Summe der Einträge des linksseitigen Eigenvektors von ${\displaystyle M}$ zum Eigenwert 1, so erhält man die Wahrscheinlichkeiten der Zustände der stationären Wahrscheinlichkeitsverteilung der Markow-Kette: hier 0,2, 0,4, 0,4.

Algorithmen

Beispiele für Markow-Chain-Monte-Carlo-Verfahren sind:

• Metropolisalgorithmus: Das lokale Verfahren erzeugt Zufallsbewegungen, die mit einer gewissen Wahrscheinlichkeit akzeptiert oder zurückgewiesen werden. Es ist einfach zu implementieren, hat jedoch den Nachteil einer hohen Autokorrelation der erzeugten Systemzustände.
• Gibbs-Sampling: Das lokale Verfahren ist ein Sonderfall des Metropolis-Hastings-Algorithmus, bei dem die Zustände entsprechend der lokalen Wahrscheinlichkeitsverteilung erzeugt werden.
• Wärmebadalgorithmus
• Hybrid-Monte-Carlo-Algorithmus: Das Verfahren stellt eine Kombination aus Molekulardynamik und Zufallsbewegung her. Die Molekulardynamik wird benutzt, um effizient neue, unabhängige Zustände zu erzeugen, die mit einer gewissen Wahrscheinlichkeit akzeptiert oder zurückgewiesen werden.
• Clusteralgorithmen: Hierbei handelt es sich um spezielle, nicht-lokale Verfahren, die die Autokorrelationzeiten und damit das so genannte Critical Slowing down, bei Phasenübergängen, reduzieren. Ein solcher Algorithmus wurde erstmals 1987 von R. Swendsen und J.-S.Wang für ein Potts-Modell benutzt, siehe Swendsen-Wang-Algorithmus. Eine populäre Variante dieser Clustermethode ist der 1989 eingeführte Wolff-Algorithmus.
• Over-Relaxation-Verfahren

Literatur

• Jeff Gill: Bayesian methods: a social and behavioral sciences approach. Second Edition Auflage. Chapman and Hall/CRC, London 2008, ISBN 978-1-58488-562-7.
• N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller und E. Teller: Equation of State Calculations by Fast Computing Machines. In: Journal of Chemical Physics. Band 21, 1953, S. 1087–1092, doi:10.1063/1.1699114.
• W. K. Hastings: Monte Carlo Sampling Methods Using Markov Chains and Their Applications. In: Biometrika. Band 57, 1970, S. 97–109.
• R. H. Swendsen, J.-S. Wang: Nonuniversal critical dynamics in Monte Carlo simulations. In: Phys. Rev. Lett. Band 58, 1987, S. 86.
• S. L. Adler: Overrelaxation method for Monte Carlo evaluation of the partition function for multiquadratic actions. In: Phys. Rev. D. Band 23, 1988, S. 2901.
• Thomas Müller-Gronbach, Erich Novak, Klaus Ritter: Monte Carlo-Algorithmen. Springer-Verlag, Berlin 2012, ISBN 978-3-540-89140-6, Kapitel 6, S. 179–244.
• Charles Geyer: Introduction to Markov Chain Monte Carlo. In: Chapman & Hall/CRC Handbooks of Modern Statistical Methods. Chapman and Hall/CRC, 2011, ISBN 978-1-4200-7941-8, doi:10.1201/b10905-2 (mcmchandbook.net [PDF]).

Einzelnachweise

1. Sudipto Banerjee, Bradley P. Carlin, Alan P. Gelfand: Hierarchical Modeling and Analysis for Spatial Data, Second. Auflage, CRC Press, 12. September 2014, ISBN 978-1-4398-1917-3, S. xix.
2. Ankur Gupta, James B. Rawlings: Comparison of Parameter Estimation Methods in Stochastic Chemical Kinetic Models: Examples in Systems Biology. In: AIChE Journal. 60, Nr. 4, April 2014, S. 1253–1268. doi:10.1002/aic.14409. PMID 27429455. PMC 4946376 (freier Volltext).