PyMC für nichtparametrisches Clustering: Der Dirichlet-Prozess zur Schätzung der Parameter der Gaußschen Mischung kann nicht geclustert werden

Problemeinrichtung

Eines der ersten Spielzeugprobleme, auf das ich PyMC anwenden wollte, ist das nichtparametrische Clustering: Modellieren Sie anhand einiger Daten diese als Gaußsche Mischung und lernen Sie die Anzahl der Cluster sowie den Mittelwert und die Kovarianz jedes Clusters. Das meiste, was ich über diese Methode weiß, stammt aus Videovorträgen von Michael Jordan und Yee Whye Teh, circa 2007 (bevor die Sparsamkeit zur Wut wurde) und den letzten Tagen, in denen Dr. Fonnesbecks und E. Chens Tutorials [fn1], [ fn2]. Das Problem ist jedoch gut untersucht und hat einige zuverlässige Implementierungen [fn3].

In diesem Spielzeugproblem generiere ich zehn Draws aus einem eindimensionalen Gaußschen und vierzig Draws aus N ( μ = 4 , σ = 2 ) . Wie Sie unten sehen können, habe ich die Ziehungen nicht gemischt, um leicht erkennen zu können, welche Proben von welcher Mischungskomponente stammen.$\mathcal{N}(\mu=0, \sigma=1)$$\mathcal{N}(\mu=4, \sigma=2)$

$y_i \sim \mathcal{N}(\mu_{z_i}, \sigma_{z_i})$$i=1,...,50$$z_i$$i$$z_i \in [1,...,N_{DP}]$$N_{DP}$$N_{DP}=50$

$z_i$$z_i \sim Categorical(p)$$p \sim Stick(\alpha)$$\alpha$$N_{DP}$$p$$N_{DP}$$\alpha$$\alpha$$\alpha \sim Uniform(0.3, 100)$

$N_{DP}$$\mu_{z_i}$$\sigma_{z_i}$$\mu_{z_i} \sim \mathcal{N}(\mu=0, \sigma=50)$$\sigma_{z_i} \sim Uniform(0, 100)$

$\mu_{z_i}$$\mu_{z_i} \sim \mathcal{N}(\mu_0, \sigma_0)$$\mu_0$$\sigma_0$

Zusammenfassend ist mein Modell:

$y_i \sim \mathcal{N}(\mu_{z_i}, \sigma_{z_i})$$i$

$z_i \sim Categorical(p)$$N_{DP}-1=49$$p \sim Stick(\alpha)$$N_{DP}$$\alpha \sim Uniform(0.3, 100)$, ein Skalar. (Ich bedaure jetzt etwas, dass die Anzahl der Datenproben der abgeschnittenen Länge des Dirichlets zuvor entspricht, aber ich hoffe, es ist klar.)

$\mu_{z_i} \sim \mathcal{N}(\mu=0, \sigma=50)$$\sigma_{z_i} \sim Uniform(0, 100)$$N_{DP}$$N_{DP}$

Hier ist das grafische Modell: Die Namen sind Variablennamen, siehe den folgenden Codeabschnitt.

Problemstellung

Trotz mehrerer Optimierungen und fehlgeschlagener Korrekturen ähneln die gelernten Parameter überhaupt nicht den wahren Werten, die die Daten generiert haben.

$z_i$$\alpha=5$

$p$$\mu_{z_i}$$i$$\mu_{z_i}$$i=1,...,20$

Unter Hinweis darauf, dass die ersten zehn Datenproben aus einem Modus stammten und der Rest aus dem anderen, kann das obige Ergebnis dies eindeutig nicht erfassen.

Wenn ich eine zufällige Initialisierung der Cluster-IDs zulasse, erhalte ich mehr als einen Cluster, aber der Cluster bedeutet, dass alle auf derselben 3.5-Ebene wandern:

$z_i$

$\alpha$

Anhang: Code

import pymc
import numpy as np

### Data generation

# Means and standard deviations of the Gaussian mixture model. The inference
# engine doesn't know these.
means = [0, 4.0]
stdevs = [1, 2.0]

