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 samples
wyglą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. centers
i sigmas
są rozkładem pierwszeństwa dla hiperparametrów reprezentujących 2 centra i 2 sigmy 2 Gaussów. alpha
jest 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 mu
i tau
są zmienne deterministyczne które określają parametry rozkładu normalnego (zależą od category
zmiennej 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ę:
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ę).
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.
źródło
proposal_distribution
iproposal_sd
dlaczego przy użyciuPrior
jest lepsze dla zmiennych kategorycznych?Nie powinieneś modelować za pomocą Normalnej, w ten sposób dopuszczasz wartości ujemne dla standardowej odmiany. Zamiast tego użyj czegoś takiego jak:σ
aktualizacja:
Podszedłem do początkowych parametrów danych, zmieniając te części twojego modelu:
i wywołując mcmc z pewnym przerzedzeniem:
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ódło
mcmc.sample(60000, 3000, 3)
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).
źródło