Skip to Tutorial Content

Bayesianische Mischverteilungsmodelle

Mischverteilungen

  • (Normal-)Verteilungsannahme basiert in der Regel auf der Annahme, dass die Beobachtungen identisch verteilt sind, es also eine Verteilung der Grundgesamtheit gibt
  • Alternative Annahme: Es gibt mehrere Gruppen in der Grundgesamheit, die alle die selbe Verteilungsklasse haben, aber unterschiedliche Parameter
  • Einfachster Fall: Mischung von Normalverteilungen

Mischung von Normalverteilungen

  • Natürlich sind auch verschiedene Varianzen und verschiedene Gewichte möglich

Beispiel

Länge von 256 Fischen (aus D. M. Titterington, A. F. M. Smith and U.E. Makov (1985) Statistical Analysis of Finite Mixture Distributions. Wiley.)

Fragen

  • Welche Parameter und welche Gewichte haben die einzelnen Verteilungen?
  • Wieviele Mischverteilungen brauchen wir?
  • Zu welcher Verteilung gehört welche Beobachtung

Latente Klassen

Seien \(i=1,\ldots,n\) die Beobachtungen.

Idee der Modellierung: Führe latente Klassen \(S_i\) ein. Latent heißt dabei, das wir die Klassenzugehörigkeit nicht beobachten können, wir sie aber im Modell brauchen.

\[ \begin{aligned} x_i | S_i & \sim N(\mu_{S_i},\sigma^2_{S_i})\\ \mu_{j} &\sim N(\mu_0, \sigma^2_0)\\ \sigma^2_j &\sim IG(a,b)\\ \end{aligned} \]

Jede Beobachtung (jeder Fisch) \(i\) hat also eine Klasse \(S_i\). Sprich \(S_i\) zeigt uns an, zu welcher Verteilung die Beobachtung \(i\) gehört. Jede Klasse hat einen eigenen Erwartungswert und eine eigene Varianz. Gehört der Fisch der Klasse \(j\) an, erwarten wir im Mittel also ein Gewicht von \(\mu_j\).

Die Klassenzugehörigkeiten \(S_i\) sind uns nicht bekannt. Wir brauchen also eine Verteilung auf \(S_i\), um \(S_i\) Bayesianische zu schätzen.

Da \(S_i\) diskret, ist einfach \(P(S_i=j) = \pi_j\)

Die Wahrscheinlichkeiten der Klassenzugehörigkeiten kennen wir aber auch nicht. Um diese zu schätzen, brauchen wir wieder eine Priori-Verteilung. Hier bietet sich die Dirichlet-Verteilung an:

\[ (\pi_1,\ldots,\pi_K) \sim Diri(\alpha,\ldots,\alpha) \]

Die Dichte der Dirichlet-Verteilung ist

\[ p(\pi_1,\ldots,\pi_K) \propto \prod \pi_i^{\alpha_i-1} \]

Für die Dirichletverteilung gilt außerdem \(0<\pi_k<1\) für alle \(k\) und \(\sum_{k=1}^K\pi_k=1\), was genau den Eigenschaften von Wahrscheinlichkeiten entspricht.

Full conditionals

  • Die full conditional von \(\mu_j\) sind Normalverteilungen, wobei die \(x_i\) mit \(S_i=j\) eingehen
  • Analog ist die full conditional von \(\sigma^2_j\) eine IG-Verteilung
  • \(p(S_i=j)\) berechnet sich aus Dichte von \(x_i\) mit \(\mu_j\) und \(\sigma^2_j\) (und Priori)
  • Full Conditional von \(\mathbf{\pi}\) ist Dirichlet mit \((\alpha+n_1,\ldots,\alpha+n_k)\), wobei \(n_j\) die aktuelle Anzahl der Beobachtungen in Klasse \(j\) ist

Das heißt, hier kann ein Gibbs-Sampler benutzt werden.

bayesmix-Paket

