PyMC dla grupowania nieparametrycznego: proces Dirichleta do oszacowania parametrów mieszanki Gaussa nie ulega zgrupowaniu

10

Konfiguracja problemu

Jednym z pierwszych problemów z zabawkami, do których chciałem zastosować PyMC, jest grupowanie nieparametryczne: biorąc pod uwagę pewne dane, zamodeluj je jako mieszaninę Gaussa i poznaj liczbę skupień oraz średnią i kowariancję każdego skupienia. Większość tego, co wiem o tej metodzie, pochodzi z wykładów wideo Michaela Jordana i Yee Whye Teh, około 2007 r. (Zanim sparsity stało się wściekłością), oraz ostatnich kilku dni, w których przeczytałem tutoriale dra Fonnesbecka i E. Chena [fn1], [ fn2]. Ale problem jest dobrze zbadany i ma pewne niezawodne implementacje [fn3].

W tym problemie z zabawkami generuję dziesięć losowań z jednowymiarowego Gaussa i czterdzieści losowań z . Jak widać poniżej, nie losowałem losowań, aby ułatwić stwierdzenie, które próbki pochodzą z którego składnika mieszaniny.N ( μ = 4 , σ = 2 )N(μ=0,σ=1)N(μ=4,σ=2)

Dane modelu mieszanki Gaussa

każdą próbkę danych dla i gdzie wskazuje klaster dla tego tego punktu danych: . tutaj jest długość zastosowanego skróconego procesu Dirichleta: dla mnie .i = 1 , . . . , 50 oo i i oo i[ 1 , . . . , N D P ] N D P N D P = 50yiN(μzi,σzi)i=1,...,50ziizi[1,...,NDP]NDPNDP=50

Rozbudowując infrastrukturę procesu Dirichleta, każdy klastra jest czerpaniem z kategorycznej zmiennej losowej, której funkcja masy prawdopodobieństwa jest dana przez konstrukcję : z dla parametr stężenia . Łamanie konstruktów konstruuje wektor , który musi sumować się do 1, najpierw uzyskując iid losowania rozproszone w wersji beta, które zależą od , patrz [fn1]. A ponieważ chciałbym, aby dane informowały moją ignorancję o , podążam za [fn1] i zakładam .z iC a t e g o r i c a l ( p ) p S t i c k ( α ) α N D P p N D P α α α U n i f o r m ( 0,3 , 100 )ziziCategorical(p)pStick(α)αNDPpNDPαααUniform(0.3,100)

Określa sposób generowania identyfikatora klastra każdej próbki danych. Każdy klaster ma powiązaną średnią i odchylenie standardowe, i . Następnie i . μ z i σ z i μ z iN ( μ = 0 , σ = 50 ) σ z iU n i f o r m ( 0 , 100 )NDPμziσziμziN(μ=0,σ=50)σziUniform(0,100)

