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)