JGU Logo JGU Logo JGU Logo JGU Logo

Institut für Informatik

Modellierung 2

Michael Wand
David Hartmann
Sommersemester 2022

Lehreinheit 3

Bayes'sche Statistik & etwas statistische Datenanalyse
Letzte Änderung: 02. May 2022, 11:00 Uhr
Abgabe: Montag, der 16.05.2021, 10 Uhr

  /  




Augabe 1: Bayes'sche Inferenz

Angenommen, wir werfen eine Münze 50 mal. Sie zeigt 20 mal "Kopf" und 30 mal "Zahl".
Die Frage die wir uns jetzt stellen ist, ob wir es mit einer fairen Münze zu tun haben.
Hierzu probieren wir drei Möglichkeiten vorzugehen, um herauszufinden, ob wir es mit einer fairen Münze zu tun gehabt haben.


  1. Frequentistische Ansicht (alternative Hypothesen):
    Es gibt gute Münzen und betrügerische Münzen. Gute Münzen haben eine 50/50 Chance für Kopf und Zahl. Eine betrügerische zeige als Beispiels ein mittleres Verhältnis von 75 Kopf zu 25 Zahl (statt 50/50).
    Berechnen Sie die Wahrscheinlichkeit für die beiden Alternativen (Betrug & gut).
  2. Frequentistische Ansicht (Signifikanztest):
    Es gibt gute Münzen und betrügerische Münzen. Gute Münzen haben eine 50/50 Chance für Kopf und Zahl.
    Wie hoch ist die Wahrscheinlichkeit, dass die Münze nicht verfälscht wurde und trotzdem mehr als 10% (5 Würfe) vom Mittelwert (25/25) abweicht?
  3. Bayes'sche Sichtweise:
    Die Wahrscheinlichkeit, dass unsere Münze "Kopf" zeigt, beträgt \(\theta \in [0,1]\) . Berechnen Sie eine Wahrscheinlichkeitsverteilung für den Wert von \(\theta\) angesichts der Daten, die Sie gesehen haben (nehmen Sie einen gleichmäßigen Prior auf \(\theta\) an, d.h., bevor Sie die Daten gesehen haben, werden alle Werte von \(\theta\) als gleich wahrscheinlich angesehen).
    Zeichnen Sie die Verteilung in ein Diagramm.

Aufgabe 2: Welford-Algorithmus


Welford-Algorithmus

Gegeben Sei \(X\) eine 1-d Zufallsvariable und seien weiter \(x_1, \dots, x_n\) beobachtete Werte dieser Zufallsvariable.
Dann ist \( \overline x_k := \frac{1}{k} \sum_i x_i \) der Mittelwert und \( \sigma_k := \frac{1}{k} \sum_i \left(x_i - \overline x_k\right)^2 \) die Varianz des (Teil-)Datensatzes \((x_1,\dots,x_k)\).


Erwartungswert:
Die folgende Umformung führt zu der Online-Variante dieser Berechnung: \[ \overline x_n := \overline x_{n-1} + \frac{x_n - \overline x_{n-1}}{n}. \]


Varianz:
Ganz ähnlich lässt sich auch die online-Berechnung der Varianz herleiten: \[ \sigma_n^2 = \frac{(n-1)\sigma_{n-1}^2+(x_n - \overline x_{n-1})(x_n-\overline x_n)}{n} = \sigma^2_{n-1} + \frac{(x_n - \overline x_{n-1})(x_n - \overline x_n) - \sigma^2_{n-1}}{n}. \] Diese Art die Varianz zu bestimmen kann jedoch zu numerischen Problemen führen — wir subtrahieren hier eine potentiell sehr kleine Zahl von einer beliebig großen Zahl. Eine numerisch stabilere Variante nutzt den linken Teil der Gleichung bzw. die Tatsache, dass die Formel der Standardabweichung die Summe aller Differenzen vom Mittelwert enthält: \[ \begin{align} M_{2,n} &= M_{2,n-1} + (x_n - \overline x_{n-1})(x_n - \overline x_n)\\ \sigma^2_n &= \frac{M_{2,n}}{n}. \end{align} \]


Aufgabe 2a

