Masalah

Saya ingin menyesuaikan parameter model populasi campuran 2-Gaussian sederhana. Mengingat semua hype di sekitar metode Bayesian saya ingin mengerti jika untuk masalah ini kesimpulan Bayesian adalah alat yang lebih baik daripada metode pemasangan tradisional.

Sejauh ini MCMC berkinerja sangat buruk dalam contoh mainan ini, tapi mungkin saya hanya mengabaikan sesuatu. Jadi mari kita lihat kodenya.

Alat-alatnya

Saya akan menggunakan python (2.7) + tumpukan scipy, lmfit 0.8 dan PyMC 2.3.

Buku catatan untuk mereproduksi analisis dapat ditemukan di sini

Hasilkan data

Pertama mari kita hasilkan data:

from scipy.stats import distributions

# Sample parameters
nsamples = 1000
mu1_true = 0.3
mu2_true = 0.55
sig1_true = 0.08
sig2_true = 0.12
a_true = 0.4

# Samples generation
np.random.seed(3)  # for repeatability
s1 = distributions.norm.rvs(mu1_true, sig1_true, size=round(a_true*nsamples))
s2 = distributions.norm.rvs(mu2_true, sig2_true, size=round((1-a_true)*nsamples))
samples = np.hstack([s1, s2])

Histogram samplesterlihat seperti ini:

"puncak luas", komponen-komponennya sulit dikenali oleh mata.

Pendekatan klasik: sesuai dengan histogram

Mari kita coba pendekatan klasiknya terlebih dahulu. Menggunakan lmfit mudah untuk mendefinisikan model 2-puncak:

import lmfit

peak1 = lmfit.models.GaussianModel(prefix='p1_')
peak2 = lmfit.models.GaussianModel(prefix='p2_')
model = peak1 + peak2

model.set_param_hint('p1_center', value=0.2, min=-1, max=2)
model.set_param_hint('p2_center', value=0.5, min=-1, max=2)
model.set_param_hint('p1_sigma', value=0.1, min=0.01, max=0.3)
model.set_param_hint('p2_sigma', value=0.1, min=0.01, max=0.3)
model.set_param_hint('p1_amplitude', value=1, min=0.0, max=1)
model.set_param_hint('p2_amplitude', expr='1 - p1_amplitude')
name = '2-gaussians'

Akhirnya kami menyesuaikan model dengan algoritma simpleks:

fit_res = model.fit(data, x=x_data, method='nelder')
print fit_res.fit_report()

Hasilnya adalah gambar berikut (garis putus-putus merah dilengkapi dengan pusat):

Sekalipun masalahnya agak sulit, dengan nilai awal dan kendala yang tepat, model-model tersebut bertemu dengan estimasi yang cukup masuk akal.

Pendekatan Bayesian: MCMC

Saya mendefinisikan model dalam PyMC secara hierarkis. centersdan sigmasadalah distribusi prior untuk hyperparameter yang mewakili 2 pusat dan 2 sigma dari 2 Gaussians. alphaadalah fraksi dari populasi pertama dan distribusi sebelumnya di sini adalah Beta.

Variabel kategorikal memilih antara dua populasi. Ini adalah pemahaman saya bahwa variabel ini harus berukuran sama dengan data ( samples).

Akhirnya mudan taumerupakan variabel deterministik yang menentukan parameter dari distribusi Normal (mereka bergantung pada categoryvariabel sehingga mereka secara acak beralih antara dua nilai untuk dua populasi).

sigmas = pm.Normal('sigmas', mu=0.1, tau=1000, size=2)
centers = pm.Normal('centers', [0.3, 0.7], [1/(0.1)**2, 1/(0.1)**2], size=2)
#centers = pm.Uniform('centers', 0, 1, size=2)

alpha  = pm.Beta('alpha', alpha=2, beta=3)
category = pm.Categorical("category", [alpha, 1 - alpha], size=nsamples)

@pm.deterministic
def mu(category=category, centers=centers):
    return centers[category]

@pm.deterministic
def tau(category=category, sigmas=sigmas):
    return 1/(sigmas[category]**2)

observations = pm.Normal('samples_model', mu=mu, tau=tau, value=samples, observed=True)
model = pm.Model([observations, mu, tau, category, alpha, sigmas, centers])

Kemudian saya menjalankan MCMC dengan jumlah iterasi yang cukup lama (1e5, ~ 60s pada mesin saya):

mcmc = pm.MCMC(model)
mcmc.sample(100000, 30000)

Namun hasilnya sangat aneh. Misalnya trace (fraksi populasi pertama) cenderung 0 untuk konvergen ke 0,4 dan memiliki autokorelasi yang sangat kuat:$\alpha$

Juga pusat-pusat Gaussians tidak bertemu juga. Sebagai contoh:

