Numerische Berechnung

Zeit: 10 min

Theorie: Wenn die Mathematik aufgibt

Das Bayes Theorem berechnet die Posterior Verteilung \(P(\theta|D)\). Die Formel erfordert im Nenner die totale Wahrscheinlichkeit der Daten \(D\). Um diesen Wert zu ermitteln, verlangt die Mathematik ein Integral über alle erdenklichen Werte des Parameters \(\theta\).

Bei simplen Modellen ist dieses Integral analytisch lösbar. Die Formel liefert ein exaktes Ergebnis. In der realen Welt bauen wir jedoch oft komplexe Modelle mit zahlreichen Parametern. Hier versagt die reine Mathematik. Das Integral ist schlicht unlösbar.

Die Lösung für dieses Problem lautet Markov Chain Monte Carlo (kurz MCMC). MCMC ist ein numerischer Algorithmus. Er berechnet die Posterior Verteilung nicht analytisch exakt, sondern approximiert sie durch intelligente Simulation.

Der Algorithmus wandert schrittweise durch den Parameterraum. Er verbringt automatisch viel Zeit an Orten mit hoher Wahrscheinlichkeit und wenig Zeit an Orten mit niedriger Wahrscheinlichkeit. Sammelt man all diese besuchten Orte am Ende auf, erhält man ein perfektes Abbild der gesuchten Posterior Verteilung \(P(\theta|D)\).

Beispiel SmartRail: Der blinde Optimierer

Du willst SmartRails wahre durchschnittliche Verspätungsreduktion \(\theta\) schätzen, basierend auf den gesammelten Testdaten \(D\). Dein Modell berücksichtigt jedoch zahlreiche Faktoren wie Tageszeit, Streckenauslastung und Wetter. Die Formel für den Posterior \(P(\theta|D)\) ist für deinen Computer zu komplex zum direkten Ausrechnen.

Stelle dir den MCMC Algorithmus wie einen blinden SmartRail-Analysten vor, der auf einem Zahlenstrahl möglicher Verspätungsreduktionen wandert. Je besser ein Wert auf dem Zahlenstrahl zu den echten Testdaten \(D\) passt, desto höher ist die Wahrscheinlichkeit. Der Analyst will den wahrscheinlichsten Bereich kartografieren.

Er startet blind bei einem Startwert, etwa \(\theta = 0\) Minuten Reduktion. Er macht einen zufälligen Schritt zu einer neuen Position, etwa \(\theta = 3\) Minuten. Er prüft die neue Position. Da 3 Minuten Reduktion deutlich besser zu den empirischen Testdaten \(D\) passen als 0 Minuten, akzeptiert er den Schritt. Das ist ein neuer Eintrag in seiner Kette.

Manchmal geht er bewusst in einen Bereich geringerer Wahrscheinlichkeit. Das verhindert, dass er sich in einem lokalen Maximum verrennt und das wahre globale Zentrum verfehlt. Lässt du diesen blinden Analysten zehntausend Schritte machen, zeigt dir die Dichte seiner gesammelten Fußspuren exakt, wo SmartRails wahre durchschnittliche Verspätungsreduktion \(\theta\) am wahrscheinlichsten liegt.

Deine Aufgabe

Die Applikation simuliert diesen blinden SmartRail-Analysten. Die rote Zielkurve repräsentiert den wahren, mathematisch unlösbaren Posterior \(P(\theta|D)\) für die Verspätungsreduktion \(\theta\).

  1. Beobachten: Klicke bei wenigen Schritten auf den Button. Der Analyst irrt kurz umher. Das untere Diagramm (seine gesammelten Fußspuren) sieht der roten Kurve noch nicht ähnlich.
  2. Iterieren: Erhöhe die Anzahl der Schritte massiv. Lass den Analysten den Parameterraum ausgiebig erkunden.
  3. Ergebnis analysieren: Das graue Histogramm der besuchten Orte nähert sich sehr gut der roten Zielkurve an. Du hast das unlösbare mathematische Integral durch reine Rechenpower und Simulation geknackt.

⏳ Die Anwendung wird geladen, dies kann bis zu 30 Sekunden dauern.

#| '!! shinylive warning !!': |
#|   shinylive does not work in self-contained HTML documents.
#|   Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 650

from shiny import App, render, ui, reactive
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

app_ui = ui.page_fluid(
    ui.card(
        ui.card_header("SmartRail: MCMC findet die wahre Verspätungsreduktion"),
        ui.layout_columns(
            ui.div(
                ui.h5("Steuerung des Analysten"),
                ui.input_slider("n_steps", "Anzahl der Schritte", 100, 5000, 200, step=100),
                ui.input_action_button("run_btn", "Wanderung starten", class_="btn-primary"),
                ui.hr(),
                ui.p("Der Analyst sucht SmartRails wahrscheinlichste Verspätungsreduktion Theta basierend auf den Testdaten D."),
            ),
            ui.output_plot("mcmc_plot"),
            col_widths=(4, 8)
        )
    )
)

def server(input, output, session):
    
    sim_data = reactive.Value(None)
    
    @reactive.Effect
    @reactive.event(input.run_btn, ignore_none=False)
    def run_mcmc():
        steps = input.n_steps()
        
        target_mu = 6.0
        target_sd = 2.5
        
        samples = np.zeros(steps)
        current_theta = 0.0 
        
        np.random.seed(input.run_btn() if input.run_btn() else 42)
        
        for i in range(steps):
            proposal_theta = current_theta + np.random.normal(0, 2.0)
            
            prob_current = norm.pdf(current_theta, target_mu, target_sd)
            prob_proposal = norm.pdf(proposal_theta, target_mu, target_sd)
            
            acceptance_ratio = prob_proposal / prob_current if prob_current > 0 else 1.0
            
            if np.random.uniform(0, 1) < acceptance_ratio:
                current_theta = proposal_theta
                
            samples[i] = current_theta
            
        sim_data.set(samples)

    @render.plot
    def mcmc_plot():
        samples = sim_data.get()
        if samples is None:
            return
            
        fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8), layout="constrained")
        
        ax1.plot(samples, color='#34495e', alpha=0.7, linewidth=0.5)
        ax1.set_title("Wanderweg der MCMC Kette (Trace Plot)")
        ax1.set_xlim(0, len(samples))
        ax1.set_ylim(-3, 15)
        ax1.set_xlabel("Schrittnummer")
        ax1.set_ylabel("Besuchter Wert Theta")
        ax1.spines['top'].set_visible(False)
        ax1.spines['right'].set_visible(False)
        
        x_vals = np.linspace(-5, 18, 200)
        target_y = norm.pdf(x_vals, 6.0, 2.5)
        
        ax2.hist(samples, bins=50, density=True, color='#95a5a6', alpha=0.7, label='Besuchte Orte (Simulation)')
        ax2.plot(x_vals, target_y, color='#e74c3c', linewidth=3, label='Wahrer Posterior P(Theta|D)')
        
        ax2.set_title("Vergleich Simulation und echtes Ziel")
        ax2.set_xlim(-5, 18)
        ax2.set_yticks([])
        ax2.set_xlabel("Durchschnittliche Verspätungsreduktion Theta (Minuten)")
        ax2.set_ylabel("Wahrscheinlichkeitsdichte")
        ax2.legend(loc="upper right")
        ax2.spines['top'].set_visible(False)
        ax2.spines['right'].set_visible(False)
        
        return fig

app = App(app_ui, server)