Thanks! Your insight about why scan is slow is helpful. Since my n is always 2, I may be able to modify the code for my use case… I’ll give it a try. I have read Saatchi’s thesis and was considering implementing it myself, but found that it was already implemented here 
and since I’m already using pymc for other models I’m comparing this to, it would be nice to continue using pymc.