PyMC untuk pengelompokan nonparametrik: Proses Dirichlet untuk memperkirakan parameter campuran Gaussian gagal dikelompokkan

10

Pengaturan masalah

Salah satu masalah mainan pertama yang ingin saya terapkan pada PyMC adalah pengelompokan nonparametrik: diberi beberapa data, modelkan sebagai campuran Gaussian, dan pelajari jumlah cluster dan rata-rata dan kovarian setiap cluster. Sebagian besar yang saya ketahui tentang metode ini berasal dari ceramah video oleh Michael Jordan dan Yee Whye Teh, sekitar 2007 (sebelum sparsity menjadi kemarahan), dan beberapa hari terakhir membaca tutorial Dr Fonnesbeck dan E. Chen [fn1], [ fn2]. Tetapi masalahnya dipelajari dengan baik dan memiliki beberapa implementasi yang dapat diandalkan [fn3].

Dalam masalah mainan ini, saya menghasilkan sepuluh undian dari Gaussian satu dimensi ( μ = 0 , σ = 1 ) dan empat puluh undian dari N ( μ = 4 , σ = 2 ) . Seperti yang dapat Anda lihat di bawah, saya tidak mengacak gambarnya, untuk memudahkan menentukan sampel mana yang berasal dari komponen campuran mana.N(μ=0,σ=1)N(μ=4,σ=2)

Data model campuran Gaussian

Saya memodelkan setiap sampel data yang , untuk i = 1 , . . . , 50 dan di mana z i menunjukkan cluster untuk ini saya th data titik: z i[ 1 , . . . , N D P ] . N D P di sini adalah panjang dari proses Dirichlet terpotong yang digunakan: bagi saya, NyiN(μzi,σzi)i=1,...,50ziizi[1,...,NDP]NDP.NDP=50

Memperluas infrastruktur proses Dirichlet, masing-masing klaster ID adalah menarik dari variabel acak kategoris, yang probabilitas massa fungsi diberikan oleh konstruk tongkat-melanggar: z i ~ C a t e g o r i c a l ( p ) dengan p S t i c k ( α ) untuk parameter konsentrasi α . Stick-breaking membangun N D P -long vector pziziCategorical(p)pStick(α)αNDPp, yang harus berjumlah 1, dengan terlebih dahulu mendapatkan iid Undian terdistribusi-beta yang bergantung pada α , lihat [fn1]. Dan karena saya ingin data untuk menginformasikan ketidaktahuan saya tentang α , saya mengikuti [fn1] dan menganggap α U n i f o r m ( 0,3 , 100 ) .NDPαααUniform(0.3,100)

Ini menentukan bagaimana ID cluster masing-masing sampel data dihasilkan. Masing-masing dari cluster memiliki mean dan deviasi standar terkait, μ z i dan σ z i . Kemudian, μ z iN ( μ = 0 , σ = 50 ) dan σ z iU n i f o r m ( 0 , 100 ) .NDPμziσziμziN(μ=0,σ=50)σziUniform(0,100)

μziμziN(μ0,σ0)μ0σ0

Singkatnya, model saya adalah:

yiN(μzi,σzi)i

ziCategorical(p)NDP1=49pStick(α)NDPαUniform(0.3,100), skalar. (Sekarang saya sedikit menyesal membuat jumlah sampel data sama dengan panjang Dirichlet yang terpotong sebelumnya, tapi saya harap ini jelas.)

μziN(μ=0,σ=50)σziUniform(0,100)NDPNDP

Berikut model grafis: nama-nama adalah nama variabel, lihat bagian kode di bawah ini.

Grafik

Pernyataan masalah

Meskipun ada beberapa penyesuaian dan perbaikan yang gagal, parameter yang dipelajari sama sekali tidak sama dengan nilai sebenarnya yang menghasilkan data.

ziα=5

pμziiμzii=1,...,20

Berarti dengan ID cluster nol yang diinisialisasi

Mengingat bahwa sepuluh sampel data pertama berasal dari satu mode dan sisanya dari yang lain, hasil di atas jelas gagal menangkapnya.

