Hi Maria,
I meant to get back to you on this - I actually started to play with the model, and see if there are ways to replace the theano.scan. I also come to the conclusion that it is possibly too difficult, and instead rewrite it slightly so that it doesnt output the whole 4D matrix (in the hope that it would make things faster). Unfortunately it is still pretty slow and the sample doesnt seems to work. Anyway, here is the code https://github.com/junpenglao/Planet_Sakaar_Data_Science/blob/master/PyMC3QnA/discourse_1735.ipynb, hopefully it can give you some inspiration. I suggest try just model 1 or 2 subjects and make sure it works first, as the sampling is too slow