I think almost everything in the continuous case is good to go, except for the randomization probability.
Sometimes the P(theta_a < theta_b | Data) < 0, so I adjusted for that. I'm not sure why this is the case.
Fixed randomization probability.
As it turns out, it is actually easy to Calculate P(X + delta < Y | data) when X and Y are independent normals, because X - Y itself is normal! No need to do numerical integration for the continuous case: the function pnorm() is key.