(Wcześniej bezmyślnie śledziłem [fn1] i umieszczałem hiperprior na , czyli z samym losowaniem z rozkład normalny o ustalonych parametrach i z munduru. Ale według https://stats.stackexchange.com/a/71932/31187 moje dane nie obsługują tego rodzaju hierarchicznego hiperpriora.) μ z iN ( μ 0 , σ 0 ) μ 0 σ 0μziμziN(μ0,σ0)μ0σ0

Podsumowując, mój model to:

yiN(μzi,σzi) gdzie działa od 1 do 50 (liczba próbek danych).i

ziCategorical(p) i może przyjmować wartości od 0 do ; , długi wektor ; i , skalar. (Teraz trochę żałuję, że wcześniej liczba próbek danych była równa skróconej długości Dirichleta, ale mam nadzieję, że to jasne).NDP1=49pStick(α)NDPαUniform(0.3,100)

μziN(μ=0,σ=50) i . Istnieje tych średnich i odchyleń standardowych (po jednym dla każdego z możliwych klastrów.)σziUniform(0,100)NDPNDP

Oto model graficzny: nazwy są nazwami zmiennych, patrz sekcja kodu poniżej.

Wykres

Opis problemu

Pomimo kilku poprawek i nieudanych poprawek wyuczone parametry wcale nie są podobne do prawdziwych wartości, które wygenerowały dane.

Obecnie inicjuję większość zmiennych losowych na stałe wartości. Zmienne średnie i odchylenie standardowe są inicjowane do wartości oczekiwanych (tj. 0 dla wartości normalnych, środek ich wsparcia dla wartości jednolitych). wszystkie identyfikatory klastrów na 0. I zainicjowałem parametr stężenia .ziα=5

Przy takich inicjalizacjach iteracje 100 000 MCMC po prostu nie mogą znaleźć drugiego klastra. Pierwszy element jest bliski 1 i prawie wszystkie losowania dla wszystkich próbek danych są takie same, około 3,5. tutaj co 100 losowanie dla pierwszych dwudziestu próbek danych, tj. dla :pμziiμzii=1,...,20

Oznacza z zainicjowanym zerem ID klastra

Przypominając, że pierwsze dziesięć próbek danych pochodziło z jednego trybu, a reszta z drugiego, powyższy wynik wyraźnie tego nie wychwycił.

Jeśli pozwolę na losową inicjalizację identyfikatorów klastra, uzyskam więcej niż jeden klaster, ale klaster oznacza, że ​​wszystkie wędrują wokół tego samego poziomu 3.5:

Oznacza z losowo inicjowanym identyfikatorem klastra

To sugeruje mi, że jest to zwykły problem z MCMC, że nie może osiągnąć innego trybu a posteriori od tego, w którym się znajduje: pamiętaj, że te różne wyniki występują po zmianie inicjalizacji identyfikatorów klastra , a nie ich priorytetów lub coś jeszcze.zi

Czy popełniam jakieś błędy modelowania? Podobne pytanie: https://stackoverflow.com/q/19114790/500207 chce użyć rozkładu Dirichleta i dopasować 3-elementową mieszaninę Gaussa i napotyka nieco podobne problemy. Czy powinienem rozważyć utworzenie modelu w pełni sprzężonego i stosowanie próbkowania Gibbs dla tego rodzaju grupowania? (Wdrożyłem próbnik Gibbsa dla parametrycznego przypadku dystrybucji Dirichleta, z wyjątkiem użycia stałej koncentracji w przeszłości i działało to dobrze, więc spodziewaj się, że PyMC będzie w stanie rozwiązać przynajmniej ten problem ręcznie.)α

Dodatek: kod

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)

Bibliografia

  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
Ahmed Fasih
źródło
Twój priorytet w stosunku do wariancji komponentów to Uniform (0,100), co może powodować poważne problemy. Tylko 2% masy tego uprzedniego obejmuje prawdziwe wariancje 1 i 2. Oczekiwana wariancja twoich składników w ramach tego wcześniejszego wynosi 50, co jest tak szerokim gaussowskim, że można łatwo uwzględnić twoje dane za pomocą jednego składnika.
jerad
Czy przeczytałeś już ten rozdział książki Probabilistic Programming and Bayesian Statistics for Hackers ? Ma przykład, który może ci pomóc!
Tim
Odpowiedź wydaje się trochę krótka. To wydaje się bardziej komentarzem. Czy możesz przynajmniej nakreślić, jakie informacje OP otrzyma po ich przeczytaniu?
Glen_b
@TimRich tak, przeczytałem go i wziąłem klasę na studiach podyplomowych i pracowałem w branży, wykonując statystyki stosowane;) To pytanie dotyczy PyMC.
Ahmed Fasih,
1
Nie odrzuciłbym hierarchicznego przeora. Powszechnie wiadomo, że stawiasz się w złej pozycji, jeśli stawiasz płaskie priory na składnikach mieszanki - szczególnie, gdy jednocześnie próbujesz poznać liczbę klastrów. To właśnie powoduje dziwne skoki na twoich wykresach śladowych. Wydaje się, że wielkie nazwy w NP-Bayes ustawiają i albo używają szacunkowych wtyczek albo umieszczając informacje na temat tych komponentów zaprojektowanych do umieść je w odpowiednich skalach. μ 0 , σ 0μZiN(μ0,σ0)μ0,σ0
facet