# Rather than randomizing between the mixands, just specify how many
# to draw from each. This makes it really easy to know which draws
# came from which mixands (the first N1 from the first, the rest from
# the secon). The inference engine doesn't know about N1 and N2, only Ndata
N1 = 10
N2 = 40
Ndata = N1+N2

# Seed both the data generator RNG  as well as the global seed (for PyMC)
RNGseed = 123
np.random.seed(RNGseed)

def generate_data(draws_per_mixand):
    """Draw samples from a two-element Gaussian mixture reproducibly.

    Input sequence indicates the number of draws from each mixand. Resulting
    draws are concantenated together.

    """
    RNG = np.random.RandomState(RNGseed)
    values = np.hstack([RNG.normal(means[i], stdevs[i], ndraws)
                        for (i,ndraws) in enumerate(draws_per_mixand)])
    return values

observed_data = generate_data([N1, N2])


### PyMC model setup, step 1: the Dirichlet process and stick-breaking

# Truncation level of the Dirichlet process
Ndp = 50

# "alpha", or the concentration of the stick-breaking construction. There exists
# some interplay between choice of Ndp and concentration: a high concentration
# value implies many clusters, in turn implying low values for the leading
# elements of the probability mass function built by stick-breaking. Since we
# enforce the resulting PMF to sum to one, the probability of the last cluster
# might be then be set artificially high. This may interfere with the Dirichlet
# process' clustering ability.
#
# An example: if Ndp===4, and concentration high enough, stick-breaking might
# yield p===[.1, .1, .1, .7], which isn't desireable. You want to initialize
# concentration so that the last element of the PMF is less than or not much
# more than the a few of the previous ones. So you'd want to initialize at a
# smaller concentration to get something more like, say, p===[.35, .3, .25, .1].
#
# A thought: maybe we can avoid this interdependency by, rather than setting the
# final value of the PMF vector, scale the entire PMF vector to sum to 1? FIXME,
# TODO.
concinit = 5.0
conclo = 0.3
conchi = 100.0
concentration = pymc.Uniform('concentration', lower=conclo, upper=conchi,
                             value=concinit)

# The stick-breaking construction: requires Ndp beta draws dependent on the
# concentration, before the probability mass function is actually constructed.
betas = pymc.Beta('betas', alpha=1, beta=concentration, size=Ndp)

@pymc.deterministic
def pmf(betas=betas):
    "Construct a probability mass function for the truncated Dirichlet process"
    # prod = lambda x: np.exp(np.sum(np.log(x))) # Slow but more accurate(?)
    prod = np.prod
    value = map(lambda (i,u): u * prod(1.0 - betas[:i]), enumerate(betas))
    value[-1] = 1.0 - sum(value[:-1]) # force value to sum to 1
    return value

# The cluster assignments: each data point's estimated cluster ID.
# Remove idinit to allow clusterid to be randomly initialized:
idinit = np.zeros(Ndata, dtype=np.int64)
clusterid = pymc.Categorical('clusterid', p=pmf, size=Ndata, value=idinit)

### PyMC model setup, step 2: clusters' means and stdevs

# An individual data sample is drawn from a Gaussian, whose mean and stdev is
# what we're seeking.

# Hyperprior on clusters' means
mu0_mean = 0.0
mu0_std = 50.0
mu0_prec = 1.0/mu0_std**2
mu0_init = np.zeros(Ndp)
clustermean = pymc.Normal('clustermean', mu=mu0_mean, tau=mu0_prec,
                          size=Ndp, value=mu0_init)

# The cluster's stdev
clustersig_lo = 0.0
clustersig_hi = 100.0
clustersig_init = 50*np.ones(Ndp) # Again, don't really care?
clustersig = pymc.Uniform('clustersig', lower=clustersig_lo,
                          upper=clustersig_hi, size=Ndp, value=clustersig_init)
clusterprec = clustersig ** -2

### PyMC model setup, step 3: data

# So now we have means and stdevs for each of the Ndp clusters. We also have a
# probability mass function over all clusters, and a cluster ID indicating which
# cluster a particular data sample belongs to.

