That makes sense in the analytic case, but you mean in the monte-carlo case? I presume just 4D Newton's method on the final integral over the pdf_X to get the value we're looking for?
It makes sense, I'm just surprised the blog never comes out and says "you have to do a brute force search to get the parameters for the likelihoods A ~ Beta(2,13) and B ~ Beta(3,11) for P(B > A) to have about the right value".