Блог

Я пытаюсь совместить отрицательную биномиальную смесь с PyMC. Кажется, я делаю что-то не так, потому что прогноз совсем не похож на входные данные. Проблема, вероятно, в предшествующих отрицательных биномиальных параметрах. Есть предложения?

     from sklearn.cluster import KMeans
    import pymc as mc
    n = 3 #Number of components of the mixture
    ndata = len(data)

    dd = mc.Dirichlet('dd', theta=(1,)*n)
    category = mc.Categorical('category', p=dd, size=ndata)

    kme = KMeans(n) # This is not needed but it is to help convergence
    kme.fit(data[:,newaxis])
    alphas = mc.TruncatedNormal('alphas', kme.cluster_centers_[:,0], 0.1, a=0. ,b=100000 ,size=n)
    means = mc.TruncatedNormal('means', kme.cluster_centers_[:,0],0.1,a=0.0 ,b=100000, size=n)

    @mc.deterministic
    def mean(category=category, means=means):
        return means[category]

    @mc.deterministic
    def alpha(category=category, alphas=alphas):
        return alphas[category]

    obs = mc.NegativeBinomial('obs', mean, alpha, value=data, observed = True)

    predictive = mc.NegativeBinomial('predictive', mean, alpha)

    model = mc.Model({'dd': dd,
                  'category': category,
                  'alphas': alphas,
                  'means': means,
                  'predictive':predictive,
                  'obs': obs})

    mcmc = mc.MCMC( model )
    mcmc.sample( iter=n_samples, burn=int(n_samples*0.7))
  

Вы правильно реализовали байесовскую оценку смеси трех распределений, но модель MCMC дает неправильные значения.

Проблема в том, что category она недостаточно быстро сходится, и параметры в means , alphas , и dd убегают от хороших значений, прежде category чем решить, какие точки принадлежат какому распределению.

 data = np.atleast_2d(list(mc.rnegative_binomial(100., 10., size=s))  
    list(mc.rnegative_binomial(200., 1000., size=s))  
    list(mc.rnegative_binomial(300., 1000., size=s))).T
nsamples = 10000
  

Вы можете видеть, что последующее значение for category неверно, визуализируя его:

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

Максимизация ожиданий — это классический подход к стабилизации скрытых переменных, но вы также можете использовать результаты быстрой и грязной подгонки k-средних для получения начальных значений для MCMC:

 category = mc.Categorical('category', p=dd, size=ndata, value=kme.labels_)
  

Тогда оценки сходятся к разумным значениям.

Для вашего предыдущего варианта alpha вы можете просто использовать один и тот же дистрибутив для всех из них:

 alphas = mc.Gamma('alphas', alpha=1, beta=.0001 ,size=n)
  

Эта проблема не специфична для отрицательного биномиального распределения; Дирихле-смеси нормальных распределений терпят неудачу таким же образом; это связано с наличием многомерного категориального распределения, которое MCMC неэффективно при оптимизации.