@pymc.deterministic
def data_cluster_mean(clusterid=clusterid, clustermean=clustermean):
    "Converts Ndata cluster IDs and Ndp cluster means to Ndata means."
    return clustermean[clusterid]

@pymc.deterministic
def data_cluster_prec(clusterid=clusterid, clusterprec=clusterprec):
    "Converts Ndata cluster IDs and Ndp cluster precs to Ndata precs."
    return clusterprec[clusterid]

data = pymc.Normal('data', mu=data_cluster_mean, tau=data_cluster_prec,
                   observed=True, value=observed_data)

Verweise

1. fn1: http://nbviewer.ipython.org/urls/raw.github.com/fonnesbeck/Bios366/master/notebooks/Section5_2-Dirichlet-Processes.ipynb
2. fn2: http://blog.echen.me/2012/03/20/infinite-mixture-models-with-nonparametric-bayes-and-the-dirichlet-process/
3. fn3: http://scikit-learn.org/stable/auto_examples/mixture/plot_gmm.html#example-mixture-plot-gmm-py

Ich bin mir nicht sicher, ob sich jemand mehr mit dieser Frage befasst, aber ich habe Ihre Frage an rjags gesendet, um Toms Gibbs-Stichprobenvorschlag zu testen, während ich Einblicke von Guy in die Wohnung vor der Standardabweichung einbeziehe.

Dieses Spielzeugproblem kann schwierig sein, da 10 und sogar 40 Datenpunkte nicht ausreichen, um die Varianz ohne einen informativen Vorgänger abzuschätzen. Die aktuelle vorherige σzi∼Uniform (0,100) ist nicht aussagekräftig. Dies könnte erklären, warum fast alle Zeichnungen von μzi das erwartete Mittel der beiden Verteilungen sind. Wenn es Ihre Frage nicht zu sehr ändert, werde ich 100 bzw. 400 Datenpunkte verwenden.

Ich habe den Stick Breaking-Prozess auch nicht direkt in meinem Code verwendet. Die Wikipedia- Seite für den Dirichlet-Prozess ließ mich denken, dass p ~ Dir (a / k) in Ordnung wäre.

Schließlich handelt es sich nur um eine semiparametrische Implementierung, da noch mehrere Cluster k erforderlich sind. Ich weiß nicht, wie man ein unendliches Mischungsmodell in rjags macht.

library("rjags")

set1 <- rnorm(100, 0, 1)
set2 <- rnorm(400, 4, 1)
data <- c(set1, set2)

plot(data, type='l', col='blue', lwd=3,
     main='gaussian mixture model data',
     xlab='data sample #', ylab='data value')
points(data, col='blue')

cpd.model.str <- 'model {
  a ~ dunif(0.3, 100)
  for (i in 1:k){
    alpha[i] <- a/k
    mu[i] ~ dnorm(0.0, 0.001)
    sigma[i] ~ dunif(0, 100)
  }
  p[1:k] ~ ddirich(alpha[1:k])
  for (i in 1:n){
    z[i] ~ dcat(p)
    y[i] ~ dnorm(mu[z[i]], pow(sigma[z[i]], -2))
  }
}' 


cpd.model <- jags.model(textConnection(cpd.model.str),
                        data=list(y=data,
                                  n=length(data),
                                  k=5))
update(cpd.model, 1000)
chain <- coda.samples(model = cpd.model, n.iter = 1000,
                      variable.names = c('p', 'mu', 'sigma'))
rchain <- as.matrix(chain)
apply(rchain, 2, mean)

Das schlechte Mischen, das Sie sehen, ist höchstwahrscheinlich auf die Art und Weise zurückzuführen, wie PyMC Samples zieht. Wie in Abschnitt 5.8.1 der PyMC-Dokumentation erläutert, werden alle Elemente einer Array-Variablen zusammen aktualisiert. In Ihrem Fall bedeutet dies, dass versucht wird, das gesamte clustermeanArray in einem Schritt zu aktualisieren clusterid. PyMC führt keine Gibbs-Abtastung durch. Es ist Metropolis, wo der Vorschlag durch einige einfache Heuristiken ausgewählt wird. Dies macht es unwahrscheinlich, einen guten Wert für ein gesamtes Array vorzuschlagen.