Seperti yang Anda lihat dalam pilihan sebelumnya, saya mencoba "membantu" algoritma MCMC menggunakan distribusi Beta untuk fraksi populasi sebelumnya . Juga distribusi sebelumnya untuk pusat dan sigma cukup masuk akal (saya pikir).$\alpha$

Jadi apa yang terjadi di sini? Apakah saya melakukan sesuatu yang salah atau MCMC tidak cocok untuk masalah ini?

Saya mengerti bahwa metode MCMC akan lebih lambat, tetapi kecocokan histogram sepele tampaknya berkinerja jauh lebih baik dalam menyelesaikan populasi.

Masalahnya disebabkan oleh cara PyMC mengambil sampel untuk model ini. Seperti dijelaskan di bagian 5.8.1 dari dokumentasi PyMC, semua elemen variabel array diperbarui bersama-sama. Untuk array kecil seperti centerini bukan masalah, tetapi untuk array besar seperti categoryitu mengarah pada tingkat penerimaan yang rendah. Anda dapat melihat tingkat penerimaan melalui

print mcmc.step_method_dict[category][0].ratio

Solusi yang disarankan dalam dokumentasi adalah dengan menggunakan array variabel skalar-dihargai. Selain itu, Anda perlu mengonfigurasi beberapa distribusi proposal karena pilihan defaultnya buruk. Berikut adalah kode yang berfungsi untuk saya:

import pymc as pm
sigmas = pm.Normal('sigmas', mu=0.1, tau=1000, size=2)
centers = pm.Normal('centers', [0.3, 0.7], [1/(0.1)**2, 1/(0.1)**2], size=2)
alpha  = pm.Beta('alpha', alpha=2, beta=3)
category = pm.Container([pm.Categorical("category%i" % i, [alpha, 1 - alpha]) 
                         for i in range(nsamples)])
observations = pm.Container([pm.Normal('samples_model%i' % i, 
                   mu=centers[category[i]], tau=1/(sigmas[category[i]]**2), 
                   value=samples[i], observed=True) for i in range(nsamples)])
model = pm.Model([observations, category, alpha, sigmas, centers])
mcmc = pm.MCMC(model)
# initialize in a good place to reduce the number of steps required
centers.value = [mu1_true, mu2_true]
# set a custom proposal for centers, since the default is bad
mcmc.use_step_method(pm.Metropolis, centers, proposal_sd=sig1_true/np.sqrt(nsamples))
# set a custom proposal for category, since the default is bad
for i in range(nsamples):
    mcmc.use_step_method(pm.DiscreteMetropolis, category[i], proposal_distribution='Prior')
mcmc.sample(100)  # beware sampling takes much longer now
# check the acceptance rates
print mcmc.step_method_dict[category[0]][0].ratio
print mcmc.step_method_dict[centers][0].ratio
print mcmc.step_method_dict[alpha][0].ratio

The proposal_sddan proposal_distributionpilihan dijelaskan dalam bagian 5.7.1 . Untuk center, saya mengatur proposal agar sesuai dengan standar deviasi posterior, yang jauh lebih kecil daripada standar karena jumlah data. PyMC memang berusaha untuk menyempurnakan lebar proposal, tetapi ini hanya berfungsi jika tingkat penerimaan Anda cukup tinggi untuk memulai. Karena category, default proposal_distribution = 'Poisson'yang memberikan hasil yang buruk (saya tidak tahu mengapa ini, tapi tentu saja tidak terdengar seperti proposal yang masuk akal untuk variabel biner).

Anda tidak boleh memodelkan dengan Normal, dengan begitu Anda mengizinkan nilai negatif untuk variasi standar. Gunakan sesuatu seperti:$\sigma$

sigmas = pm.Exponential('sigmas', 0.1, size=2)

memperbarui:

Saya mendekati parameter awal data dengan mengubah bagian-bagian dari model Anda:

sigmas = pm.Exponential('sigmas', 0.1, size=2)
alpha  = pm.Beta('alpha', alpha=1, beta=1)

dan dengan memanggil mcmc dengan beberapa penipisan:

mcmc.sample(200000, 3000, 10)

hasil:

Posisinya tidak terlalu bagus ... Di bagian 11.6 dari BUGS Book mereka membahas jenis model ini dan menyatakan bahwa ada masalah konvergensi tanpa solusi yang jelas. Periksa juga di sini .

Selain itu, non-identifikasi adalah masalah besar untuk menggunakan MCMC untuk model campuran. Pada dasarnya, jika Anda beralih label pada sarana dan tugas gugus kluster Anda, kemungkinan tidak berubah, dan ini dapat membingungkan sampler (antara rantai atau dalam rantai). Satu hal yang dapat Anda coba untuk mengurangi ini adalah menggunakan Potensi di PyMC3. Implementasi GMM dengan Potensi yang baik ada di sini . Posterior untuk masalah-masalah semacam ini juga umumnya sangat multimodal, yang juga menghadirkan masalah besar. Anda mungkin ingin menggunakan EM (atau Variational inference).