Jika saya mengizinkan inisialisasi acak dari ID cluster, maka saya mendapatkan lebih dari satu cluster tetapi cluster berarti semua berjalan di sekitar level 3,5 yang sama:

Berarti dengan ID cluster inisialisasi acak

zi

α

Lampiran: kode

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)

Referensi

  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
sumber
Varians komponen sebelumnya Anda adalah Seragam (0,100) yang dapat menyebabkan Anda mengalami masalah besar. Hanya 2% dari massa sebelumnya ini mencakup varian sebenarnya dari 1 dan 2. Varian yang diharapkan dari komponen Anda di bawah ini adalah 50, yang merupakan Gaussian luas sehingga dapat dengan mudah memperhitungkan data Anda dengan satu komponen.
jerad
Sudahkah Anda membaca bab Pemrograman Probabilistik dan Statistik Bayesian untuk Hacker ini? Ini memiliki contoh yang dapat membantu Anda!
Tim
Ini sepertinya sedikit singkat untuk dijawab. Tampaknya lebih dari komentar. Bisakah Anda menguraikan informasi OP apa yang akan didapat dengan membacanya?
Glen_b -Reinstate Monica
@TimRich ya saya sudah membacanya, dan mengambil kelas di sekolah pascasarjana, dan bekerja di industri melakukan statistik terapan;) ini adalah pertanyaan khusus PyMC.
Ahmed Fasih
1
μZiN(μ0,σ0)μ0,σ0

Jawaban:

4

Saya tidak yakin apakah ada yang melihat pertanyaan ini lagi, tetapi saya mengajukan pertanyaan Anda ke rjags untuk menguji saran sampling Gibbs Tom sambil memasukkan wawasan dari Guy tentang flat sebelum standar deviasi.

Masalah mainan ini mungkin sulit karena 10 dan bahkan 40 titik data tidak cukup untuk memperkirakan variasi tanpa informasi sebelumnya. Σzi∼Uniform sebelumnya (0,100) sebelumnya tidak informatif. Ini mungkin menjelaskan mengapa hampir semua gambar μzi adalah rata-rata yang diharapkan dari dua distribusi. Jika tidak terlalu banyak mengubah pertanyaan Anda, saya akan menggunakan masing-masing 100 dan 400 poin data.

Saya juga tidak menggunakan proses pemecahan tongkat secara langsung dalam kode saya. The wikipedia halaman untuk proses Dirichlet membuat saya berpikir p ~ Dir (a / k) akan ok.

Akhirnya itu hanya implementasi semi-parametrik karena masih membutuhkan sejumlah cluster k. Saya tidak tahu cara membuat model campuran tak terbatas dalam rjags.

markov chain mu cluster 1

markov chain mu cluster 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)
pengguna1068430
sumber
1
KKαi=a/K500K=25JAGSJAGSp
1

Pencampuran yang buruk yang Anda lihat kemungkinan besar karena cara PyMC menarik sampel. Seperti dijelaskan di bagian 5.8.1 dari dokumentasi PyMC, semua elemen variabel array diperbarui bersama-sama. Dalam kasus Anda, itu berarti akan mencoba memperbarui seluruh clustermeanarray dalam satu langkah, dan demikian pula untuk clusterid. PyMC tidak melakukan sampling Gibbs; itu tidak Metropolis di mana proposal dipilih oleh beberapa heuristik sederhana. Ini membuatnya tidak mungkin untuk mengusulkan nilai yang baik untuk seluruh array.

Tom Minka
sumber
Segera setelah Anda berkata, "ini akan mencoba memperbarui seluruh array dalam satu langkah", saya memahami kelemahan Metropolis (dalam hal ini) versus Gibbs. Apakah ada sesuatu yang istimewa tentang STAN atau JAGS yang mungkin membuat mereka melakukan yang lebih baik dalam hal ini? Dalam kedua kasus, saya akan menghabiskan waktu untuk mengimplementasikan Gibbs di PyMC. Terima kasih! (Saya sudah menjadi penggemar pekerjaan Anda sejak lightspeed, jadi terima kasih dua kali lipat!)
Ahmed Fasih
1
STAN tidak menangani variabel diskrit, tetapi JAGS patut dicoba.
Tom Minka