Note: this topic is very similar to this thread (How do we predict on new unseen groups in a hierarchical model in PyMC3?) which I read and adopted to my particular model:
I have a hierarchical model with e.g. 2 levels of hierarchy: clusters and instances.
For instance, I define the model like this:
def model_factory(y, cluster, instance, n_clusters, n_instances): with pm.Model() as model: a = pm.Normal('a', mu= 0., sigma=3., shape=n_clusters) sigma_a = pm.Exponential('sigma_a', 1., shape=n_clusters) mu_tilde = pm.Normal('mu_t', mu=0., sigma=1., shape=n_instances) mu = pm.Deterministic('mu', mu_tilde * sigma_a[cluster] + a[cluster]) sigma = pm.Exponential('sigma', 1., shape=n_instances) data_obs = pm.Normal('data', mu=mu[instance], sigma=sigma[instance], observed=y) return model
You can see that I have a hyperprior for the clusters. In each clusters, there can be multiple instances which use the prior of the cluster.
For instance, I create toy data like this:
np.random.seed(123) N_clusters = 3 # Number of clusters N_instances = [10, 5, 1] # Number of instances per cluster total_instances = sum(N_instances) N = 100 # Number of samples per instance cluster_means = [1., 2., 3.] # Mean of means within cluster cluster_means_std = [0.3, 0.1, 0.1] # Std of means within cluster std = 0.5 data =  true_means =  for i in range(N_clusters): if N_instances[i] > 0: means = np.random.normal(loc=cluster_means[i], scale=cluster_means_std[i], size=N_instances[i]) true_means = np.append(true_means, means) data.append(np.array([np.random.normal(means[j], std, N) for j in range(N_instances[i])])) data = np.vstack(data).reshape(-1) clusters =  for i in range(N_clusters): clusters += [i] * N_samples[i] instance = np.repeat(np.arange(sum(N_samples)), N)
i.e., I have 3 clusters and each cluster has a certain number of instances (and each instances has 100 data points).
Now I am fitting the model to this data:
with model_factory(data, clusters, instance, N_clusters, total_instances): trace = pm.sample(tune=2000, target_accept=0.99) idata = az.from_pymc3(trace)
Now let’s assume I want to predict the data for an unseen instance of cluster 1. Following the thread mentioned above, I am trying to do this like this:
df = pm.trace_to_dataframe(trace, varnames=['a', 'sigma_a'], include_transformed=True) with model_factory(0, , , 3, 1): post_pred = pm.sample_posterior_predictive(df.to_dict('records'))
, i.e., I am only taking the parameters at the cluster level (
sigma_a) from the posterior trace and specify that my instance (0) belongs to cluster 1. However, this leads to the following error:
ValueError: 0-dimensional argument does not have enough dimensions for all core dimensions ('i_0_0',)
Note that if I change the number of instances to
total_instances in the call above, i.e.:
with model_factory(0, , , 3, total_instances): post_pred = pm.sample_posterior_predictive(df.to_dict('records'))
I am not getting any error, but the posterior prediction seems to ignore the posterior samples for
a, sigma_a and instead samples from the hyperprior (mean=0, sigma=3.) specified above. Moreover, it doesn’t make sense to me to specify the number of total instances here, because I want to create a new model with an arbitrary number of instances, which just uses the learned parameter for the clusters from the fit.