PS - I tried a much simpler thing and it seemed to work pretty well actually.
Just wrapped the mu=pm.Normal(… in a tt.sort() (was also talked about in the old discussion). That seems to have caused some improvement. Interestingly, if I also wrap the tau=pm.Gamma(… in the tt.sort(), then the performance is severely downgrade.
Regardless, the above issue is still at hand because the error occurred in the line for applying transform to mu line.