Monte-Carlo Methode - Grundlagen Der Messunsicherheit Für ChemikerInnen

Die Monte-Carlo Methode erlaubt uns, die komninierte Messunsicherheit von Größen mit nicht-normalverteilten Messunsicherheiten zu ermitteln. Dabei gehen wir wie folgt vor:

wir spezifizieren die Verteilungen der Messunsicherheiten aller Größen, die in unsere Messfunktion eingehen
wir berechnen wiederholt den Wert der Messfunktion, wobei wir jeweils Zufallswerte aus den Verteilungen der Messunsicherheiten sampeln
die Verteilung (Histogramm) der berechneten Messfunktionen nähert sich mit der Zeit immer mehr der Verteilung der kombinierten Messunsicherheit an

(In den Beispielen im Grundpraktikum werden diese Methoden im Allgemeinen nicht benötigt.)

import matplotlib.pyplot as plt
import numpy as np
from IPython.display import HTML, Markdown, display
from matplotlib.animation import FuncAnimation

Beispiel: Messkolben¶

In unserem Messkolben Beispiel haben wir eine dreiecksförmige Messunsicherheit durch eine Normalverteilung angenähert. Es stellt sich die Frage, ob diese Annäherung das Ergebnis beeinflusst. Wir verwenden eine Monte Carlo Methode um das zu überprüfen.

Zuerst definieren wir die Verteilungen der Messgrößen:

from scipy.stats import norm, triang

V_dist = triang(loc=50 - 0.05, scale=2 * 0.05, c=0.5)
m_dist = norm(loc=248.6, scale=3.1)

Dann die Messfunktion:

def meas_function(m, V):
    return m / V

Schließlich berechnen wir die Messfunktion mit Samples aus der Verteilung der Messunsicherheiten. N gibt hierbei die Anzahl der zufälligen Samples an. Viele Samples (> 10⁴) werden benötigt, eine gute Annäherung der Verteilung der kombiniierten Messunsicherheit zu erreichen:

N = int(1e5)
beta = meas_function(m_dist.rvs(size=N), V_dist.rvs(size=N))

Der Erwartungswert und Messunsicherheit werden nun als Mittelwert und Standardabweichung der zufälligen Messwerte ermittelt:

display(Markdown(rf"{np.mean(beta):.2f} mg/mL ± {np.std(beta):.2f} mg/mL"))

Animation Monte Carlo¶

Die zunehmende Annäherung der Funktion an die Verteilung der kombinierten Messunsicherheit ist der in der Animation weiter unten dargestellt:

%%capture
from numpy import histogram

fig, ax = plt.subplots()

x = np.linspace(4.8, 5.2, 1000)
ax.plot(x, norm(loc=4.97, scale=0.062).pdf(x), label="aus der Fehlerfortpflanzung")
max_norm = norm(loc=5, scale=0.06).pdf(10)
lin = ax.plot(x, x * 0, label="aus Monte-Carlo")[0]
text = ax.annotate(
    "N = 0", (0.05, 0.8), ha="left", xycoords="axes fraction", size="large"
)


m_smps = m_dist.rvs(size=int(1e8))
V_smps = V_dist.rvs(size=int(1e8))


def update_lin(N):
    N = int(N)
    samples = meas_function(m_smps[:N], V_smps[:N])
    hist, bins = histogram(samples, bins=min(N, 500), range=(4.8, 5.2), density=True)
    y = hist
    x = (bins[:-1] + bins[1:]) / 2
    lin.set_xdata(x)
    lin.set_ydata(y)
    text.set_text(f"N = {N}")
    return fig, lin, text


ax.set_xlabel(r"$\beta$ / g/L")
ax.legend()


anim = FuncAnimation(
    fig,
    update_lin,
    frames=10 ** np.linspace(1, 8, 24 * 8),
    interval=1000 / 24,
    blit=True,
)