Thanks! This function works for me. Meanwhile, it leads to one other problem. I want to construct a matrix whose diagonal is the result of pm.math.erf(x) . For example, just like what we can do in numpy np.diag(scipy.special.erf(x)) . How can we implement this in pymc?