Gegeben Sei ein zufälliger 1-d Datensatz. Der Einfachheit halber nehmen wir hier eine normalverteilte Zufallsvariable \(X \sim \mathcal{N}_{0,1}\) an.

  1. Implementieren den obigen Welford Algorithmus zur Bestimmung der Varianz in PyTorch.
    Testen Sie an einem einfachen Beispiel, dass die Berechnung korrekt ist.
  2. Wir leiten im folgenden analog zum Welfordalgorithmus einen Algorithmus her, der zusätzlich auch die nächsten zwei statistischen Momente laufend bestimmt, also die Schiefe \(\mu_3/\sigma^3\) und die Wölbung \(\mu_4/\sigma^4\). Wir definieren diese als \[ \frac{\mu_3}{\sigma^3} = E\left[\left(\frac{X - \mu}{\sigma}\right)^3 \right],\qquad \frac{\mu_4}{\sigma^4} = E\left[\left(\frac{X - \mu}{\sigma}\right)^4 \right]. \] Ganz ähnlich zur Herleitung im Fall der Varianz, kann eine Hilfvariable für alle weiteren Momente rekursiv angegeben werden:
    Sei also \(p > 1\) eine natürliche Zahl. Dann gilt für \(M_{p,n} = \sum_i^n (x - \mu)^p\) \[ M_{p,n} = M_{p,n-1} + \sum_{k=1}^{p-2} \binom{p}{k} M_{p-k,n-1} \left(-\frac{x_n - \mu_{n-1}}{n}\right)^k + \left(\frac{(n-1) \cdot (x_n - \mu_{n-1})}{n}\right)^p\left[1 - \left(\frac{-1}{n-1} \right)^{p-1} \right] \] Optional: Zeigen Sie, dass die Gleichung korrekt ist.
    (Mit dem Wissen wie die Formel aussieht ist die Herleitung weder sehr schwer noch lang, da aber die folgenden Aufgaben wichtiger sind, kann diese Aussage gerne als wahr angenommen werden).

    Leiten Sie nun eine Online-Formel zur Berechnung der Wölbung und der Schiefe her und vervollständigen Sie ihre Implementierung.
  3. Generieren Sie einen normalverteilten Datensatz & messen Sie zuerst nach, dass Mittelwert und Standardabweichung mit den Parametern der gegebenen CDF übereinstimmen.
    Verändern Sie die Zufallsvariable so, dass Sie die vier statistischen Momente Erwartungswert, Varianz, Schiefe bzw. die Wölbung seperat steuern können.
  4. Zuletzt noch etwas anderes praktisch relevantes: Es ist stets möglich sich die letzten \(m\) Beispiele zu merken und so eine lokale Statistik zu berechnen. Durch eine kleine Änderung in der obigen Formel des Mittelwertes ist es möglich dieses ebenfalls online zu formulieren. \[ \operatorname{EMA}_\alpha(\overline x_n) := (1-\alpha)\cdot \operatorname{EMA}_\alpha(\overline x_{n-1}) + \alpha \cdot x_n. \] Bis auf den Faktor \((1-\alpha)\) ist diese Formel identisch mit der oberen Version. Diese Änderung lässt jedoch vergangene Werte langsam verfallen.

    Wie muss \(\alpha\) gewählt werden, dass der so definierte lokale Mittelwert zu 99.9% aus den letzten \(100\) Beispielen besteht?

Aufgabe 2b


  1. Bestimmen Sie so den Mittelwert und die Varianz des MNIST-Standarddatensatzes. Der Datensatz wird in Pytorch im Grunde bereits mitgeliefert. (Sie benötigen lediglich zusätzlich das Python-Paket torchvison.)

    Der folgende Code lädt den Datensatz herunter (~50 Mb) und zeigt einige Bilder in Matplotlib an. (Quelle: https://pytorch.org/tutorials/beginner/blitz/cifar10_tutorial.html) 
    import matplotlib.pyplot as plt
import numpy as np

# functions to show an image


def imshow(img):
    img = img / 2 + 0.5     # unnormalize
    npimg = img.numpy()
    plt.imshow(np.transpose(npimg, (1, 2, 0)))
    plt.show()


batch_size = 64 # number of images to get for every iteration of trainloader
trainset = torchvision.datasets.MNIST(root='./data', train=True, download=True, transform=transform)
trainloader = torch.utils.data.DataLoader(trainset, batch_size=batch_size, shuffle=True, num_workers=2)
classes = ('class-%i' % i for i in range(10))

# get some random training images
dataiter = iter(trainloader)
images, labels = dataiter.next()

# show images
imshow(torchvision.utils.make_grid(images))
# print labels
print(' '.join('%5s' % classes[labels[j]] for j in range(batch_size)))
  2. Bestimmen Sie nun die Mittelwerte der jeweiligen Klassen \(0,1,2,3,4,5,6,7,8\) und \(9\) aus diesem Datensatz.
  3. Können Sie dem Mittelwerten entnehmen, wie Sie prüfen können ob ein ungesehenes Bild (dazu gibt es den Test-Datensatz) zu einer bestimmten Klasse gehört?

    Hinweis: Sehen Sie Sich dazu insbesondere das mittlere Bild der \(0\) an.
    Hinweis: Sie können ihre Überlegung prüfen, indem Sie die noch nicht gesehenen Daten (a.k.a. test set) nutzen und angeben wie viel Prozent der Beispiele ihr Klassifizierer richtig eingeordnet hat.