Odpowiedzi:

4

Nie jestem pewien, czy ktoś już patrzy na to pytanie, ale włożyłem to pytanie do rjags, aby przetestować sugestię Toma dotyczącą próbkowania Gibbsa, jednocześnie uwzględniając wgląd od Guy'a o mieszkanie przed standardowym odchyleniem.

Ten problem z zabawkami może być trudny, ponieważ 10, a nawet 40 punktów danych nie jest wystarczających do oszacowania wariancji bez informacyjnego uprzedniego. Obecny poprzedni σzi∼Uniform (0,100) nie ma charakteru informacyjnego. To może wyjaśniać, dlaczego prawie wszystkie losowania μzi są oczekiwaną średnią z dwóch rozkładów. Jeśli nie zmieni to zbytnio twojego pytania, użyję odpowiednio 100 i 400 punktów danych.

Nie użyłem również procesu łamania pamięci bezpośrednio w kodzie. Strona Wikipedii dla procesu dirichleta sprawiła, że ​​pomyślałem, że p ~ Dir (a / k) będzie w porządku.

Wreszcie jest to tylko implementacja półparametryczna, ponieważ nadal wymaga wielu klastrów k. Nie wiem, jak zrobić model nieskończonej mieszanki w rjagach.

klaster mu łańcucha markowa 1

klaster mu łańcucha markowa 2

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)
użytkownik1068430
źródło
1
Jeśli chodzi o użycie skończonej liczby klastrów , możesz po prostu wziąć duże (i ustawić , jak masz) i uzyskać model, który jest faktycznie identyczny z procesem Dirichleta. Zbieżność tego przed procesem Dirichleta przedtem jest bardzo szybka; na obserwacji coś w rodzaju powinno wystarczyć. Jest to także mądry, podczas używania , aby nie używać budowę stick-breaking, bo nie wykryje aktualizację bloku na ; to, co tutaj zrobiłeś, jest lepsze. K α i = a / K 500 K = 25 pKKαi=a/K500K=25JAGSJAGSp
facet
1

Słabe miksowanie, które widzisz, jest najprawdopodobniej spowodowane sposobem, w jaki PyMC pobiera próbki. Jak wyjaśniono w sekcji 5.8.1 dokumentacji PyMC, wszystkie elementy zmiennej tablicowej są aktualizowane razem. W twoim przypadku oznacza to, że spróbuje zaktualizować całą clustermeantablicę w jednym kroku i podobnie clusterid. PyMC nie wykonuje próbkowania Gibbsa; robi to w Metropolis, gdzie propozycję wybiera prosta heurystyka. To sprawia, że ​​jest mało prawdopodobne, aby zaproponować dobrą wartość dla całej tablicy.

Tom Minka
źródło
Jak tylko powiedziałeś: „spróbuje zaktualizować całą tablicę w jednym kroku”, zrozumiałem wady Metropolis (w tym przypadku) w porównaniu z Gibbsem. Czy jest coś specjalnego w STAN lub JAGS, co może im pomóc w tym lepiej? W obu przypadkach zamierzam poświęcić trochę czasu na wdrożenie Gibbs w PyMC. Dziękuję Ci! (Byłem fanem twojej pracy od prędkości światła, więc podwójnie dziękuję!)
Ahmed Fasih
1
STAN nie obsługuje zmiennych dyskretnych, ale JAGS jest warty wypróbowania.
Tom Minka