Первое использование PyMc не работает

Question

Я новичок в PyMc и хотел бы знать, почему этот код не работает. Я уже потратил на это несколько часов, но я что-то упустил. Может ли кто-нибудь мне помочь?

вопрос, который я хочу обратиться:

• У меня есть набор мер НПЦ, которые показывают 3 ударов, поэтому я хочу, чтобы смоделировать это как сумма 3 гауссиан (предполагается, что меры, шумные и гауссовы approx ok) ==> Я хочу оценить 8 параметров: относительные веса ударов (т. е. 2 параметра), их 3 средства и их 3 отклонения.

• Я хочу, чтобы этот подход был достаточно широким, чтобы быть применимым на других наборах, которые могут не иметь одинаковых ударов, поэтому я беру свободные плоские приставки.

проблема: Мой код ниже дает мне дрянные оценки. что не так ? thx

""" 
hypothesis: multimodal distrib sum of 3 gaussian distributions 

model description: 
* p1, p2, p3 are the probabilities for a point to belong to gaussian 1, 2 or 3 
==> p1, p2, p3 are the relative weights of the 3 gaussians 

* once a point is associated with a gaussian, 
it is distributed normally according to the parameters mu_i, sigma_i of the gaussian 
but instead of considering sigma, pymc prefers considering tau=1/sigma**2 

* thus, PyMc must guess 8 parameters: p1, p2, mu1, mu2, mu3, tau1, tau2, tau3 

* priors on p1, p2 are flat between 0.1 and 0.9 ==> 'pm.Uniform' variables 
with the constraint p2<=1-p1. p3 is deterministic ==1-p1-p2 

* the 'assignment' variable assigns each point to a gaussian, according to probabilities p1, p2, p3 

* priors on mu1, mu2, mu3 are flat between 40 and 120 ==> 'pm.Uniform' variables 

* priors on sigma1, sigma2, sigma3 are flat between 4 and 12 ==> 'pm.Uniform' variables 
""" 

    import numpy as np 
    import pymc as pm 

    data = np.loadtxt('distrib.txt') 
    Npts = len(data) 

    mumin = 40 
    mumax = 120 
    sigmamin=4 
    sigmamax=12 

    p1 = pm.Uniform("p1",0.1,0.9) 
    p2 = pm.Uniform("p2",0.1,1-p1) 
    p3 = 1-p1-p2 
    assignment = pm.Categorical('assignment',[p1,p2,p3],size=Npts) 
    mu = pm.Uniform('mu',[mumin,mumin,mumin],[mumax,mumax,mumax]) 
    sigma = pm.Uniform('sigma',[sigmamin,sigmamin,sigmamin], 
         [sigmamax,sigmamax,sigmamax]) 
    tau = 1/sigma**2 

    @pm.deterministic 
    def assign_mu(assi=assignment,mu=mu): 
     return mu[assi] 

    @pm.deterministic 
    def assign_tau(assi=assignment,sig=tau): 
     return sig[assi] 

    hypothesis = pm.Normal("obs", assign_mu, assign_tau, value=data, observed=True) 

    model = pm.Model([hypothesis, p1, p2, tau, mu]) 
    test = pm.MCMC(model) 
    test.sample(50000,burn=20000) # conservative values, let's take a coffee... 

    print('\nguess\n* p1, p2 = ', 
      np.mean(test.trace('p1')[:]),' ; ', 
      np.mean(test.trace('p2')[:]),' ==> p3 = ', 
      1-np.mean(test.trace('p1')[:])-np.mean(test.trace('p2')[:]), 
      '\n* mu = ', 
      np.mean(test.trace('mu')[:,0]),' ; ', 
      np.mean(test.trace('mu')[:,1]),' ; ', 
      np.mean(test.trace('mu')[:,2])) 

    print('why does this guess suck ???!!!')  

Я могу отправить файл данных 'distrib.txt'. Это ~ 500 кб, и данные приведены ниже. Например последний прогон дал мне:

p1, p2 = 0.366913192214 ; 0.583816452532 ==> p3 = 0.04927035525400003 
mu = 77.541619286 ; 75.3371615466 ; 77.2427165073 

в то время, очевидно, натыкается около ~ 55, ~ 75 и ~ 90, с вероятностью около ~ 0,2 ~ 0,5 и ~ 0,3

Answer 1

вас есть проблема, описанная здесь: Negative Binomial Mixture in PyMC

Проблема: категориальная переменная слишком медленно сходится, чтобы распределить три компонента, чтобы они были близки.

Во-первых, мы создаем ваши тестовые данные:

data1 = np.random.normal(55,5,2000) 
data2 = np.random.normal(75,5,5000) 
data3 = np.random.normal(90,5,3000) 
data=np.concatenate([data1, data2, data3]) 
np.savetxt("distrib.txt", data) 

Затем построим гистограмму, окрашенную уступке задней группы:

tablebyassignment = [data[np.nonzero(np.round(test.trace("assignment")[:].mean(axis=0)) == i)] for i in range(0,3) ] 
plt.hist(tablebyassingment, bins=30, stacked = True) 

Это в конечном счете сходятся, но не достаточно быстро быть полезным для вас.

Вы можете решить эту проблему путем подбора значения уступки перед началом MCMC:

from sklearn.cluster import KMeans 
kme = KMeans(3) 
kme.fit(np.atleast_2d(data).T) 
assignment = pm.Categorical('assignment',[p1,p2,p3],size=Npts, value=kme.labels_) 

Что дает: Использование K-средства для инициализации категорический не может работать все время, но это лучше, чем не сходится.