How to verify that uncertainty (estimated from pymc3) is accurate?

Thanks, good find.

Took a look through the code and there’s a bunch of obfuscation because of class objects but as far as I can tell it boils down to the exact equation I wrote above, just without the scaling to (0,1):

def order_stat(prior_val, posterior_samples):
        return np.sum(prior_val < posterior_samples, axis=0)