Problem

Chcę dopasować parametry modelu prostej populacji mieszanki 2-Gaussa. Biorąc pod uwagę cały szum wokół metod bayesowskich, chcę zrozumieć, czy w przypadku tego problemu wnioskowanie bayesowskie jest lepszym narzędziem niż tradycyjne metody dopasowywania.

Do tej pory MCMC radzi sobie bardzo słabo w tym przykładzie z zabawkami, ale może coś przeoczyłem. Zobaczmy więc kod.

Narzędzia

Użyję Pythona (2.7) + Scipy Stack, Lmfit 0.8 i PyMC 2.3.

Notatnik do odtworzenia analizy można znaleźć tutaj

Wygeneruj dane

Najpierw wygenerujmy dane:

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 sampleswygląda następująco:

„szeroki pik”, elementy trudno dostrzec wzrokiem.

Klasyczne podejście: dopasuj histogram

Spróbujmy najpierw klasycznego podejścia. Za pomocą lmfit łatwo jest zdefiniować model z dwoma pikami:

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'

Wreszcie dopasowujemy model do algorytmu simpleks:

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

Rezultat to następujący obraz (czerwone przerywane linie są dopasowanymi środkami):

Nawet jeśli problem jest dość trudny, przy odpowiednich wartościach początkowych i ograniczeniach modele zbiegły się do całkiem rozsądnego oszacowania.

Podejście bayesowskie: MCMC

Model w PyMC definiuję w sposób hierarchiczny. centersi sigmassą rozkładem pierwszeństwa dla hiperparametrów reprezentujących 2 centra i 2 sigmy 2 Gaussów. alphajest ułamkiem pierwszej populacji, a wcześniejsza dystrybucja to Beta.

Zmienna kategoryczna wybiera między dwiema populacjami. Rozumiem, że ta zmienna musi mieć taki sam rozmiar jak data ( samples).

Na koniec mui tausą zmienne deterministyczne które określają parametry rozkładu normalnego (zależą od categoryzmiennej więc losowo przełączanie pomiędzy dwiema wartościami dla dwóch populacji).

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])

Następnie uruchamiam MCMC z dość dużą liczbą iteracji (1e5, ~ 60s na moim komputerze):

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

Jednak wyniki są bardzo dziwne. Na przykład trace (część pierwszej populacji) ma tendencję do 0 zamiast zbieżności do 0,4 i ma bardzo silną autokorelację:$\alpha$

Także centra Gaussów również się nie zbiegają. Na przykład:

Jak widać w poprzednim wyborze, próbowałem „pomóc” algorytmowi MCMC przy użyciu rozkładu Beta dla poprzedniej frakcji populacji . Również wcześniejsze rozkłady dla centrów i sigm są dość rozsądne (tak myślę).$\alpha$

Co się tu dzieje? Czy robię coś źle, czy MCMC nie jest odpowiednie dla tego problemu?

Rozumiem, że metoda MCMC będzie wolniejsza, ale trywialne dopasowanie histogramu wydaje się działać znacznie lepiej w rozwiązywaniu populacji.

Problem jest spowodowany sposobem, w jaki PyMC pobiera próbki dla tego modelu. Jak wyjaśniono w sekcji 5.8.1 dokumentacji PyMC, wszystkie elementy zmiennej tablicowej są aktualizowane razem. W przypadku takich małych tablic centernie stanowi to problemu, ale w przypadku takiej dużej tablicy categoryprowadzi do niskiej częstotliwości akceptacji. Możesz zobaczyć wskaźnik akceptacji przez

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

Rozwiązaniem zaproponowanym w dokumentacji jest użycie tablicy zmiennych o wartości skalarnej. Ponadto należy skonfigurować niektóre rozkłady propozycji, ponieważ domyślne opcje są złe. Oto kod, który działa dla mnie:

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

proposal_sdI proposal_distributionopcje są opisane w rozdziale 5.7.1 . W przypadku centrów ustawiłem propozycję, aby z grubsza pasowała do standardowego odchylenia tylnej, która jest znacznie mniejsza niż domyślna ze względu na ilość danych. PyMC próbuje dostroić szerokość propozycji, ale działa to tylko wtedy, gdy wskaźnik akceptacji jest wystarczająco wysoki na początek. Na categoryprzykład wartość domyślna, proposal_distribution = 'Poisson'która daje słabe wyniki (nie wiem, dlaczego tak jest, ale z pewnością nie brzmi to jak rozsądna propozycja zmiennej binarnej).

Nie powinieneś modelować za pomocą Normalnej, w ten sposób dopuszczasz wartości ujemne dla standardowej odmiany. Zamiast tego użyj czegoś takiego jak:$\sigma$

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

aktualizacja:

Podszedłem do początkowych parametrów danych, zmieniając te części twojego modelu:

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

i wywołując mcmc z pewnym przerzedzeniem:

mcmc.sample(200000, 3000, 10)

wyniki:

Tylniaki nie są zbyt mili ... W rozdziale 11.6 książki o BŁĘDACH omawiają ten typ modelu i stwierdzają, że istnieją problemy z konwergencją bez oczywistego rozwiązania. Sprawdź także tutaj .

Również brak możliwości identyfikacji jest dużym problemem przy stosowaniu MCMC w modelach mieszanin. Zasadniczo, jeśli zmienisz etykiety środków klastra i przypisań klastra, prawdopodobieństwo się nie zmieni, a to może pomylić próbnik (między łańcuchami lub wewnątrz łańcuchów). Jedną z rzeczy, które możesz spróbować złagodzić, jest wykorzystanie potencjałów w PyMC3. Dobra implementacja GMM z potencjałem jest tutaj . Późniejsze dla tego rodzaju problemów jest również ogólnie wysoce multimodalne, co również stanowi duży problem. Zamiast tego możesz użyć EM (lub wnioskowania wariacyjnego).