I am trying to implement a "hierarchical" model in PYMC in which the membership of observations to groups is not static (similar to the latent assignment of words to topics in Latent Dirichlet Allocation).

Problem Setting

I have a set of c coins. The probability of heads $p_i$ for coin $i$ is observed for all coins. The coins come from $g$ groups, where a group contains coins with similar $p_i$ values. The value of $g$ is known. However, the true assignment of coin $i$ to group $g_i$ is not known apriori, and I wish to fit a hierarchical model to my observations to identify the underlying groups and infer the distribution of coins in each group.



$p_i\sim\text{Logit-Normal}(\mu_{g_i}, \sigma_{g_i})$

$\mu_{g_i}\sim\text{Uniform}(-4, 4)$

$\sigma_{g_i}\sim\text{Uniform}(0.05, 0.3)$

I do realize that this might not qualify as a "fully Bayesian" hierarchical model, since I have not taken into consideration that uncertainty in the parameters of my priors for $\mu_{g_i}$ and $\sigma_{g_i}$. I am trying to keep things simple in my question. Moreover, going by the plots of Logit-normal for different values of $\mu_{g_i}$ and $\sigma_{g_i}$, I have a very strong belief that $\mu_{g_i}\in[-4,4]$ and $\sigma_{g_i}<0.3$. Hence, the choice of uniform priors. Moreover, as you would see below, the number of observations are pretty less, so I believe that the addition of hyperparameters would only worsen my situation.

I wish to infer the posterior distributions using the observed $p_i$ values. A crucial thing to note is that in general, $c$ will be of the order of $50-100$, which means I do not have many observations at my disposal.

I have implemented the model in PYMC. Here is the relevant code snippet:

num_groups = 5
num_coins = 50

def inv_logit(self, x):
    expo = numpy.exp(x)
    return expo/(1.0 + expo)

def logit(self, x):
    return numpy.log(x/(1.0 - x))

# Dirichlet prior for the parameters of Categorical
theta = pymc.Container([pymc.CompletedDirichlet('theta_%s' % i, 
                            pymc.Dirichlet('dir_%s' % i, 
                                theta = [0.5]*num_groups))
                        for i in xrange(0, num_coins)])

# Group assignment for a classifier is sampled from a Categorical
group = pymc.Container([pymc.Categorical('categ_%s' % i, 
                                        p = theta[i],
                                        value = numpy.random.randint(0,
                         for i in xrange(0, num_coins)])

# Uniform prior for the mu parameter of Logit-Normal
mu = pymc.Container([pymc.Uniform('mu_%s' % j, 
                                    lower = -4, 
                                    upper = 4, 
                                    value = 0)
                        for j in xrange(0, num_groups)])

# Uniform prior for the sigma parameter of Logit-Normal
sigma = pymc.Container([pymc.Uniform('sigma_%s' % j, 
                                    lower = 0.05, 
                                    upper = 0.3, 
                                    value = 0.15)
                        for j in xrange(0, num_groups)])

# Converting sigma to tau because of PYMC parameterization 
tau = pymc.Container([pymc.Lambda('tau_%s' % j, 
                                lambda s = sigma[j]: 1.0/(s**2))
                        for j in xrange(0, len(sigma))])

# p value for a classifier follows Logit-Normal distribution
obs = pymc.Container([pymc.Normal('obs_%s' % i,
                                mu = pymc.Lambda('omu_%s' % i, 
                                              lambda g=group[i]: mu[g]),
                                tau = pymc.Lambda('otau_%s' % i, 
                                              lambda g=group[i]: tau[g]),
                                value = logit(p),
                                observed = True)
                        for i, p in enumerate(observations)])

# Predictive distribution
# Note that if x ~ LogitNormal, then logit(x)~Normal
pred = pymc.Container(pymc.Lambda('predictive_%s' % j,
                lambda p = pymc.Normal('log_pred_%s' % j,
                                        mu = mu[j],
                                        tau = tau[j]) : \

# Adding Stochastics to model dictionary
for j in xrange(0, num_groups):
    model_dict['predictive_%s' % j] = pred[j]
    model_dict['mu_%s' % j] = mu[j]
    model_dict['sigma_%s' % j] = sigma[j]

for i in xrange(0, num_coins)
    model_dict['obs_%s' % i] = self.obs[i]

model = pymc.Model(model_dict)

# commented because MAP estimate with Categorical generates warning
# map_estimate = pymc.MAP(model)
# map_estimate.fit(method='fmin_l_bfgs_b')

# Using a special sampler and proposal for Categorical
mcmc = pymc.MCMC(model)
for i in xrange(0, num_coins):

mcmc.sample(iter = 50000,
            burn = 25000,
            thin = 25)

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

The use of a different step method for Categorical is motivated by this answer.

I am not getting even mildly decent results from this model. The Categorical part of the model seems to be the issue. If I remove the Categorical part from my model, and provide observations to the rest of the model with the correct coin-group assignment, the posterior distribution nicely fits the data. posterior predictive, true histogram without Categorical

However, with the Categorical part, the traces look like the following (showing only one plot for $\mu$ and $\sigma$ each; others look quite similar). posterior predictive, true histogram with Categorical mu, sigma group0, group1 group2,group3

Can someone point out what I might be doing wrong here? Some of the potential issues I can see are:

  1. The acceptance rate for Categorical is very low (5-7%).
  2. Number of parameters are too high, especially because of the number of $\theta$ parameters are $c*g$ . In fact, the number of parameters are always going to be several times larger than the number of observations, no matter how many observations I have.
  3. MCMC has not yet converged. Some of the categorical plots seem to hint at this, but I am not sure. With 50,000 iterations, sampling took about 700sec.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.