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.
Model
$g_i\sim\text{Categorical}(\theta^i_1,\ldots,\theta^i_g);~~~\vec{\theta^i}=(\theta^i_1,\ldots,\theta^i_g)\sim\text{Dirichlet}(0.5,\ldots,0.5)$
$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,
num_groups))
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]) : \
inv_logit(p)))
# 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.use_step_method(pymc.DiscreteMetropolis,
group[i],
proposal_distribution='Prior')
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.
However, with the Categorical part, the traces look like the following (showing only one plot for $\mu$ and $\sigma$ each; others look quite similar).
Can someone point out what I might be doing wrong here? Some of the potential issues I can see are:
- The acceptance rate for Categorical is very low (5-7%).
- 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.
- 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.
