I am a bit confused by your model. In my opinion, you should not do
shifts = pm.Deterministic('shifts', tt.mean(batches-f.T, axis=1))
pm.Normal('infer', mu, sigma, observed = shifts, total_size= n_subjects)
I cant quite wrap my head around why you are taking the deterministic shifts as the observed… What are you intend to do here?
Also, are you sure you want to compute the mean over axis=1? If you are modelling the shift over n_obs you should compute the mean across axis=0