Übungsblatt 4: MCMC

Aufgabe 1
Aufgabe 2



In [1]:

    
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
from functools import partial

plt.style.use('ggplot')

Aufgabe 1

Erzeugen Sie Paare von Zufallszahlen (x, y) die gemäß der zweidimensionalen Dichte

\begin{equation} f(x, y) \propto \exp\left(-((x-2)x + y^2)/10\right)\cdot\sin^2(x\cdot y + x) \end{equation}

verteilt sind. Bestimmen Sie daraus die Verteilung $\rho(r)$ mit $r^2 = x^2 + y^2$.



In [2]:

    
def mcmc_step(pdf, xs):
    ds = np.random.uniform(-3, 3, size=len(xs))
    r = np.random.uniform()
    nxs = xs + ds
    p0 = pdf(xs)
    p1 = pdf(nxs)
    if p1 > p0 or r < (p1 / p0):
        return nxs
    return xs


def mcmc_sample(pdf, xs, size=100000):
    xss = []
    for _ in range(size):
        xs = mcmc_step(pdf, xs)
        xss.append(xs)
    return np.array(xss).T


def mcmc_generate(pdf, xs, size=100000):
    ds = np.random.uniform(-3, 3, size=(size, len(xs)))
    r = np.random.uniform(size=size)
    p0 = pdf(xs)
    xss = []
    for i in range(size):
        nxs = xs + ds[i]
        p1 = pdf(nxs)
        if p1 > p0 or r[i] < (p1 / p0):
            xs = nxs
            p0 = p1
        xss.append(xs)
    return np.array(xss).T


def pdf(xs):
    x, y = xs
    return np.exp(-((x - 2) * x + y ** 2) / 10) * np.sin(x * y + x) ** 2


%time xs, ys = mcmc_sample(pdf, np.array([0, 0]))
%time xs, ys = mcmc_generate(pdf, np.array([0, 0]))









    



CPU times: user 1.34 s, sys: 33.3 ms, total: 1.37 s
Wall time: 1.41 s
CPU times: user 633 ms, sys: 18.3 ms, total: 652 ms
Wall time: 687 ms



In [3]:

    
plt.hist2d(xs, ys, bins=40)
plt.show()

plt.hist(np.sqrt(xs ** 2 + ys ** 2), bins=40)
plt.show()

Aufgabe 2

Ein Experiment beobachtet $n$ Ereignisse. Die Nachweiswahrscheinlichkeit für ein Ereignis sei $ε$. Die wahre Anzahl ist Poisson-verteilt mit einem Mittelwert $μ$, der eine Funktion zweier Theorieparameter $a$ und $b$ ist, $μ = a^b$. Aus Kalibrationsmessungen ist bekannt dass $ε$ Gauss-verteilt ist mit $ε = 0.75 \pm 0.05$. Vom Parameter $b$ ist bekannt dass er nahe bei $b = 1$ liegt. Für ihn wird eine gleichförmige Verteilung im Bereich $0.9 < b < 1.1$ angenommen. Der Parameter $a$ muss positiv sein. Verwenden Sie bayes’sche Statistik und MCMC, und bestimmen Sie für $n = 1$ und $n = 10$ die posterior Verteilung von $a$ nach Integration über die Nuissance-Parameter $ε$ und $b$.

Zunächst müssen wir berechnen, wie die Effizienz die Verteilung der $n$ beeinflusst. Es sei $q_n$ die Wahrscheinlichkeit $n$ Ereignisse zu beobachten, wenn die $n$ um einen Mittelwert $\mu$ Poisson-verteilt sind und mit einer Effizienz $\varepsilon$ beobachtet werden. Es kann also sein, dass genau $n$ Ereignisse auftreten und alle beobachtet werden, oder $n+1$ beobachtet werden und eins nicht usw. Das bedeutet

\begin{align} q_n &= \sum_{k=n}^\infty \mathrm{Poisson}(k, \mu) \binom{k}{n} \varepsilon^n (1-\varepsilon)^{k-n} \\ &= \sum_{k=n}^\infty \exp(-\mu) \frac{\mu^k}{k!} \frac{k!}{n!(k-n)!} \varepsilon^n (1-\varepsilon)^{k-n} \\ &= \exp(-\mu) \frac{\mu^n}{n!} \sum_{k=n}^\infty \frac{\mu^k (1-\varepsilon)^{k-n}}{(k-n)!} \\ &= \exp(-\mu) \frac{(\varepsilon\mu)^n}{n!} \sum_{k=n}^\infty \frac{(\mu(1-\varepsilon)^{k-n}}{(k-n)!} \\ &= \exp(-\mu) \frac{(\varepsilon\mu)^n}{n!} \sum_{k=0}^\infty \frac{(\mu(1-\varepsilon)^{k}}{(k)!} \\ &= \exp(-\mu) \frac{(\varepsilon\mu)^n}{n!} \exp(\mu(1-\varepsilon)) \\ &= \exp(-\varepsilon\mu) \frac{(\varepsilon\mu)^n}{n!} \\ &= \mathrm{Poisson}(n, \varepsilon\mu) \\ \end{align}

Damit ist

\begin{equation} p(n\,|\,a,b,\varepsilon) \propto \exp\left(\mu\varepsilon\right)\cdot\left(\mu\varepsilon\right)^n \quad\text{mit}\quad \mu = a^b \,. \end{equation}

Mit Hilfe von Bayes Theorem können wir nun die Posteriorverteilung der unbekannten Parameter $\left\{a,b,\varepsilon\right\}$ ausdrücken. Es ist

\begin{equation} p(a, b, \varepsilon\,|\,n) \propto p(n\,|\,a, b, \varepsilon) \cdot p(a, b, \varepsilon) \end{equation}

mit der Priorwahrscheinlichkeitsdichte

\begin{equation} p(a, b, \varepsilon) \propto \Theta(a) \cdot \Theta(b - 0.9) \cdot \Theta(1.1 - b) \cdot \exp\left(-200(\varepsilon - 0.75)^2\right) \,. \end{equation}

Dieser Ausdruck wird jetzt per MCMC über $\varepsilon$ und $b$ integriert um die Verteilung von $p(a\,|\,n)$ zu erhalten.



In [4]:

    
uniform_b = scipy.stats.uniform(0.9, 1.1).pdf
gaussian_eps = scipy.stats.norm(0.75, 0.05).pdf

def p(a, n, eps, b):
    if a <= 0:
        return 0
    epsab = eps * a**b
    return scipy.stats.poisson.pmf(n, epsab) * uniform_b(b) * gaussian_eps(eps)

Diese können wir dann bequem mit unseren, in Aufgabe 1 formulierten, MCMC-Funktionen verwenden. Wir müssen uns nur einen kleinen Wrapper schreiben.



In [ ]:

    
def p_n(n, pars):
    a, eps, b = pars
    return p(a, n, eps, b)


p_1 = partial(p_n, 1)
as1, _, _ = mcmc_generate(p_1, np.array([0.01, 0.75, 1]), size=1000000)

plt.hist(as1, bins=40)
plt.title('$n = 1$')
plt.xlabel('$a$')
plt.ylabel('Absolute Häufigkeit')
plt.show()



In [ ]:

    
p_10 = partial(p_n, 10)
as10, _, _ = mcmc_generate(p_10, np.array([0.01, 0.75, 1]), size=1000000)

plt.hist(as10, bins=40)
plt.title('$n = 10$')
plt.xlabel('$a$')
plt.ylabel('Absolute Häufigkeit')
plt.show()