python - Mezcla Binomial Negativa en PyMC

mcmc python

Question

Estoy tratando de encajar una mezcla binomial negativa con PyMC. Parece que hago algo mal, porque el predictivo no se parece en nada a los datos de entrada. El problema probablemente esté en el previo de los parámetros binomiales Negativos. ¿Alguna sugerencia?

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

Answer 1

Ha implementado correctamente una estimación bayesiana de una mezcla de tres distribuciones, pero el modelo MCMC arroja valores erróneos.

El problema es que la category no converge lo suficientemente rápido, y los parámetros en los means , alphas y dd se escapan de los valores buenos antes de que la category decida qué puntos pertenecen a qué distribución.

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

Puede ver que la category posterior para la category es incorrecta al visualizarla:

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)

La maximización de expectativas es el enfoque clásico para estabilizar las variables latentes, pero también puede usar los resultados del ajuste rápido y sucio de k-means para proporcionar valores iniciales para la MCMC:

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

Entonces las estimaciones convergen a valores de aspecto razonable.

Para su alfa anterior, puede usar la misma distribución para todos ellos:

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

Este problema no es específico de la distribución binomial negativa; Las mezclas de Dirichlet de distribuciones normales fallan de la misma manera; resulta de tener una distribución categórica de alta dimensión que MCMC no es eficiente para optimizar.