Fitting a hierarchical model with different sample size


I am trying to define a hierarchical model formulated as:
mu ~ N(50,10^2)
theta_i ~ N(mu, 10^2)
y_i ~ N(theta_i, 1)

Here let’s say we have 100 y_i and theta_i, but for each y_i is a vector but their dimension is not the same. From the tutorial I learned to define the parameters as:
with model:
mu = pm.Normal(‘mu’, 50, 10)
theta = pm.Normal(‘theta’, mu, 10, shape=100)

but I am not sure how to define the likelihood in this case. Could anyone provide an example? Thanks!

Concatenate all your y into one long 1-D array, and create a design matrix X such that the dot product of X and \theta can be used as your prior means on y.

Thank you so much!