I have a data set where I am trying to hierarchically estimate the response distribution but running into indexing and array broadcast problems. It is a repeated measure design each participant does some number of trials in each condition and the parameters need to be multiplied by an additional variable before the likelihood is estimated.
I am new to pymc3 and bayesian modeling in general. Any help will be appreciated!!
What I tried is as follows:
sub_cond = np.concatenate((sub_idx,cond_idx),axis=1)
with pm.Model() as hierarchical_model:
# Hyperpriors
mu_a = pm.Normal('mu_a', mu=0., sigma=100,shape=n_cond)
sigma_a = pm.HalfNormal('sigma_a', 5.)
mu_b = pm.Normal('mu_b', mu=0., sigma=100,shape=n_cond)
sigma_b = pm.HalfNormal('sigma_b', 5.)
# mean for each subject, distributed around group mean mu_a
a = pm.Normal('a', mu=mu_a, sigma=sigma_a, shape=[n_sub,n_cond])
# beta for each each subject, distributed around group mean mu_b
b = pm.Normal('b', mu=mu_b, sigma=sigma_b, shape=[n_sub,n_cond])
transformed_a = a[sub_cond]*loc_idx
transformed_b = b[sub_cond]*loc_idx
# Data likelihood with sigma without random error
y = pm.TruncatedNormal( 'y',mu=transformed_a, sigma=transformed_b,observed=data)
hierarchical_trace = pm.sample(2000, tune=2000, target_accept=.9)
after_trace_burned = hierarchical_trace[500:]
pm.model_to_graphviz(hierarchical_model)