Das bayesmix-Paket erstellt ein entsprechendes Mischverteilungsmodell und benutzt zur Durchführung des MCMC-Algorithmus JAGS.

library(bayesmix)
model <- BMMmodel(fish, k = 4,
  initialValues = list(S0 = 2),
  priors = list(kind = "independence",
  parameter = "priorsFish", 
  hierarchical = "tau"))
## Data for nodes: b0, B0inv, nu0Half, g0Half, g0G0Half, k, N, e, y
## Initial values for nodes: eta, mu, tau, S0
## 
## Model specification in BUGS language:
## 
## var 
##  b0,
##      B0inv,
##      nu0Half,
##      g0Half,
##      g0G0Half,
##      k,
##      N,
##      eta[4], 
##  mu[4], 
##  tau[4], 
##  nu0S0Half,
##      S0,
##      e[4], 
##  y[256], 
##  S[256]; 
## 
## model    {
##  for (i in 1:N) {
##      y[i] ~ dnorm(mu[S[i]],tau[S[i]]);
##      S[i] ~ dcat(eta[]);
##  }
##  for (j in 1:k) {
##      mu[j] ~ dnorm(b0,B0inv);
##      tau[j] ~ dgamma(nu0Half,nu0S0Half);
##  }
##  S0 ~ dgamma(g0Half,g0G0Half);
##  nu0S0Half <- nu0Half * S0;
## 
##  eta[] ~ ddirch(e[]);
## }

## Compiling model graph
##    Declaring variables
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 256
##    Unobserved stochastic nodes: 266
##    Total graph size: 1047
## 
## Initializing model

Wir sortieren die vier Klassen noch nach ihrem Erwartungswert:

## 
## Call:
## JAGSrun(y = fish, model = model, control = control)
## 
## Markov Chain Monte Carlo (MCMC) output:
## Start = 1001 
## End = 6000 
## Thinning interval = 1 
## 
##  Empirical mean, standard deviation and 95% CI for eta 
##          Mean      SD    2.5%  97.5%
## eta[1] 0.1289 0.07699 0.07061 0.3425
## eta[2] 0.5003 0.07249 0.34015 0.5962
## eta[3] 0.2702 0.05194 0.15888 0.3669
## eta[4] 0.1007 0.03242 0.05115 0.1751
## 
##  Empirical mean, standard deviation and 95% CI for mu 
##        Mean     SD  2.5%  97.5%
## mu[1] 3.444 0.3656 3.128  5.005
## mu[2] 5.299 0.1868 5.129  5.481
## mu[3] 7.436 0.2296 7.090  7.720
## mu[4] 9.954 0.4480 8.819 10.658
## 
##  Empirical mean, standard deviation and 95% CI for sigma2 
##             Mean     SD   2.5% 97.5%
## sigma2[1] 0.3477 0.3012 0.1485 1.405
## sigma2[2] 0.4055 0.3024 0.2403 1.377
## sigma2[3] 0.5327 0.4037 0.2389 1.400
## sigma2[4] 1.0431 0.4775 0.5389 2.298
  • eta entspricht hier \(\pi_k\). Hier tritt die Klasse 2 am häufigsten auf.
  • mu ist der Erwartungswert pro Klasse
  • sigma2 ist die Varianz. Diese steigt hier mit dem Erwartungswert.

Medianmodell

  • Wir berechnen Punktschätzer aus den Medianen der Ziehungen.
  • Daraus können wir die (Median)-Misch-Verteilung wieder zusammensetzen

Die ersten 256 Parameter sind dabei die \(S_i\), also die Klassenzugehörigkeiten. Danach kommen eta, mu und sigma2.

Posteriori-Verteilung der Klassen

  • Schließlich können wir uns noch für jede Beobachtung die Posteriori-Wahrscheinlichkeit für die Zugehörigkeit zu den Klassen 1 bis 4 anzeigen lassen.
  • Rot entspricht Klasse 1, Grün Klasse 2, Blau Klasse 3 und Orange Klasse 4.

Weiter

Bayesianische Mischverteilungsmodelle