How to model a multivariate multilevel linear model

Questions
2

This is a little bit cleaner. You can subsequently use a for loop to build up mu:

import theano
import theano.tensor as tt

with pm.Model() as tight2:
    beta_mu_hp = pm.Normal('beta_mu_hp', 0., 1., shape=7)
    beta_sd_hp = pm.HalfCauchy('beta_sd_hp', 1., shape=7)
    
    alpha = pm.Normal('alpha', mu=beta_mu_hp[0], sd=beta_sd_hp[0], shape=n_site)
    
    beta_mat = [pm.Normal('_beta_{}'.format(i), mu=beta_mu_hp[i+1], sd=beta_sd_hp[i+1], shape=n_site) 
                for i in range(6)]
    beta = pm.Deterministic('beta', tt.stack(beta_mat, axis=1))
    
    sigma = pm.HalfCauchy('sigma', 1)
    
    tr_tight = pm.sample_prior_predictive(100)

@junpenglao

I feel like this should work:

with pm.Model() as tight:
    beta_mu_hp = pm.Normal('beta_mu_hp', 0., 1., shape=7)
    beta_sd_hp = pm.HalfCauchy('beta_sd_hp', 1., shape=7)
    
    alpha = pm.Normal('alpha', mu=beta_mu_hp[0], sd=beta_sd_hp[0], shape=n_site)
    
    beta = pm.Normal('beta', mu=beta_mu_hp[1:].reshape((6,1)), sd=beta_sd_hp[1:].reshape((6,1)), shape=(6,n_site))
    sigma = pm.HalfCauchy('sigma', 1)
    
    tr_tight = pm.sample_prior_predictive(100)

but it generates an error during prior predictive sampling:

~/anaconda3/lib/python3.7/site-packages/pymc3-3.6-py3.7.egg/pymc3/distributions/distribution.py in <listcomp>(.0)
    670     elif broadcast_shape[:len(size_tup)] == size_tup:
    671         suffix = broadcast_shape[len(size_tup):] + dist_shape
--> 672         samples = [generator(*args, **kwargs).reshape(size_tup + (1,)) for _ in range(np.prod(suffix, dtype=int))]
    673         samples = np.hstack(samples).reshape(size_tup + suffix)
    674     else:

ValueError: cannot reshape array of size 600 into shape (100,1)